text:
abstract type DataBase end

"""
`dbPDBe` stores the residue `number` and `name` in PDBe as strings.
"""
@auto_hash_equals struct dbPDBe <: DataBase
    number::String # Cross referenced residue number
    name::String # Cross referenced residue name
end

"""
`dbInterPro` stores the residue `id`, `number`, `name` and `evidence` in InterPro as strings.
"""
@auto_hash_equals struct dbInterPro <: DataBase
    id::String
    number::String # Cross referenced residue number
    name::String # Cross referenced residue name
    evidence::String
end

"""
`dbEnsembl` stores the residue (gene) accession `id`, the `transcript`,
`translation` and `exon` ids in Ensembl as strings, together with the residue
`number` and `name` using the UniProt coordinates.
"""
@auto_hash_equals struct dbEnsembl <: DataBase
    id::String # (gene) accession id
    number::String # Cross referenced residue number
    name::String # Cross referenced residue name
    transcript::String
    translation::String
    exon::String
end

for ref_type in [:dbUniProt, :dbPfam, :dbNCBI]
    @eval begin

        @auto_hash_equals struct $(ref_type) <: DataBase
            id::String # The cross reference database identifier
            number::String # Cross referenced residue number
            name::String # Cross referenced residue name
        end

    end
end

@doc """
`dbUniProt` stores the residue `id`, `number` and `name` in UniProt as strings.
""" dbUniProt

@doc """
`dbPfam` stores the residue `id`, `number` and `name` in Pfam as strings.
""" dbPfam

@doc """
`dbNCBI` stores the residue `id`, `number` and `name` in NCBI as strings.
""" dbNCBI

for ref_type in [:dbPDB, :dbCATH, :dbSCOP, :dbSCOP2, :dbSCOP2B]
    @eval begin

        @auto_hash_equals struct $(ref_type) <: DataBase
            id::String
            number::String
            name::String
            chain::String
        end

    end
end

@doc """
`dbPDB` stores the residue `id`, `number`, `name` and `chain` in PDB as strings.
""" dbPDB

@doc """
`dbCATH` stores the residue `id`, `number`, `name` and `chain` in CATH as
strings.
""" dbCATH

@doc """
`dbSCOP` stores the residue `id`, `number`, `name` and `chain` in SCOP as
strings.
""" dbSCOP

@doc """
`dbSCOP2` stores the residue `id`, `number`, `name` and `chain` in SCOP2 as
strings.
""" dbSCOP2

@doc """
`dbSCOP2B` stores the residue `id`, `number`, `name` and `chain` in SCOP2B as
strings. *SCOP2B* is expansion of *SCOP2* domain annotations at superfamily
level to every *PDB* with same *UniProt* accession having at least 80% *SCOP2*
domain coverage.
""" dbSCOP2B

"""
Returns "" if the attributte is missing
"""
function _get_attribute(elem::LightXML.XMLElement, attr::String)
    text = attribute(elem, attr)
    if text === nothing || text == "None"
        return("")
    else
        return(text)
    end
end

"""
Returns `missing` if the attributte is missing
"""
function _get_nullable_attribute(elem::LightXML.XMLElement,
                                 attr::String)::Union{String,Missing}
    text = attribute(elem, attr)
    (text === nothing || text == "None") ? missing : text
end

for ref_type in [:dbPDB, :dbCATH, :dbSCOP, :dbSCOP2, :dbSCOP2B]
    @eval begin

        function $(ref_type)(map::LightXML.XMLElement)
            $(ref_type)(
                _get_attribute(map, "dbAccessionId"),
                _get_attribute(map, "dbResNum"),
                _get_attribute(map, "dbResName"),
                _get_attribute(map, "dbChainId")
            )
        end
    end
end

for ref_type in [:dbUniProt, :dbPfam, :dbNCBI]
    @eval begin

        function $(ref_type)(map::LightXML.XMLElement)
            $(ref_type)(
                _get_attribute(map, "dbAccessionId"),
                _get_attribute(map, "dbResNum"),
                _get_attribute(map, "dbResName")
            )
        end
    end
end

function dbEnsembl(map::LightXML.XMLElement)
    dbEnsembl(
        _get_attribute(map, "dbAccessionId"),
        _get_attribute(map, "dbResNum"),
        _get_attribute(map, "dbResName"),
        _get_attribute(map, "dbTranscriptId"),
        _get_attribute(map, "dbTranslationId"),
        _get_attribute(map, "dbExonId")
    )
end

function dbInterPro(map::LightXML.XMLElement)
    dbInterPro(
        _get_attribute(map, "dbAccessionId"),
        _get_attribute(map, "dbResNum"),
        _get_attribute(map, "dbResName"),
        _get_attribute(map, "dbEvidence")
    )
end

function dbPDBe(map::LightXML.XMLElement)
    dbPDBe(
        _get_attribute(map, "dbResNum"),
        _get_attribute(map, "dbResName")
      )
end

"""
A `SIFTSResidue` object stores the SIFTS residue level mapping for a residue. It has the
following fields that you can access at any moment for query purposes:

    - `PDBe` : A `dbPDBe` object, it's present in all the `SIFTSResidue`s.
    - `UniProt` : A `dbUniProt` object or `missing`.
    - `Pfam` : A `dbPfam` object or `missing`.
    - `NCBI` : A `dbNCBI` object or `missing`.
    - `InterPro` : An array of `dbInterPro` objects.
    - `PDB` : A `dbPDB` object or `missing`.
    - `SCOP` : A `dbSCOP` object or `missing`.
    - `SCOP2` : An array of `dbSCOP2` objects.
    - `SCOP2B` : A `dbSCOP2B` object or `missing`.
    - `CATH` : A `dbCATH` object or `missing`.
    - `Ensembl` : An array of `dbEnsembl` objects.
    - `missing` : It's `true` if the residue is missing, i.e. not observed, in the structure.
    - `sscode` : A string with the secondary structure code of the residue.
    - `ssname` : A string with the secondary structure name of the residue.
"""
@auto_hash_equals struct SIFTSResidue
    PDBe::dbPDBe
    # crossRefDb
    UniProt::Union{dbUniProt,Missing}
    Pfam::Union{dbPfam,Missing}
    NCBI::Union{dbNCBI,Missing}
    InterPro::Array{dbInterPro,1}
    PDB::Union{dbPDB,Missing}
    SCOP::Union{dbSCOP,Missing}
    SCOP2::Array{dbSCOP2,1}
    SCOP2B::Union{dbSCOP2B,Missing}
    CATH::Union{dbCATH,Missing}
    Ensembl::Array{dbEnsembl,1}
    # residueDetail
    missing::Bool  # XML: <residueDetail dbSource="PDBe" property="Annotation" ...
    sscode::String # XML: <residueDetail dbSource="PDBe" property="codeSecondaryStructure"...
    ssname::String # XML: <residueDetail dbSource="PDBe" property="nameSecondaryStructure"...
end

# Getters
# -------

@inline _name(::Type{dbPDBe})    = "PDBe"
@inline _name(::Type{dbUniProt}) = "UniProt"
@inline _name(::Type{dbPfam})    = "Pfam"
@inline _name(::Type{dbNCBI})    = "NCBI"
@inline _name(::Type{dbInterPro}) = "InterPro"
@inline _name(::Type{dbPDB})     = "PDB"
@inline _name(::Type{dbSCOP})    = "SCOP"
@inline _name(::Type{dbSCOP2})   = "SCOP2"
@inline _name(::Type{dbSCOP2B})  = "SCOP2B"
@inline _name(::Type{dbCATH})    = "CATH"
@inline _name(::Type{dbEnsembl}) = "Ensembl"

@inline Base.get(res::SIFTSResidue, db::Type{dbPDBe})    = res.PDBe
@inline Base.get(res::SIFTSResidue, db::Type{dbUniProt}) = res.UniProt
@inline Base.get(res::SIFTSResidue, db::Type{dbPfam})    = res.Pfam
@inline Base.get(res::SIFTSResidue, db::Type{dbNCBI})    = res.NCBI
@inline Base.get(res::SIFTSResidue, db::Type{dbInterPro}) = res.InterPro
@inline Base.get(res::SIFTSResidue, db::Type{dbPDB})     = res.PDB
@inline Base.get(res::SIFTSResidue, db::Type{dbSCOP})    = res.SCOP
@inline Base.get(res::SIFTSResidue, db::Type{dbSCOP2})   = res.SCOP2
@inline Base.get(res::SIFTSResidue, db::Type{dbSCOP2B})  = res.SCOP2B
@inline Base.get(res::SIFTSResidue, db::Type{dbCATH})    = res.CATH
@inline Base.get(res::SIFTSResidue, db::Type{dbEnsembl}) = res.Ensembl

function Base.get(res::SIFTSResidue,
                  db::Type{T},
                  field::Symbol,
                  default::Union{String,Missing}=missing) where T <: Union{dbUniProt,dbPfam,dbNCBI,dbPDB,dbSCOP,dbSCOP2B,dbCATH}
    database = get(res, db)
    ismissing(database) ? default : getfield(database, field)
end

# Print
# -----

function Base.show(io::IO, res::SIFTSResidue)
    if res.missing
        println(io, "SIFTSResidue (missing)")
    else
        println(io, "SIFTSResidue with secondary structure code (sscode): \"",
                res.sscode, "\" and name (ssname): \"", res.ssname, "\"")
    end
    println(io, "  PDBe:")
    println(io, "    number: ", res.PDBe.number)
    println(io, "    name: ", res.PDBe.name)
    for dbname in [:UniProt, :Pfam, :NCBI, :PDB, :SCOP, :SCOP2B, :CATH]
        dbfield = getfield(res, dbname)
        if !ismissing(dbfield)
            println(io, "  ", dbname, " :")
            for f in fieldnames(typeof(dbfield))
                println(io, "    ", f, ": ",  getfield(dbfield, f))
            end
        end
    end
    length(res.SCOP2) > 0 && println(io, "  SCOP2: ", res.SCOP2)
    length(res.InterPro) > 0 && println(io, "  InterPro: ",  res.InterPro)
    length(res.Ensembl) > 0 && println(io, "  Ensembl: ",  res.Ensembl)
end

# Creation
# --------

function SIFTSResidue(residue::LightXML.XMLElement, missing_residue::Bool,
                      sscode::String, ssname::String)
    PDBe = dbPDBe(residue)
    UniProt = missing
    Pfam = missing
    NCBI = missing
    InterPro = dbInterPro[]
    PDB = missing
    SCOP = missing
    SCOP2 = dbSCOP2[]
    SCOP2B = missing
    CATH = missing
    Ensembl = dbEnsembl[]
    for crossref in get_elements_by_tagname(residue, "crossRefDb")
        db = attribute(crossref, "dbSource")
        if db == "UniProt"
            UniProt = dbUniProt(crossref)
        elseif db == "Pfam"
            Pfam = dbPfam(crossref)
        elseif db == "NCBI"
            NCBI = dbNCBI(crossref)
        elseif db == "InterPro"
            push!(InterPro, dbInterPro(crossref))
        elseif db == "PDB"
            PDB = dbPDB(crossref)
        elseif db == "SCOP"
            SCOP = dbSCOP(crossref)
        elseif db == "SCOP2"
            push!(SCOP2, dbSCOP2(crossref))
        elseif db == "SCOP2B"
            SCOP2B = dbSCOP2B(crossref)
        elseif db == "CATH"
            CATH = dbCATH(crossref)
        elseif db == "Ensembl"
            push!(Ensembl, dbEnsembl(crossref))
        else
            @warn(string(db, " is not in the MIToS' DataBases."))
        end
    end
    SIFTSResidue(PDBe,
                 UniProt,
                 Pfam,
                 NCBI,
                 InterPro,
                 PDB,
                 SCOP,
                 SCOP2,
                 SCOP2B,
                 CATH,
                 Ensembl,
                 missing_residue,
                 sscode,
                 ssname)
end

function SIFTSResidue(residue::LightXML.XMLElement)
    SIFTSResidue(residue, _get_details(residue)...)
end

# Mapping Functions
# =================

_is_All(::Any) = false
_is_All(::Type{All}) = true

"""
Parses a SIFTS XML file and returns a `OrderedDict` between residue numbers of
two `DataBase`s with the given identifiers. A `chain` could be specified
(`All` by default). If `missings` is `true` (default) all the residues are
used, even if they haven’t coordinates in the PDB file.
"""
function siftsmapping(filename::String,
                      db_from::Type{F},
                      id_from::String,
                      db_to::Type{T},
                      id_to::String;
                      chain::Union{Type{All},String}=All,
                      missings::Bool=true) where {F,T}
    mapping = OrderedDict{String,String}()
    xdoc = parse_file(filename)
    try
        for entity in _get_entities(xdoc)
            segments = _get_segments(entity)
            for segment in segments
                residues = _get_residues(segment)
                for residue in residues
                    in_chain = _is_All(chain)
                    key_data = _name(db_from) == "PDBe" ? attribute(residue, "dbResNum") : missing
                    value_data = _name(db_to) == "PDBe" ? attribute(residue, "dbResNum") : missing
                    if missings || !_is_missing(residue)
                        crossref = get_elements_by_tagname(residue, "crossRefDb")
                        for ref in crossref
                            source = attribute(ref, "dbSource")
                            if source == _name(db_from) && attribute(ref, "dbAccessionId") == id_from
                                key_data = _get_nullable_attribute(ref, "dbResNum")
                            end
                            if source == _name(db_to) && attribute(ref, "dbAccessionId") == id_to
                                value_data = _get_nullable_attribute(ref, "dbResNum")
                            end
                            if !in_chain && source == "PDB" # XML: <crossRefDb dbSource="PDB" ... dbChainId="E"/>
                                in_chain = attribute(ref, "dbChainId") == chain
                            end
                        end
                        if !ismissing(key_data) && !ismissing(value_data) && in_chain
                            key = key_data
                            if haskey(mapping, key)
                                @warn string("$key is already in the mapping with the value ",
                                             mapping[key],". The value is replaced by ",
                                             value_data)
                            end
                            mapping[key] = value_data
                        end
                    end
                end
            end
        end
    finally
        free(xdoc)
    end
    sizehint!(mapping, length(mapping))
end

"""
`parse(document::LightXML.XMLDocument, ::Type{SIFTSXML}; chain=All, missings::Bool=true)`

Returns a `Vector{SIFTSResidue}` parsed from a `SIFTSXML` file.
By default, parses all the `chain`s and includes missing residues.
"""
function Base.parse(document::LightXML.XMLDocument, ::Type{SIFTSXML};
                    chain::Union{Type{All},String}=All, missings::Bool=true)
    vector = SIFTSResidue[]
    for entity in _get_entities(document)
        for segment in _get_segments(entity)
            residues = _get_residues(segment)
            for residue in residues
                missing_residue, sscode, ssname = _get_details(residue)
                if missings || !missing_residue
                    sifts_res = SIFTSResidue(residue, missing_residue, sscode, ssname)
                    if _is_All(chain) || (!ismissing(sifts_res.PDB) && sifts_res.PDB.chain == chain)
                        push!(vector, sifts_res)
                    end
                end
            end
        end
    end
    vector
end

# Find SIFTSResidue
# -----------------

for F in (:findall, :filter!, :filter)
    @eval begin
        function Base.$(F)(f::Function, list::AbstractVector{SIFTSResidue}, db::Type{T}) where T <: DataBase
            $(F)(list) do res
                database = get(res, db)
                if !ismissing(database)
                    f(database)
                end
            end
        end
    end
end


text:
@show "ha" ^ 5

# The ^ operator is really just call to the `repeat` function
@show repeat("ha", 5)

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type AbstractThresholdDither <: AbstractDither end

"""
    WhiteNoiseThreshold()

Use white noise as a threshold map.
"""
struct WhiteNoiseThreshold <: AbstractThresholdDither end

function binarydither!(::WhiteNoiseThreshold, out::GenericGrayImage, img::GenericGrayImage)
    tmap = rand(eltype(img), size(img))
    return out .= img .> tmap
end

"""
    ConstantThreshold(threshold)

Use a constant threshold map. Defaults to 0.5 if `threshold` isn't specified.
"""
struct ConstantThreshold{T<:Real} <: AbstractThresholdDither
    threshold::T

    function ConstantThreshold(; threshold=0.5)
        if threshold < 0 || threshold > 1
            throw(DomainError(threshold, "Threshold needs to be between 0 and 1."))
        end
        return new{typeof(threshold)}(threshold)
    end
end

function binarydither!(alg::ConstantThreshold, out::GenericGrayImage, img::GenericGrayImage)
    threshold = eltype(img)(alg.threshold)
    return out .= img .> threshold
end


text:
"""
Calculates the mean quality of the items in the given registry
"""
function benchmarkRegistryMean(registry, agents)
    mean(registry)
end

"""
Calculates the mean quality of the top 10 items in the given registry
"""
function top10Mean(history)
    sorted = mean(sort(history, rev=true)[1:10])
end

"""
Calculates the mean quality of the top X(count) items in the given registry
"""
function topMean(history, count)
    sorted = mean(sort(history, rev=true)[1:count])
end

"""
Calculates the mean accurancy of a group of agents
"""
function meanAccuracy(registry, agents)
    accs = []
    for agent in agents
        push!(accs, agent.accuracy)
    end
    meanAccuracy = mean(accs)
    return meanAccuracy
end

"""
Calculates the number of agents possesing tokens thus being able to stake
them to participate in challenges and votings
"""
function activeAgents(registry, agents)
    length(filter(a -> a.balance > 0, agents))
end

"""
Calculates the difference in accurancy of agents who posses tokens and those
who does not hold any tokens
"""
function benchmarkAccuracyBoost(registry, agents)
    eAcc = []
    count = 0
    for agent in agents
        if agent.balance > 0
            push!(eAcc, agent.accuracy)
            count = count + 1
        end
    end
    effective = mean(eAcc)
    normal = meanAccuracy(registry, agents)
    boost = effective-normal
end

"""
Calculates the difference in accurancy of agents weighted by the number of tokens
they posses
"""
function benchmarkWeightedAccuracyBoost(registry, agents)
    count = 0
    eAcc = 0;
    for agent in agents
        if agent.balance > 0
            eAcc += agent.balance * agent.accuracy
            count += agent.balance
        end
    end
    effective = eAcc/count
    normal = meanAccuracy(registry, agents)
    boost = effective-normal
end

"""
Calculates the correlation between agents accuracy and balance
"""
function benchmarkAccuracyBalanceCorrelation(registry, agents)
    a = [agent.accuracy for agent in agents]
    b = [agent.balance for agent in agents]
    corspearman(a, b)
end

"""
Calculates the gini coefficient - how equally tokens are distributed among agents
"""
function gini(registry, agents)
    x = [agent.balance for agent in agents]
    n = length(x)
    xx = sort(x)
    2*(sum(collect(1:n).*xx))/(n*sum(xx))-1
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: condenser.mso 555 2008-07-18 19:01:13Z rafael $
#*--------------------------------------------------------------------
type condenserSteady
	condenserSteady()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop in the condenser",
				:Default=>0,
				:Symbol=>"\\Delta _P"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour inlet stream",
				:PosX=>0.16,
				:PosY=>0,
				:Symbol=>"_{in}^{Vapour}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.53,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat Duty",
				:PosX=>1,
				:PosY=>0.08,
				:Symbol=>"Q_{in}",
				:Protected=>true
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Bubble Temperature",
				:Protected=>true,
				:Symbol=>"T_{bubble}"
			)),
			temp_delta(Dict{Symbol,Any}(
				:Brief=>"Degrees subcooled",
				:Symbol=>"\\Delta T_{subcooled}"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.50,
				:PosY=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Pressure  Indicator of Condenser",
				:Protected=>true,
				:PosX=>0.32,
				:PosY=>0
			)),
			[
				:(InletVapour.F = OutletLiquid.F),
				:(InletVapour.z = OutletLiquid.z),
				:(InletVapour.F*InletVapour.h + InletQ = OutletLiquid.F*OutletLiquid.h),
				:(OutletLiquid.P = InletVapour.P - Pdrop),
				:(Tbubble = PP.BubbleT(OutletLiquid.P,OutletLiquid.z)),
				:(OutletLiquid.T = Tbubble-Deg_Subcooled),
				:(TI * "K" = OutletLiquid.T),
				:(PI * "atm" = OutletLiquid.P),
			],
			[
				"Molar Flow Balance","Molar Composition Balance","Energy Balance","Pressure Drop","Bubble Temperature","Temperature","Temperature indicator","Pressure indicator",
			],
			[:PP,:NComp,:Pdrop,],
			[:InletVapour,:OutletLiquid,:InletQ,:Tbubble,:Deg_Subcooled,:TI,:PI,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Pdrop::press_delta
	InletVapour::stream
	OutletLiquid::liquid_stream
	InletQ::power
	Tbubble::temperature
	Deg_Subcooled::temp_delta
	TI::control_signal
	PI::control_signal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export condenserSteady
function setEquationFlow(in::condenserSteady)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
end
function atributes(in::condenserSteady,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/CondenserSteady"
	fields[:Brief]="Model of a  Steady State condenser with no thermodynamics equilibrium."
	fields[:Info]="== ASSUMPTIONS ==
* perfect mixing of both phases;
* no thermodynamics equilibrium.

== SET ==
* the pressure drop in the condenser;

== SPECIFY ==
* the InletVapour stream;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the condenser model has two control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
"
	drive!(fields,_)
	return fields
end
condenserSteady(_::Dict{Symbol,Any})=begin
	newModel=condenserSteady()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

module GeoStats

using Reexport

# reexport project modules
@reexport using GeoStatsBase
@reexport using Variography
@reexport using KrigingEstimators

end


text:
println("loading packages... ")

using ArgParse
using OrthographicProjection

println("packages OK")

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table! s begin
	"source"
		help = "A text file with Potree directory list"
		arg_type = String
		required = true
	"--projectname", "-p"
		help = "Project name"
		arg_type = String
		required = true
	"--output", "-o"
		help = "Output folder"
		arg_type = String
		required = true
	"--bbin"
        help = "Bounding box as 'x_min y_min z_min x_max y_max z_max' or Potree JSON volume model"
		required = true
	"--listpoints"
		help = "A text file containing a list of two points."
		arg_type = String
		required = true
	"--axis"
		help = "A vector in plane"
		arg_type = String
		default = "0 0 1"
	"--thickness"
		help = "Section thickness"
		arg_type = Float64
		default = 0.1
	end

	return parse_args(s)
end

function main()
	args = parse_commandline()

	OrthographicProjection.flushprintln("== params ==")
	for (arg,val) in args
		OrthographicProjection.flushprintln("$arg  =>  $val")
	end

	txtpotreedirs = args["source"]
	project_name = args["projectname"]
	output_folder = args["output"]
	bbin = args["bbin"]
	file_listpoints = args["listpoints"]
	axis_ = args["axis"]
	thickness = args["thickness"]

	b = tryparse.(Float64,split(bbin, " "))
	if length(b) == 6
		bbin = OrthographicProjection.AABB(b[4],b[1],b[5],b[2],b[6],b[3])
	end

	axis_y = tryparse.(Float64,split(axis_, " "))
	@assert length(axis_y) == 3 "a 3D axis needed"

	proj_folder, models = OrthographicProjection.preprocess(project_name, output_folder, bbin, file_listpoints, axis_y, thickness)
	OrthographicProjection.extract_models(txtpotreedirs, project_name, proj_folder, bbin, models)

end

@time main()


text:
using ProgressMeter: @showprogress

test = "389125467"
puzzle = "962713854"

parse_input(input::AbstractString) = parse.(Int, split(input, ""))

mutable struct Ring
    linked::AbstractVector
    current::Integer

    function Ring(v::AbstractVector)::Ring
        linked = similar(v)
        linked[v] .= [v[2:end]; v[1]]
        return new(linked, v[1])
    end
end

Base.getindex(ring::Ring, i::Integer) = ring.linked[i]
Base.setindex!(ring::Ring, v, i::Integer) = ring.linked[i] = v

function Base.convert(::Type{<: Vector},
                      ring::Ring;
                      from = ring.current
                     )::Vector
    v = similar(ring.linked)
    foldl(keys(v); init = from) do i, j
        v[j] = i
        return ring.linked[i]
    end
    return v
end

function move!(ring::Ring)::Ring
    pick = reduce(1:3; init = ([], ring.current)) do (pick, prev), _
        return let next = ring[prev]; ([pick; next], next) end
    end |> first
    dest = setdiff(predecessors(ring), pick) |> first

    ring[ring.current] = ring[last(pick)]
    ring[last(pick)] = ring[dest]
    ring[dest] = first(pick)
    next!(ring)
    return ring
end

predecessors(ring::Ring) = rem.(ring.current - 2:-1:ring.current - 5,
                                length(ring.linked),
                                RoundDown) .+ 1

next!(ring::Ring)::Ring = (ring.current = ring[ring.current]; ring)

function play(v::Vector, rounds::Integer)::Vector
    ring = Ring(v)
    @showprogress for _ in 1:rounds
        move!(ring)
    end
    return convert(Vector, ring; from = 1)
end

function q1(v::Vector, rounds::Integer = 100)::Int
    v = play(v, rounds)
    return parse(Int, join(v[2:end]))
end

function q2(v::Vector)::Int
    v = [v; range(maximum(v) + 1; step = 1, length = 1_000_000 - length(v))]
    v = play(v, 10_000_000)
    return v[2] * v[3]
end

let v = parse_input(test)
    @assert q1(v, 10) == 92658374
    @assert q1(v) == 67384529
    @assert q2(v) == 149245887792
end

v = parse_input(puzzle)
println("Q1: ", q1(v))
println("Q2: ", q2(v))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
module FrostHelperBlueBoosterModule

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/BlueBooster" FrostHelperBlueBoosterPlacement(x::Integer, y::Integer, boostTime::Number=0.25, respawnTime::Number=1.0, particleColor::String="LightSkyBlue", directory="objects/FrostHelper/blueBooster/", reappearSfx::String="event:/game/04_cliffside/greenbooster_reappear", enterSfx::String="event:/game/04_cliffside/greenbooster_enter", boostSfx::String="event:/game/04_cliffside/greenbooster_dash", releaseSfx::String="event:/game/04_cliffside/greenbooster_end", red::Bool=false, preserveSpeed::Bool=false)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

const placements = Ahorn.PlacementDict(
	"Blue Booster (Frost Helper)" => Ahorn.EntityPlacement(
		FrostHelperBlueBoosterPlacement
	)
)

Ahorn.editingOptions(entity::FrostHelperBlueBoosterPlacement) = Dict{String, Any}(
    "particleColor" => colors
)

sprite = "objects/booster/booster00"

function Ahorn.selection(entity::FrostHelperBlueBoosterPlacement)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FrostHelperBlueBoosterPlacement, room::Maple.Room)
	texture = strip(get(entity.data, "directory", "objects/FrostHelper/blueBooster/"), '/')
	realTexture = "$texture/booster00.png"

	Ahorn.drawSprite(ctx, realTexture, 0, 0)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using FFTW
using Plots

import Plots: default
default(show=false, reuse=true)

ENV["GKSwstype"]="nul"

function signal_plot(i, signal1, signal2, out; filebase = "out",
                     out_name = "convolution")
    println(i)
    n = length(signal1)
    shift_signal = real([signal2[i:n-1];signal2[1:i]])
    conv_sum = real(signal1).*shift_signal
    plt = plot(real(signal1); label = "signal")
    plt = plot!(plt, shift_signal; label = "filter")
    plt = plot!(plt, zeros(n); label = "area", ribbon=(zeros(n),conv_sum))
    plt = plot!(plt, real(out); label = out_name)
    savefig(filebase*lpad(string(i),5,string(0))*".png")
end

function is_norm(signal)
    if isapprox(sum(signal), 1.0)
        return true
    else
        return false
    end
end

function norm(signal)
    return signal ./ sum(signal)
end

# TODO: better name than temp_signal
function expectation_value(signal)
    temp_signal = signal
    if !is_norm(signal)
        temp_signal = norm(signal)
    end
    expectation_value = 0
    for i = 1:length(signal)
        expectation_value += i*temp_signal[i]
    end
    return expectation_value
end

function covariance(signal1, signal2)
    return expectation_value((signal1 .- expectation_value(signal1)).*
                             (signal2 .- expectation_value(signal2)))
end

function conv_lin(signal1::Array{Complex{Float64},1},
                  signal2::Array{Complex{Float64},1})
    n = length(signal1)
    out = Array{Complex{Float64},1}(undef,n)
    rsum = 0

    # time domain
    for i = 1:n-1
        # inner domain between signals
        for j = 1:length(signal2)
            if (j < length(signal1))
                index = i-j
                if i-j < 1
                    index = length(signal2)+(i-j)
                end
                rsum += signal1[j] * signal2[(index)%(length(signal2))+1]
            end
        end
        out[i] = rsum
        rsum = 0
    end

    return out
end

function conv_fft(signal1::Array{Complex{Float64},1},
                  signal2::Array{Complex{Float64},1})
    return ifft(fft(signal1).*fft(signal2))
end

function conv_plot(signal1::Array{Complex{Float64},1},
                   signal2::Array{Complex{Float64},1}; norm_factor = 1)
    n = length(signal1)
    out = zeros(Complex, n)
    rsum = 0

    # time domain
    for i = 1:n-1
        # inner domain between signals
        for j = 1:length(signal2)
            if (j < length(signal1))
                index = j-i
                if j-i < 1
                    index = length(signal2)+(j-i)
                end
                rsum += signal1[j] * signal2[(index)%(length(signal2))+1]
            end
        end

        out[i] = rsum / norm_factor
        rsum = 0

        signal_plot(i, signal1, signal2, out)

    end

    return out
end

function corr_lin(signal1::Array{Complex{Float64},1},
                  signal2::Array{Complex{Float64},1})
    n = length(signal1)
    out = Array{Complex{Float64},1}(undef,n)
    rsum = 0

    for i = 1:n-1
        for j = 1:length(signal2)
            if j < length(signal1)
                rsum += conj(signal1[j]) * signal2[(j+i)%(length(signal2))+1]
            end
        end
        out[i] = rsum
        rsum = 0
    end

    return out
end

function corr_plot(signal1::Array{Complex{Float64},1},
                   signal2::Array{Complex{Float64},1}; norm_factor = 1)
    n = length(signal1)
    out = zeros(Complex, n)
    rsum = 0

    for i = 1:n-1
        for j = 1:length(signal2)
            if j < length(signal1)
                rsum += conj(signal1[j]) * signal2[(j+i)%(length(signal2))+1]
            end

        end

        out[i] = rsum / norm_factor
        rsum = 0

        signal_plot(i, signal1, signal2, out)

    end

    return out
end

function main()
    x = [exp(-((i-50)/100)^2/.01) + 0im for i = 1:100]
    out_conv = conv_plot(x,x, norm_factor=20)

#=
    x = zeros(Complex{Float64},100)
    x[40:60] .= 1

    y = copy(x)
    y[40:60] .= [float(1-i*0.05) for i = 0:20]

    # Make plots
    out_corr = corr_plot(x,y;norm_factor=20; filebase = "corr_out")
    out_conv = conv_plot(x,y, norm_factor=20)
=#
end


text:
using TrajectoryOptimization
using Plots
# Pendulum
n = 2
m = 1
model, = Dynamics.pendulum # same as UDP pendulum dynamics

Qf = 30.0*Matrix(I,n,n)
Q = 0.3*Matrix(I,n,n)
R = 0.3*Matrix(I,m,m)
dt = 0.1
tf = 5.0
x0 = [0.0; 0.0]
xf = [pi; 0.0]

obj = LQRObjective(Q,R,Qf,tf,x0,xf)

opts = SolverOptions()
opts.cost_tolerance = 1e-5
opts.gradient_norm_tolerance = 1e-5
opts.verbose = true
opts.cache = true
opts.c2 = 100.0

solver = TrajectoryOptimization.Solver(model,obj,dt=dt,opts=opts)
U = zeros(solver.model.m, solver.N)
results, stats = TrajectoryOptimization.solve(solver,U)

# Cartpole
n = 4
m = 1
model = Model(Dynamics.cartpole_dynamic_udp!,n,m)

dt = .1; # time step
# quadratic costs
Qf = 1000*Matrix(I,n,n);
Q = .1*Matrix(I,n,n)
R = .01*Matrix(I,m,m);

tf = 5.0
x0 = [0.0; 0.0; 0.0; 0.0]
xf = [0.0; pi; 0.0; 0.0]

obj = LQRObjective(Q,R,Qf,tf,x0,xf)

opts = SolverOptions()
opts.cost_tolerance = 1e-5
opts.gradient_norm_tolerance = 1e-5
opts.verbose = true
opts.cache = true

solver = TrajectoryOptimization.Solver(model,obj,dt=dt,opts=opts)
U = zeros(solver.model.m, solver.N)
results, stats = TrajectoryOptimization.solve(solver,U)
println(results.X[:,end])
# plot(results.cost[1:results.termination_index])


text:
#!/usr/bin/env julia6
using ReadGlobal, InplaceRealFFT, Derivatives

function tplus!(r::Array{T,N},u::Array{T,N},v::Array{T,N}) where {T,N}
    Threads.@threads for i::Int in 1:length(u)
      @inbounds r[i] = 0.5*(u[i] + v[i])
    end
    return r
end

function tminus!(r::Array{T,N},u::Array{T,N},v::Array{T,N}) where {T,N}
    Threads.@threads for i::Int in 1:length(u)
      @inbounds r[i] = 0.5*(u[i] - v[i])
    end
    return r
end

function sym!(r,u::A,v::A) where {A<:PaddedArray{T,N} where {T,N}}
    urr = data(u)
    vrr = data(v)
    tplus!(data(r),urr,vrr)
    return r
end

function antisym!(r,u::A,v::A) where {A<:PaddedArray{T,N} where {T,N}}
    urr = data(u)
    vrr = data(v)
    tminus!(data(r),urr,vrr)
    return r
end

function MakeD()

  N,Fil = getnfilter()
  nx, ny, nz, lx, ly, lz = getdimsize()
  dim = (nx,ny,nz)
  isfile("fftw_wisdom") && FFTW.import_wisdom("fftw_wisdom")

  aux1 = PaddedArray("$(Fil)u1_N$N",(nx,ny,nz),false)
  write("$(Fil)D11_N$N",dx!(aux1,lx))
  read!("$(Fil)u1_N$N",aux1,false)

  aux2 = PaddedArray("$(Fil)u2_N$N",(nx,ny,nz),false)
  write("$(Fil)D22_N$N",dy!(aux2,ly))
  read!("$(Fil)u2_N$N",aux2,false)

  dy!(aux1,ly)
  dx!(aux2,lx)

  aux3 = similar(aux1)
  write("$(Fil)D12_N$N", real(sym!(aux3,aux2,aux1)))
  write("$(Fil)W12_N$N", real(antisym!(aux3,aux2,aux1)))

  read!("$(Fil)u3_N$N",aux2,false)
  write("$(Fil)D33_N$N",dz!(aux2,lz))
  read!("$(Fil)u3_N$N",aux2,false)

  read!("$(Fil)u1_N$N",aux1,false)

  dz!(aux1,lz)
  dx!(aux2,lx)
  write("$(Fil)D13_N$N", real(sym!(aux3,aux2,aux1)))
  write("$(Fil)W13_N$N", real(antisym!(aux3,aux2,aux1)))

  read!("$(Fil)u2_N$N",aux1,false)
  read!("$(Fil)u3_N$N",aux2,false)
  dz!(aux1,lz)
  dy!(aux2,ly)
  write("$(Fil)D23_N$N", real(sym!(aux3,aux2,aux1)))
  write("$(Fil)W23_N$N", real(antisym!(aux3,aux2,aux1)))

  return 0
end

FFTW.set_num_threads(Threads.nthreads())
MakeD()


text:
using EpistemicNetworkAnalysis
using UMAP
using Plots
using LinearAlgebra
using DataFrames
using Random
using Statistics
using CSV
using Distances
using Dates
using Colors
using GLM
using Clustering
using HypothesisTests

cd(Base.source_dir())
include("helpers.jl")

let # create a non-global scope

# Data
data = DataFrame(CSV.File("data/data.csv", normalizenames=true))

# Config
codes = [
    :NegativeVE,
    :PositiveVE,
    :Normalcy,
    :Safety,
    :SharedSpace,
    :ComingOfAge,
    :Friends,
    :Learning,
    :Mental,
    :Technology,
    :Policy,
    :Management
]

conversations = [:ThreadNumber]
units = [:ThreadNumber, :CommentPosition]
dropEmpty=true
sphereNormalize=true

# K-Means clustering approach
data[!, :LABEL] = repeat(["No Label"], nrow(data))
colorMap = Dict("No Label" => colorant"black")
ena = ENAModel(data, codes, conversations, units, dropEmpty=dropEmpty, sphereNormalize=sphereNormalize)
function autocluster!(ena, data, colorMap, k, dims)

    ### Prepare data
    X = Matrix{Float64}(ena.accumModel[!, dims])
    X = Matrix{Float64}(transpose(X))

    ### Using k-means to detect labels
    results = kmeans(X, k)

    ### Applying the labels
    ena.metadata[!, :LABEL] .= "No Label"
    for (i, cluster) in enumerate(results.assignments)
        mylabel = "Auto Cluster #$(cluster)"
        ena.metadata[i, :LABEL] = mylabel
        if !haskey(colorMap, mylabel)
            if i <= length(EpistemicNetworkAnalysis.DEFAULT_EXTRA_COLORS)
                colorMap[mylabel] = EpistemicNetworkAnalysis.DEFAULT_EXTRA_COLORS[i]
            else
                colorMap[mylabel] = RGB(rand(0.4:0.01:0.9), rand(0.4:0.01:0.9), rand(0.4:0.01:0.9))
            end
        end
    end

    ### Saving labels back into original dataframe
    data[!, :LABEL] .= "No Label"
    for row in eachrow(ena.metadata)
        dataRows = [join(dataRow[ena.units], ".") == row[:ENA_UNIT] for dataRow in eachrow(data)]
        data[dataRows, :LABEL] = row[:LABEL]
    end

    ### Return score
    SSE = 0
    for (i, cluster) in enumerate(results.assignments)
        SSE += evaluate(Euclidean(), X[:, i], results.centers[:, cluster]) ^ 2
    end

    return SSE
end

K = 5
for k in 1:K
    autocluster!(ena, data, colorMap, k, ena.networkModel[!, :relationship])
    p = plot(ena, groupBy=:LABEL, showNetworks=true, weakLinks=false, leg=false)
    savefig(p, "images/SVD-$(k).png")
    p = plot(ena, groupBy=:LABEL, showNetworks=false, leg=false)
    savefig(p, "images/SVD-$(k)-NoNetwork.png")
end

# # DBSCAN approach
# epsval = 0.07
# min_cluster_size=10
# min_neighbors=2
# autocluster2!(ena, data, colorMap, epsval, min_cluster_size, min_neighbors, [:pos_x, :pos_y])
# p = plot(ena, groupBy=:LABEL, showNetworks=true, weakLinks=false, leg=false)
# savefig(p, "images/SVD-DBSCAN.png")
# p = plot(ena, groupBy=:LABEL, showNetworks=false, leg=false)
# savefig(p, "images/SVD-DBSCAN-NoNetwork.png")

# Quandrant approach
ena.metadata[!, :LABEL] = map(eachrow(ena.accumModel)) do row
    if row[:pos_y] > 0
        if row[:pos_x] > 0
            return "I"
        else
            return "II"
        end
    else
        if row[:pos_x] > 0
            return "IV"
        else
            return "III"
        end
    end
end

data[!, :LABEL] .= "No Label"
for row in eachrow(ena.metadata)
    dataRows = [join(dataRow[ena.units], ".") == row[:ENA_UNIT] for dataRow in eachrow(data)]
    data[dataRows, :LABEL] = row[:LABEL]
end

p = plot(ena, groupBy=:LABEL, showNetworks=true, weakLinks=false, leg=false)
savefig(p, "images/SVD-Quadrant.png")
p = plot(ena, groupBy=:LABEL, showNetworks=false, leg=false)
savefig(p, "images/SVD-Quadrant-NoNetwork.png")

# rotation = LDARotation(:LABEL)
# ena = ENAModel(data, codes, conversations, units, dropEmpty=dropEmpty, sphereNormalize=sphereNormalize, rotateBy=rotation)
# p = plot(ena, groupBy=:LABEL, showNetworks=true, weakLinks=false, leg=false, lims=0.25)
# savefig(p, "images/LDA-Quadrant.png")
# p = plot(ena, groupBy=:LABEL, showNetworks=false, leg=false, lims=0.25)
# savefig(p, "images/LDA-Quadrant-NoNetwork.png")

# # The Nonlinear Gamut
# seed = 4321
# knn = 35
# min_cluster_size=10
# min_neighbors=2
# colorMap = Dict("No Label" => colorant"black")
# function gamut(epsval, w)

#     ## ENA
#     ena = ENAModel(data, codes, conversations, units, dropEmpty=dropEmpty, sphereNormalize=sphereNormalize)

#     ## UMAP
#     model = embedUnits!(ena, :WordCount, knn, w, seed)
#     embedNetwork!(ena, model, seed)

#     ## DBSCAN
#     autocluster2!(ena, data, colorMap, epsval, min_cluster_size, min_neighbors)

#     ## Plotting
#     p = plot(ena, weakLinks=false, groupBy=:LABEL)
#     savefig(p, "images/SVD_UMAP.png")

#     p = plotUMAP(ena, colorMap, :WordCount, colormode=:label)
#     savefig(p, "images/LabelUMAP.png")

#     for (i, group) in enumerate(sort(unique(ena.metadata[!, :LABEL])))
#         if group != "No Label"
#             p = plotUMAP(ena, colorMap, :WordCount, colormode=:label, group=group)
#             savefig(p, "images/GroupedUMAP_$(i).png")
#         end
#     end
# end

# gamut(0.33, 1 / (nrow(ena.networkModel) + 1))

# cluster1rows = [row[:LABEL] == "Auto Cluster #1" for row in eachrow(ena.metadata)]
# cluster2rows = [row[:LABEL] == "Auto Cluster #2" for row in eachrow(ena.metadata)]
# cluster5rows = [row[:LABEL] == "Auto Cluster #5" for row in eachrow(ena.metadata)]
# cluster6rows = [row[:LABEL] == "Auto Cluster #6" for row in eachrow(ena.metadata)]
# result12 = MannWhitneyUTest(ena.accumModel[cluster1rows, :pos_y], ena.accumModel[cluster2rows, :pos_y])
# result15 = MannWhitneyUTest(ena.accumModel[cluster5rows, :pos_x], ena.accumModel[cluster1rows, :pos_x])
# result56 = MannWhitneyUTest(ena.accumModel[cluster5rows, :pos_x], ena.accumModel[cluster6rows, :pos_x])

println("Done!")

end # let


text:
using Zygote, Test
using Zygote: gradient, ZygoteRuleConfig
using CUDA
using CUDA: has_cuda

@testset "all" begin  # Overall testset ensures it keeps running after failure

  if has_cuda()
    @testset "CUDA tests" begin
      include("cuda.jl")
    end
    @info "CUDA tests have run"
  else
    @warn "CUDA not found - Skipping CUDA Tests"
  end

  @testset "Interface" begin
    include("interface.jl")
  end

  @testset "Tools" begin
    include("tools.jl")
  end

  @testset "Utils" begin
    include("utils.jl")
  end

  @testset "lib" begin
    include("lib/number.jl")
    include("lib/lib.jl")
    include("lib/array.jl")
  end

  @testset "Features" begin
    include("features.jl")
    @info "features.jl done"
  end

  @testset "Forward" begin
    include("forward/forward.jl")
  end

  @testset "Data Structures" begin
    include("structures.jl")
  end

  @testset "ChainRules" begin
    include("chainrules.jl")
    @info "chainrules.jl done"
  end

  @testset "Gradients" begin
    include("gradcheck.jl")
  end

  @testset "Complex" begin
    include("complex.jl")
  end

  @testset "Compiler" begin
    include("compiler.jl")
  end

end # @testset "all"


text:
# Julia wrapper for header: /usr/local/include/libavcodec/vorbis_parser.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_vorbis_parse_init,
    av_vorbis_parse_free,
    av_vorbis_parse_frame_flags,
    av_vorbis_parse_frame,
    av_vorbis_parse_reset

function av_vorbis_parse_init(extradata,extradata_size::Integer)
    ccall((:av_vorbis_parse_init,libavcodec),Ptr{AVVorbisParseContext},(Ptr{UInt8},Cint),extradata,extradata_size)
end

function av_vorbis_parse_free(s)
    ccall((:av_vorbis_parse_free,libavcodec),Void,(Ptr{Ptr{AVVorbisParseContext}},),s)
end

function av_vorbis_parse_frame_flags(s,buf,buf_size::Integer,flags)
    ccall((:av_vorbis_parse_frame_flags,libavcodec),Cint,(Ptr{AVVorbisParseContext},Ptr{UInt8},Cint,Ptr{Cint}),s,buf,buf_size,flags)
end

function av_vorbis_parse_frame(s,buf,buf_size::Integer)
    ccall((:av_vorbis_parse_frame,libavcodec),Cint,(Ptr{AVVorbisParseContext},Ptr{UInt8},Cint),s,buf,buf_size)
end

function av_vorbis_parse_reset(s)
    ccall((:av_vorbis_parse_reset,libavcodec),Void,(Ptr{AVVorbisParseContext},),s)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
"""
Package: Forecast

    summarize(x; varnames = nothing)

Return statistical summary for x

The values returned are dividen in three sections, the first one shows Minimum, 1st Quantile, Median, Mean, 3rd Quantile, Maxixum and the p-value for the Jarque-Bera Normality Test. The second one show the first four moment; Mean, Variance, Skewness and Kurtosis, an finally a summary with the different types contained in the Array.

# Arguments
- `x`: Array or DataFrame of data.
- `varnames`: Names for the columns to be summarized, it defaults to automatic naming
            or the existing names in when a DataFrame.

# Returns
A SUMMARIZE struct

# Examples
```julia-repl
julia> summarize(rand(100,3); varnames = ["a","b","c"])
┌──────────┬────────────┬──────────┬──────────┬──────────┬──────────┬──────────┬───────────┐
│ Variable │        Min │       1Q │   Median │     Mean │       3Q │      Max │ H0 Normal │
├──────────┼────────────┼──────────┼──────────┼──────────┼──────────┼──────────┼───────────┤
│        a │ 0.00520465 │ 0.205712 │ 0.462199 │ 0.465784 │   0.6913 │  0.97946 │ 0.0593599 │
│        b │ 0.00218787 │ 0.247344 │ 0.485465 │ 0.498587 │ 0.723371 │ 0.985226 │ 0.0562301 │
│        c │  0.0244256 │ 0.247598 │ 0.530821 │ 0.498689 │ 0.722731 │ 0.967952 │ 0.0356495 │
└──────────┴────────────┴──────────┴──────────┴──────────┴──────────┴──────────┴───────────┘
┌──────────┬──────────┬───────────┬───────────┬───────────┐
│ Variable │     Mean │  Variance │  Skewness │  Kurtosis │
├──────────┼──────────┼───────────┼───────────┼───────────┤
│        a │ 0.465784 │ 0.0823949 │ 0.0823949 │ 0.0823949 │
│        b │ 0.498587 │ 0.0854883 │ 0.0854883 │ 0.0854883 │
│        c │ 0.498689 │ 0.0790597 │ 0.0790597 │ 0.0790597 │
└──────────┴──────────┴───────────┴───────────┴───────────┘
┌──────────┬─────────┐
│ Variable │ Float64 │
├──────────┼─────────┤
│        a │     100 │
│        b │     100 │
│        c │     100 │
└──────────┴─────────┘
```
"""
function summarize(x::DataFrame)
    summarize(Array(tots(x)[:,2:end]), varnames = names(x)[2:end])
end

function summarize(x::AbstractVector{<:Real}; varnames = nothing)
    summarize(reshape(x,length(x),1), varnames = varnames)
end

function summarize(x::Real; varnames = nothing)
    summarize([x,x], varnames = varnames)
end

function summarize(x::AbstractArray{T}; varnames = nothing) where T<:Real

    x = compact(x)
    @assert (nd = ndims(x)) <= 2 "Data should have two dimensions at most."
    n,m = nd == 1 ? (length(x),1) : size(x)

    varnames = isnothing(varnames) ? ["x"*string(i) for i in 1:m] : varnames

    # Quantiles and Mean
    sk = ["Variable", "Min","1Q","Median","Mean","3Q","Max","H0 Normal"]
    sv = zeros(m,7)
    s = Array{Any,2}(undef,m,7)

    minimum_Number(x) = minimum(filter(x -> x isa Number, x))
    maximum_Number(x) = maximum(filter(x -> x isa Number, x))
    mean_Number(x) = mean(filter(x -> x isa Number, x))
    quantile_Number(x,q) = quantile(filter(x -> x isa Number, x),q)
    normality(x) = pvalue(JarqueBeraTest(x))

    sv[nd <= 1 ? 1 : 1:end,7] = nd <= 1 ? normality(x) : mapslices(normality,x,dims=[1])
    sv[nd <= 1 ? 1 : 1:end,1] = nd <= 1 ? minimum(x) : mapslices(minimum_Number,x,dims=[1])
    sv[nd <= 1 ? 1 : 1:end,6] = nd <= 1 ? maximum(x) : mapslices(maximum_Number,x,dims=[1])
    sv[nd <= 1 ? 1 : 1:end,4] = nd <= 1 ? mean(x) : mapslices(mean_Number,x,dims=[1])
    sv[nd <= 1 ? 1 : 1:end,[2,3,5]] = nd <= 1 ? quantile_Number(x,[.25,.5,.75]) :
                             mapslices(x -> quantile_Number(x,[.25,.5,.75]),x,dims=[1])'
    qua = DataFrame(hcat(varnames,sv), :auto)
    DataFrames.rename!(qua, sk)

    # Moments
    sk = ["Variable", "Mean", "Variance", "Skewness", "Kurtosis"]
    sv = zeros(m,4)
    s = Array{Any,2}(undef,m,4)

    var_Number(x) = var(filter(x -> x isa Number, x))
    skewness_Number(x) = var(filter(x -> x isa Number, x))
    kurtosis_Number(x) = var(filter(x -> x isa Number, x))

    sv[nd <= 1 ? 1 : 1:end,1] = nd <= 1 ? mean(x) : mapslices(mean_Number,x,dims=[1])
    sv[nd <= 1 ? 1 : 1:end,2] = nd <= 1 ? var(x) : mapslices(var_Number,x,dims=[1])
    sv[nd <= 1 ? 1 : 1:end,3] = nd <= 1 ? skewness(x) : mapslices(skewness_Number,x,dims=[1])
    sv[nd <= 1 ? 1 : 1:end,4] = nd <= 1 ? kurtosis(x) : mapslices(kurtosis_Number,x,dims=[1])

    mo = DataFrame(hcat(varnames,sv), :auto)
    DataFrames.rename!(mo, sk)

    # Format Check
    tox = string.(typeof.(x))
    tx = sort(unique(tox))
    nt = length(tx)
    mv = length(varnames)
    vt = Array{Any}(undef,mv,nt)
    vt[:] .= 0

    for (i,ty) in enumerate(tx)
        vt[:,i] = count(==(ty), tox, dims=1)
    end

    tdf = DataFrame(hcat(varnames,vt),:auto)
    fc = DataFrames.rename!(tdf,vcat("Variable",  tx))

    SUMMARIZE(qua,mo,fc)

end


text:
#!/usr/bin/env julia

import Pkg3: TOML

prefix = joinpath(homedir(), ".julia", "registries", "Stdlib")

# TODO: use Sys.STDLIBDIR instead once implemented
let vers = "v$(VERSION.major).$(VERSION.minor)"
global stdlibdir = realpath(abspath(Sys.BINDIR, "..", "share", "julia", "stdlib", vers))
isdir(stdlibdir) || error("stdlib directory does not exist: $stdlibdir")
end
juliadir = dirname(stdlibdir)

stdlib_uuids = Dict{String,String}()
stdlib_trees = Dict{String,String}()
stdlib_deps = Dict{String,Vector{String}}()

for pkg in readdir(stdlibdir)
    project_file = joinpath(stdlibdir, pkg, "Project.toml")
    isfile(project_file) || continue
    project = TOML.parsefile(project_file)
    stdlib_uuids[pkg] = project["uuid"]
    stdlib_trees[pkg] = split(readchomp(`git -C $juliadir ls-tree HEAD -- stdlib/$pkg`))[3]
    stdlib_deps[pkg] = String[]
    haskey(project, "deps") || continue
    append!(stdlib_deps[pkg], sort!(collect(keys(project["deps"]))))
end

#=
write_toml(prefix, "Registry") do io
    repo = "https://github.com/JuliaRegistries/Stdlib.git"
    uuid = string(uuid5(uuid_registry, repo))
    println(io, "name = ", repr("Stdlib"))
    println(io, "uuid = ", repr(uuid))
    println(io, "repo = ", repr(repo))
    println(io, "\ndescription = \"\"\"")
    print(io, """
        Official Julia Standard Library registry where versions
        of packages that ship with Julia itself are registered.
        """)
    println(io, "\"\"\"")
    println(io, "\n[packages]")
    for pkg in sort!(collect(keys(stdlib_deps)), by=pkg->stdlib_uuids[pkg])
        println(io, stdlib_uuids[pkg], " = { name = ", repr(pkg), ", path = ", repr(pkg), " }")
    end
end

for (pkg, uuid) in stdlib_uuids
    url = "https://github.com/JuliaLang/julia.git"
    tree = stdlib_trees[pkg]
    deps = stdlib_deps[pkg]

    # Package.toml
    write_toml(prefix, pkg, "Package") do io
        println(io, "name = ", repr(pkg))
        println(io, "uuid = ", repr(uuid))
        println(io, "repo = ", repr(url))
    end

    # Versions.toml
    write_toml(prefix, pkg, "Versions") do io
        println(io, "[", toml_key("0.7.0-DEV+r$(tree[1:8])"), "]")
        println(io, "tree.git.sha1 = ", repr(tree))
    end

    # Deps.toml
    if isempty(deps)
        rm(joinpath(prefix, pkg, "Deps.toml"), force=true)
    else
        write_toml(prefix, pkg, "Deps") do io
            println(io, "[", toml_key("0.7"), "]")
            for dep in sort!(deps, by=lowercase)
                println(io, "$dep = ", repr(stdlib_uuids[dep]))
            end
        end
    end

    # Compat.toml
    write_toml(prefix, pkg, "Compat") do io
        println(io, "[", toml_key("0.7"), "]")
        println(io, "julia = \"0.7\"")
    end
end
=#

max_stars_repo_name: ederag/julia

max_issues_repo_name: ederag/julia

max_forks_repo_name: ederag/julia


text:
struct WendlandC8{T} <: SPHKernel
    n_neighbours::Int64
    norm_2D::T
    norm_3D::T
end

"""
    WendlandC8(T::DataType=Float64, n_neighbours::Integer=395)

Set up a `WendlandC8` kernel for a given DataType `T`.
"""
WendlandC8(T::DataType=Float64, n_neighbours::Integer=395) = WendlandC8{T}(n_neighbours, 8.0/(3π), 357.0/(64π))

"""
    kernel_value_2D(kernel::WendlandC8{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_2D(kernel::WendlandC8{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^2
        t1 = 1 - u
        t9 = t1*t1*t1*t1*t1*t1*t1*t1*t1  # (1.0 - u)^9
        u2 = u*u
        return ( t1 * t9 * (5 + 50u + 210u2 + 450u2 * u + 429u2 * u2)) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_2D(kernel::WendlandC8{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_2D(kernel::WendlandC8{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^3
        t1 = 1.0 - u
        t9 = t1*t1*t1*t1*t1*t1*t1*t1*t1  # (1.0 - u)^9
        u2 = u*u
        return -26t9 * u * (231u2 * u + 159u2 + 45u + 5 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_2D(kernel::WendlandC8, density::Real, m::Real, h_inv::Real)

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
function bias_correction_2D(kernel::WendlandC8{T}, density::Real, m::Real, h_inv::Real) where T
    return density |> T
end

"""
    kernel_value_3D(kernel::WendlandC8{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC8 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_3D(kernel::WendlandC8{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^3
        t1 = (1 - u)
        t9 = t1*t1*t1*t1*t1*t1*t1*t1*t1  # (1.0 - u)^9
        u2 = u*u
        return ( t1 * t9 * (5 + 50u + 210u2 + 450u2 * u + 429u2 * u2)) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_3D(kernel::WendlandC8, u::Real, h_inv::Real)

Evaluate the derivative of the WendlandC8 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_3D(kernel::WendlandC8{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^2 * h_inv^2
        t1 = 1 - u
        t9 = t1*t1*t1*t1*t1*t1*t1*t1*t1  # (1.0 - u)^9
        u2 = u*u
        return -26t9 * u * (231u2 * u + 159u2 + 45u + 5 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_3D(kernel::WendlandC8{T}, density::Real, m::Real, h_inv::Real) where T

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
function bias_correction_3D(kernel::WendlandC8{T}, density::Real, m::Real, h_inv::Real) where T
    return density |> T
end


text:
import AQFED.Math: norminv
import Random: rand!, randn!, rand, randn
using Statistics
import AQFED.TermStructure:
    LocalVolatilityModel, VarianceSurfaceBySection, localVarianceByLogmoneyness, discountFactor, logForward

    function ndims(model::LocalVolatilityModel, specificTimes::Vector{Float64}, timestepSize::Float64)
        genTimes = pathgenTimes(model, specificTimes, timestepSize)
        return (length(genTimes) - 1) #0.0 does not count
    end

function simulate(
    rng,
    model::LocalVolatilityModel{S},
    spot::Float64,
    payoff::VanillaOption,
    start::Int,
    nSim::Int,
    timestepSize::Float64;
    withBB = false
) where {S}
    specTimes = specificTimes(payoff)
    tte = specTimes[end]
    df = discountFactor(model, tte)
    genTimes = pathgenTimes(model, specTimes, timestepSize) #LinRange(0.0, tte, ceil(Int, nSteps * tte) + 1)
    local bb, cache
    if withBB
        cacheSize = ceil(Int, log2(length(genTimes))) * 4
        bb = BrownianBridgeConstruction(genTimes[2:end])
        cache = BBCache{Int,Vector{Float64}}(cacheSize)
    end
    logpathValues = Vector{Float64}(undef, nSim)
    z = Vector{Float64}(undef, nSim)
    t0 = genTimes[1]
    lnspot =  log(spot)
    lnf0 = logForward(model, lnspot, t0)
    lnforward = lnf0
    logpathValues .= lnf0
    local payoffValues
    #println("genTimes ",genTimes)
    for (dim,t1) in enumerate(genTimes[2:end])
        h = t1 - t0
        sqrth = sqrt(h)
        if withBB
            transformByDim!(bb, rng, start, dim, z, cache)
        else
            skipTo(rng, dim, start)
            nextn!(rng, dim, z)
            @. z *= sqrth
        end
        #z = randn(rng, Float64, nSim)
        lnf1 = logForward(model, lnspot, t1)
        lnforward += lnf1 - lnf0
        # indexT0 = searchsortedlast(model.surface.expiries, t0)
        # indexT1 = searchsortedlast(model.surface.expiries, t1)

        # lvh = map(
        #     x -> localVarianceByLogmoneyness(model.surface, x - lnforward, t0, t1, indexT0, indexT1) * h,
        #     logpathValues,
        # )
        logmoneyness = @. logpathValues - lnforward
        lv = localVarianceByLogmoneyness(model.surface, logmoneyness, t0, t1)
        @. logpathValues += lnf1 - lnf0 - lv / 2 * h  + sqrt(lv) * z
        if t1 == tte
            pathValues = z #reuse Var
            @. pathValues = exp(logpathValues)
            payoffValues = map(x -> evaluatePayoff(payoff, x, df), pathValues)
        end
        t0 = t1
        lnf0 = lnf1
    end
    payoffMean = mean(payoffValues)
    return payoffMean, stdm(payoffValues, payoffMean) / sqrt(length(payoffValues))
end


text:
#== import modules here ==#
using PyCall
psi4 = pyimport("psi4")

using HDF5
#using MPI

#== psi4 setup ==#
psi4.set_memory("500 MB")
psi4.core.set_output_file("atom_energy.dat")

#===============#
#== MPI setup ==#
#===============#
#MPI.Init()
#comm = MPI.COMM_WORLD

#mpi_rank = MPI.Comm_rank(comm)
#mpi_size = MPI.Comm_size(comm)

#==============================#
#== get atom/basis set pairs ==#
#==============================#

pairs = Vector{String}([])
status_bsed = h5open("../records/bsed.h5","r") do bsed
  #== loop over atoms ==#
  for symbol in keys(bsed)
    #== loop over basis sets in atom ==#
    for basis in keys(bsed["$symbol"])
      #== add atom/basis-set pair to pair list ==#
      pair = "$symbol/$basis"
      push!(pairs, pair)
    end
  end
end

status_eatom = h5open("eatom.h5","w") do eatom
  #if mpi_rank == 0
    flush(eatom)
  #end

  #symbol = "O"
  #basis = "6-31G"
  #pair = "$symbol/$basis"

  #MPI.Barrier(comm)
  #== mp2 energies ==#
  for (pair_idx, pair) in enumerate(pairs)
    #if mpi_rank != (pair_idx-1)%mpi_size continue end

    pair_regex = match(r"(.*)/(.*)", pair)
    symbol = pair_regex.captures[1]
    basis = pair_regex.captures[2]

    if symbol == "H"
    println("$symbol/$basis")

    try
      geom = psi4.geometry("""
        $symbol  0.0 0.0 0.0
        symmetry=c1
      """)

      options = Dict(
        "REFERENCE" => "ROHF",
        "BASIS" => basis,
        "SCF_TYPE"  => "PK",
        "GUESS" => "SAD",
	      "PUREAM" => false,
        #"SAD_SCF_TYPE" => "DIRECT",
        "E_CONVERGENCE" => 1e-10,
        "D_CONVERGENCE" => 1e-10,
        #"MAXITER" => 1,
        "FAIL_ON_MAXITER" => false,
        "MP2_TYPE" => "CONV",
        "QC_MODULE" => "OCC"
      )
      psi4.set_options(options)

      mp2_name = "mp2/"*options["BASIS"]
      mp2_e, mp2_wfn = psi4.energy(mp2_name, molecule=geom, return_wfn=true)
      write(eatom, "RIMP2/"*pair, mp2_e)
    catch e
      bt = catch_backtrace()
      msg = sprint(showerror, e, bt)
      println(msg)
      continue
    end
    end
  end

  #== rhf energies ==#
  for (pair_idx, pair) in enumerate(pairs)
    #if mpi_rank != (pair_idx-1)%mpi_size continue end

    pair_regex = match(r"(.*)/(.*)", pair)
    symbol = pair_regex.captures[1]
    basis = pair_regex.captures[2]

    println("$symbol/$basis")

    if symbol == "H"
    try
      geom = psi4.geometry("""
        $symbol  0.0 0.0 0.0
        symmetry=c1
      """)

      options = Dict(
        "REFERENCE" => "ROHF",
        "BASIS" => basis,
        "SCF_TYPE"  => "PK",
        "GUESS" => "SAD",
	      "PUREAM" => false,
        #"SAD_SCF_TYPE" => "DIRECT",
        "E_CONVERGENCE" => 1e-10,
        "D_CONVERGENCE" => 1e-10,
        #"MAXITER" => 1,
        "FAIL_ON_MAXITER" => false
      )
      psi4.set_options(options)

      scf_name = "scf/"*options["BASIS"]
      scf_e, scf_wfn = psi4.energy(scf_name, molecule=geom, return_wfn=true)
      write(eatom, "RHF/"*pair, scf_e)
    catch e
      bt = catch_backtrace()
      msg = sprint(showerror, e, bt)
      println(msg)
      continue
    end
    end
  end
end

#MPI.Finalize()


text:
using Bridge
using Test, Statistics, LinearAlgebra, Random

# tests with alpha 0.01
# test fail 1% of the time
r = 2.576
n = 1000
chiupper = 1118.95 #upper 0.005 percentile for n = 1000
chilower = 888.56 #lower 0.005 percentile for n = 1000

# varmu(x, mu) = dot(x-mu, x-mu)/length(x)
var0(x) = dot(x, x)/length(x)

brown1(s, t, n) = sample(range(s, stop=t, length=n), Wiener())
bb(u,v,t,n) = sample(range(0, stop=t, length=n), WienerBridge{Float64}(t,v), u)
X = SamplePath([0.0, 1.0, 3.0], zeros(3))
@test abs(mean(Bridge.sample!(Bridge.TransitionProb(), X,  Wiener(), 0.0).yy[end] for i in 1:n) -
    mean(Bridge.transitionprob(0.0, 0.0, 3.0, Wiener()))) <
    3*std(Bridge.transitionprob(0.0, 0.0, 3.0, Wiener()))/sqrt(n)

X = sample(range(0.0, stop=2.0, length=1000), Wiener())
@test bracket(X)[end] == bracket(X,X)[end] == quvar(X)

#quadratic variation of Brownian motion is proportional to time plus sampling bias
quv = [quvar(sample(range(0.0, stop=2.0, length=1000), Wiener())) for j in 1:1000]
s2 = var(quv)
@test abs(mean(quv)-2)/sqrt(s2)*sqrt(1000) < 3

# int0^T w_t dw_t = w_T^2/2 - T/2
@test abs((b -> (ito(b, b).yy[end] - (0.5b.yy[end]^2 - 1)))(brown1(0, 2, 10000))) < 0.1

struct Diff <: Bridge.ContinuousTimeProcess{Float64}
end
import Bridge: b, σ
Bridge.b(t,x, ::Diff) = -5x
Bridge.σ(t, x, ::Diff) = 1.

X = [solve(EulerMaruyama(), 0.0, brown1(0.0,1.0,1000), WienerBridge(2.0, 1.0)).yy[end] for i in 1:1000]
@test abs(mean(X)-0.5) < r*sqrt(2/n)

#test (roughly) the quadratic variance of Euler approximation
quv = [quvar(solve(EulerMaruyama(), 0.0, brown1(0.0,1.0,1000), Diff())) for j in 1: 1000]
s2 = var(quv)
@test abs(mean(quv)-1)/sqrt(s2)*sqrt(1000) < 5.

#  B(t) = (1-t) W(t/(1-t)).

@test abs(mean([bb(0,1,1,3).yy[2] for i in 1:1000])-0.5)/sqrt(0.25)*sqrt(1000) < r
@test abs(mean([bb(0,1,1,5).yy[3] for i in 1:1000])-0.5)/sqrt(0.25)*sqrt(1000) < r
@test chilower < 1000*abs(var0([bb(0,1,1,3).yy[2] for i in 1:1000] .- 0.5))/.25 < chiupper
@test chilower < 1000*abs(var0([bb(0,1,1,5).yy[3] for i in 1:1000] .- 0.5))/.25 < chiupper

# Covariance of a Brownian bridge from t_1 to t_2 (t_2-t)(s-t_1)/(t_2-t_1)
# here (2-1)(1)/2 = 1/2

@test chilower < 1000*abs(var0([bb(0,1,2,5).yy[3] for i in 1:1000] .- 0.5))*2 < chiupper


text:
global bulk_fields = ["IEST_BRAND_PRODUCT_LIST","BLOOMBERG_PEERS",
                      "HIST_TRR_MONTHLY","BEST_ANALYST_RECS_BULK",
                      "ERN_ANN_DT_AND_PER","EARN_ANN_DT_TIME_HIST_WITH_EPS",
                      "INDX_MEMBERS","INDX_MWEIGHT","INDX_MWEIGHT_HIST",
                      "PORTFOLIO_MEMBERS","PORTFOLIO_MPOSITION",
                      "PORTFOLIO_MWEIGHT","PORTFOLIO_DATA"]

function flatten_all(a)
    while any(x-> !(typeof(x)==String) && (typeof(x)<:Array || typeof(x)<:Tuple), a)
        a = vcat([typeof(x)==String ? [x] : collect(Base.Iterators.flatten([x])) for x in a]...)
    end
    a
end

function date(val)
    # BLPAPI_DATETIME_YEAR_PART | ..._MONTH_PART | ..._DAY_PART
    @assert val.parts == 0x01 | 0x02 | 0x04
    Dates.DateTime(val.year,val.month,val.day)
end

function datetime(val)
    # BLPAPI_DATETIME_YEAR_PART | ..._MONTH_PART | ..._DAY_PART | ..._HOURS_PART | ..._MINUTES_PART | ..._SECONDS_PART | ..._MILLISECONDS_PART
    @assert val.parts == 0x01 | 0x02 | 0x04 | 0x10 | 0x20 | 0x40 | 0x80
    Dates.DateTime(val.year,val.month,val.day,val.hours,val.minutes,val.seconds,val.milliSeconds)
end

function decode_value(element,datatype,i)
    @Match.match datatype begin
        1 => begin val=Ref{Int32}(0); blpapi_Element_getValueAsBool(element,val,i); Bool(val[]) end
        2 => error("BLPAPI_DATATYPE_CHAR")
        3 => error("BLPAPI_DATATYPE_BYTE")
        4 => begin val=Ref{Int32}(0); blpapi_Element_getValueAsInt32(element,val,i); val[] end
        5 => begin val=Ref{Int64}(0); blpapi_Element_getValueAsInt64(element,val,i); val[] end
        6 => begin val=Ref{Float32}(0); blpapi_Element_getValueAsFloat32(element,val,i); val[] end
        7 => begin val=Ref{Float64}(0); blpapi_Element_getValueAsFloat64(element,val,i); val[] end
        8 => begin val=Ref{Ptr{UInt8}}(0); blpapi_Element_getValueAsString(element,val,i); unsafe_string(val[]) end
        9 => error("BLPAPI_DATATYPE_BYTEARRAY")
        10 => begin val=blpapi_Datetime(0,0,0,0,0,0,0,0,0); blpapi_Element_getValueAsDatetime(element,val,i); date(val); end
        11 => error("BLPAPI_DATATYPE_TIME")
        12 => error("BLPAPI_DATATYPE_DECIMAL")
        13 => begin val=blpapi_Datetime(0,0,0,0,0,0,0,0,0); blpapi_Element_getValueAsDatetime(element,val,i); datetime(val); end
    end
end

function decode_element(element)
    name=blpapi_Element_nameString(element)
    datatype = blpapi_Element_datatype(element)
    if datatype <= 13
        (name,[decode_value(element,datatype,i) for i in 0:(blpapi_Element_numValues(element)-1)])
    elseif datatype == 14
        error("BLPAPI_DATATYPE_ENUMERATION")
    elseif datatype == 15
        if Bool(blpapi_Element_isArray(element))
            (name,[decode_element(@with_pointer val blpapi_Element_getValueAsElement(element,val,i)) for i in 0:(blpapi_Element_numValues(element)-1)])
        else
            (name,[decode_element(@with_pointer val blpapi_Element_getElementAt(element,val,i)) for i in 0:(blpapi_Element_numElements(element)-1)])
        end
    elseif datatype == 16
        decode_element(@with_pointer choice blpapi_Element_getChoice(element,choice))
    elseif datatype == 17
        error("BLPAPI_CORRELATION_ID")
    end
end

function process_message(msg)
    messageType = blpapi_Message_typeString(msg)
    contents = blpapi_Message_elements(msg)
    #@show AbstractTrees.Tree(decode_element(contents))
    if blpapi_Element_hasElement(contents,"responseError",C_NULL)
        response_error = process_tree(decode_element(@with_pointer err blpapi_Element_getElement(contents,err,"responseError",C_NULL)))
        error(response_error["message"])
    end
    return decode_element(contents)
end

function process_event(event)
    iterator = blpapi_MessageIterator_create(event)
    done = false
    out = []
    try
        while !done
            msg = @with_pointer msg blpapi_MessageIterator_next(iterator,msg)
            if msg == C_NULL
                done = true
            else
                try
                    out = vcat(out,process_message(msg))
                finally
                    blpapi_Message_release(msg)
                end
            end
        end
    finally
        blpapi_MessageIterator_destroy(iterator)
    end
    out
end

function retrieve_response(timeout=10)
    done = false
    out = []
    while !done
        event = @with_pointer event blpapi_Session_nextEvent(session,event,1000*timeout)
        try
            event_type = @Match.match blpapi_Event_eventType(event) begin
                1 => :ADMIN
                2 => :SESSION_STATUS
                3 => :SUBSCRIPTION_STATUS
                4 => :REQUEST_STATUS
                5 => :RESPONSE
                6 => :PARTIAL_RESPONSE
                8 => :SUBSCRIPTION_DATA
                9 => :SERVICE_STATUS
                10 => :TIMEOUT
                11 => :AUTHORIZATION_STATUS
                12 => :RESOLUTION_STATUS
                13 => :TOPIC_STATUS
                14 => :TOKEN_STATUS
                15 => :REQUEST
            end
            if event_type == :TIMEOUT || event_type == :RESPONSE
                done = true
            end
            if event_type == :RESPONSE || event_type == :PARTIAL_RESPONSE
                out = vcat(out,process_event(event))
            elseif event_type == :TIMEOUT
                @warn (event_type,)
            else
                @warn (event_type,flatten_all(process_event(event)))
            end
        finally
            blpapi_Event_release(event)
        end
    end
    out
end

function security_data(x)
    global bulk_fields
    #@show AbstractTrees.Tree(x)
    if length(x["fieldData"]) == 0
        (x["security"][1],nothing)
    elseif x["fieldData"][1][1] in bulk_fields
        if length(x["fieldData"]) > 1
            error("Bulk data fields have to be retrieved alone")
        else
            (x["security"][1],
             DataFrames.DataFrame(map(a->Dict([Symbol(k)=>process_tree(v) for (k,v) in a[2]]), x["fieldData"][1][2])))
        end
    else
        merge(Dict(:security=>x["security"][1]),
              Dict([Symbol(k)=> process_tree(v) for (k,v) in x["fieldData"]]))
    end
end

function security_data_hist(x)
    df = DataFrames.DataFrame(map(row->Dict([Symbol(k)=> process_tree(v) for (k,v) in row[2]]),x["fieldData"]))
    if size(df,1) > 0
        (x["security"][1],TimeSeries.TimeArray(df,timestamp=:date))
    else
        (x["security"][1],nothing)
    end
end

function bar_data(x)
    rows = map(x->x[2],x["barTickData"])
    DataFrames.DataFrame(map(row->Dict([Symbol(k)=>process_tree(v) for (k,v) in row]),rows))
end

function process_tree(x::Any,context="")
    if typeof(x) == Int32 && x == -2147483648
        NaN
    elseif typeof(x) == Float64 && x == -2.4245362661989844e-14
        NaN
    else
        x
    end
end

function process_tree(x::Tuple,context="")
    @Match.match x[1] begin
        "securityData" => context=="H" ? security_data_hist(Dict(x[2])) : map(x->security_data(Dict(x[2])),x[2])
        "barData" => bar_data(Dict(x[2]))
        "responseError" => Dict([k=>process_tree(v) for (k,v) in x[2]])
        "data" => map(a->Dict([Symbol(k)=>process_tree(v) for (k,v) in a[2][3][2]]),x[2][2][2])
        _ => @warn ["cannot handle",x]
    end
end

function process_tree(x::Array,context="")
    if length(x) == 1
        process_tree(x[1],context)
    else
        map(x->process_tree(x,context),x)
    end
end


text:
#cauchymp.jl
#Algorithm 4.1 of Random Eigenvalues by Alan Edelman

#Experiment:    Calculate Cauchy transform as trace of resolvent of Wishart matrix
#Theory:        Marcenko-Pastur law

m = 2000            # larger dimension of matrix
n = 1000            # smaller dimension of matrix
r = n/m             # aspect ratio
a = (1-sqrt(r))^2
b = (1+sqrt(r))^2
z = 3               # should be outside [a,b]
X = randn(m,n)
W = X'*X/m          # Generate Wishart matrix
println(( trace(inv(z*eye(n)-W))/n,
          (z-1+r-sqrt((z-a)*(z-b)))/(2*r*z) ))


text:
"""
    Labor demand by A firms.
"""
function labor_demand_A(k_I, par)
    @unpack_Pars_Simple par
    ((γ * α_hat * v * C^α_C * k_I^δ_hat) / w)^(1 / (1 - α_hat))
end

"""
    Labor demand by B firms.
"""
function labor_demand_B(k_I, par)
    @unpack_Pars_Simple par
    (((1 - γ) * α_hat * v * C^α_C * k_I^δ_hat) / w)^(1 / (1 - α_hat))
end

"""
    How much capital A firms have in place given their production and selling decision.
"""
total_capital_A(k_I_b, k_I_P, par) = (1 - par.ϕ) * k_I_b + k_I_P

"""
    How much capital B firms have in place given their production and selling decision.
"""
total_capital_B(k_I_b, k_I_P) = k_I_b + k_I_P

"""
    Market clearing condition for traded capital. For any unit shipped, only 1-τ units arrive. Needs to be scaled by sector size (γ for A and 1-γ for B).
"""
market_clearing(k_I_b_B, k_I_b_A, par) = (1 - par.γ) * k_I_b_B + (1 - par.τ) * par.γ * k_I_b_A

"""
    Marginal cost of capital production by A firms.
"""
mc_A(k_I, par) = par.ca_A * par.cb_A * k_I^(par.cb_A - 1)

"""
    Marginal cost of capital production by B firms.
"""
mc_B(k_I, par) = par.ca_B * par.cb_B * k_I^(par.cb_B - 1)

"""
    Cost of capital production by A firms.
"""
c_A(k_I, par) = par.ca_A * k_I^par.cb_A
"""
    Cost of capital production by B firms.
"""
c_B(k_I, par) = par.ca_B * k_I^par.cb_B

"""
    How much capital is produced by A firms if they are indifferent between
    producing one more unit or selling one unit less.

        AND

    The firm does not sell all of its capital! Otherwise producing one more unit has a different payoff (sell and keep the rest instead of keeping it fully).
"""
k_I_prod_A(par) = ((par.p_I / (1 - par.ϕ)) / (par.ca_A * par.cb_A))^(1 / (par.cb_A - 1))

"""
    How much capital is produced by B firms if they are indifferent between
    producing one more unit or selling one unit less.

        AND

    The firm does not sell all of its capital! Otherwise producing one more unit has a different payoff (sell and keep the rest instead of keeping it fully).
"""
k_I_prod_B(par) = ((par.p_I / (1 - par.τ)) / (par.ca_B * par.cb_B))^(1 / (par.cb_B - 1))

"""
    The price for the B-good. C_A is the quantity of A-goods.
"""
function price_A(C_A, par::Pars_Simple)
    @unpack_Pars_Simple par
    γ * v * C^α_C * C_A^(-1 / σ)
end

"""
    The price for the B-good. C_B is the quantity of B-goods.
"""
function price_B(C_B, par::Pars_Simple)
    @unpack_Pars_Simple par
    (1 - γ) * v * C^α_C * C_B^(-1 / σ)
end

"""
    The wage rate that leads to market clearing (l_A + l_B = 1)
"""
function wage_rate(k_I_A, k_I_B; par)
    @unpack_Pars_Simple par
    c = (α_hat * v * C^α_C)
    c * ((γ * k_I_A^δ_hat)^(1 / (1 - α_hat)) + ((1 - γ) * k_I_B^δ_hat)^(1 / (1 - α_hat)))^(1 - α_hat)
end

"""
    Profit for both type of firms.
"""
function profit(par)
    @unpack_Pars_Simple par
    k_I_b_A, k_I_b_B = solve_firm_problems(par)

    if par.ϕ > 0 && loss_A(k_I_b_A, par) > 1e-8
        println("")
        k_I_A, l_A = capital_and_labor_con_A(k_I_b_A, par)
    else
        k_I_A, l_A = capital_and_labor_unc_A(k_I_b_A, par)
    end
    k_I_B, l_B = capital_and_labor_B(k_I_b_B, par)

    k_I_P_A = k_I_A - (1 - ϕ) * k_I_b_A # Sold capital is only partially lost.
    k_I_P_B = k_I_B - k_I_b_B

    profit_A = γ * v * C^α_C * l_A^α_hat * k_I_A^δ_hat - w * l_A - p_I * k_I_b_A - c_A(k_I_P_A, par)
    profit_B = (1 - γ) * v * C^α_C * l_B^α_hat * k_I_B^δ_hat - w * l_B - p_I / (1 - τ) * k_I_b_B - c_B(k_I_P_B, par)

    return profit_A, profit_B
end

# function labor_demand(k_I_j, par::Pars_Simple)
#     @unpack_Pars_Simple par
#     (α_hat * γ * C^α_C * (k_I_j)^δ_hat / w)^(1 / (1 - α_hat))
# end


text:
# Stacks

mutable struct Stack{T}
  idx::Int
  data::Vector{T}
end

Stack(data::Vector{T}) where T =
  Stack{T}(length(data), data)

function Base.pop!(stk::Stack)
  i = stk.idx
  stk.idx = i == 1 ? length(stk.data) : i-1
  @inbounds return stk.data[i]
end

function _push!(a::Vector{T}, x::T) where T
  Base._growend!(a, 1)
  @inbounds a[end] = x
  return
end

# Emit

xstack(T) = stmt(Expr(:call, Vector{T}), type = Vector{T})

function alphauses(b)
  us = Set{Alpha}()
  postwalk(x -> x isa Alpha && push!(us, x), b)
  return us
end

xtuple(xs...) = xcall(:tuple, xs...)

concrete(T::DataType) = T
concrete(::Type{Type{T}}) where T = typeof(T)
concrete(T) = Any

runonce(b) = b.id in (1, length(b.ir.blocks))

function forward_stacks!(adj, F)
  stks, recs = [], []
  pr = adj.primal
  for b in blocks(pr), α in alphauses(block(adj.adjoint, b.id))
    if runonce(b)
      push!(recs, Variable(α))
    else
      T = exprtype(pr, Variable(α))
      stk = pushfirst!(pr, xstack(T))
      push!(recs, stk)
      push!(b, xcall(PartialP, :_push!, stk, Variable(α)))
    end
    push!(stks, (b.id, alpha(α)))
  end
  args = arguments(pr)[3:end]
  T = Tuple{concrete.(exprtype.((pr,), recs))...}
  isconcretetype(T) || (T = Any)
  rec = push!(pr, xtuple(recs...))
  if usetyped && length(pr.blocks) > 1
    rec = push!(pr, Expr(:call, Pullback{F,T}, rec))
  else
    P = length(pr.blocks) == 1 ? Pullback{F} : Pullback{F,Any}
    rec = push!(pr, Expr(:call, P, rec))
  end
  ret = xtuple(pr.blocks[end].branches[end].args[1], rec)
  ret = push!(pr, ret)
  pr.blocks[end].branches[end].args[1] = ret
  return pr, stks
end

function reverse_stacks!(adj, stks)
  ir = adj.adjoint
  entry = blocks(ir)[end]
  self = argument!(entry, at = 1)
  t = pushfirst!(blocks(ir)[end], xcall(:getfield, self, QuoteNode(:t)))
  repl = Dict()
  runonce(b) = b.id in (1, length(ir.blocks))
  for b in blocks(ir)
    for (i, (b′, α)) in enumerate(stks)
      b.id == b′ || continue
      if runonce(b)
        val = insertafter!(ir, t, xcall(:getindex, t, i))
      else
        stk = push!(entry, xcall(:getindex, t, i))
        stk = push!(entry, xcall(PartialP, :Stack, stk))
        val = pushfirst!(b, xcall(:pop!, stk))
      end
      repl[α] = val
    end
  end
  return IRTools.prewalk!(x -> get(repl, x, x), ir)
end

function stacks!(adj, T)
  forw, stks = forward_stacks!(adj, T)
  back = reverse_stacks!(adj, stks)
  permute!(back, length(back.blocks):-1:1)
  IRTools.domorder!(back)
  return forw, back
end

varargs(m::Method, n) = m.isva ? n - m.nargs + 1 : nothing

meta(T) = (usetyped ? IRTools.typed_meta : IRTools.meta)(T)

function getmeta(T)
  m = meta(T)
  (usetyped && m != nothing) || return m
  any(x -> isexpr(x, :goto, :gotoifnot), m.code.code) || return IRTools.meta(T)
  return m
end

function _lookup_grad(T)
  (m = getmeta(T)) == nothing && return
  m isa IRTools.TypedMeta && m.ret == Union{} && return
  va = varargs(m.method, length(T.parameters))
  forw, back = stacks!(Adjoint(IR(m), varargs = va, normalise = false), T)
  m, forw, back
end

function stacklines(T::Type)
  adj = Adjoint(IR(meta(T)), normalise = false)
  recs = []
  for b in blocks(adj.adjoint), α in alphauses(b)
    push!(recs, IRTools.exprline(adj.primal, Variable(α)))
  end
  return recs
end


text:
"""
Modia module with electric component models (inspired from Modelica Standard Library.

* Developer: Hilding Elmqvist, Mogram AB
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)

"""
module Electric

using ..Instantiation
using ..Blocks
using ..Synchronous: positive
using Unitful

using Modia

export Pin, Ground, OnePort, Resistor, Capacitor, Inductor,
  ConstantVoltage, StepVoltage, SignalVoltage, SineVoltage, IdealOpAmp3Pin, IdealDiode,
  Voltage, Current, Resistance, Capacitance

Voltage(; args...) = Variable(;T=Unitful.V, size=(), start=0.0, args...)
Current(; args...) = Variable(;T=Unitful.A, size=(), start=0.0, args...)
Resistance(; args...) = Variable(;T=Unitful.Ω, size=(), args...)
Capacitance(; args...) = Variable(;T=Unitful.F, size=(), args...)

@model Pin1 begin
  v=Float()
  i=Float(flow=true)
end

@model Pin begin
  v=Voltage()
  i=Current(flow=true)
end

@model Ground begin
  p=Pin()
@equations begin
  p.v = 0
  end
end

@model OnePort1 begin
  v=Float()
  i=Float()
  p=Pin()
  n=Pin()
@equations begin
  v = p.v - n.v
  0 = p.i + n.i
  i = p.i
  end
end

@model OnePort begin
#  v=Voltage()
#  i=Current()
  v=Voltage()
  i=Current()
  p=Pin()
  n=Pin()
@equations begin
  v = p.v - n.v
  0 = p.i + n.i
  i = p.i
  end
end

@model Resistor begin # Ideal linear electrical resistor
  @extends OnePort()
  @inherits i, v
  R=1 # Parameter(start=1.0) # undefined # Resistance
@equations begin
  R*i = v
  end
end

@model Resistor1 begin
  p=Pin()
  n=Pin()
  v=Float()
  i=Float()
  R=Parameter(description="Resistance")
@equations begin
  v = p.v - n.v # Voltage drop
  0 = p.i + n.i # Kirchhoff's current law within component
  i = p.i
  R*i = v
  end
end

# Setting state=false for v does not work with extends.
@model Capacitor begin
  @extends OnePort(v=Float(start=0.0))
  @inherits i, v
#  C=Capacitance() # undefined
  C=undefined
@equations begin
  C*der(v) = i
#  der(v) = i/C
  end
end

@model Capacitor1 begin
  p=Pin()
  n=Pin()
  v=Float()
  i=Float()
  C=undefined
@equations begin
  v = p.v - n.v # Voltage drop
  0 = p.i + n.i # Kirchhoff's current law within component
  i = p.i
  C*der(v) = i
  end
end

@model Inductor begin
    # Ideal linear electrical inductor
  @extends OnePort()
  @inherits i, v
  L=Parameter()   # Inductance
  @equations begin
  L*der(i) = v
  end
end

@model ConstantVoltage begin
  V=1u"V"
  @extends OnePort()
  @inherits v
@equations begin
  v = V
  end
end

@model StepVoltage begin
  V=1u"V"
  startTime = 0*Seconds
  t = Var(start=0.0)
  @extends OnePort()
  @inherits v
@equations begin
  v = if t < startTime; 0 else V end
  der(t) = 1
  end
end

@model SignalVoltage begin
    # Generic voltage source using the input signal as source voltage
  p=Pin()
  n=Pin()
  v=Float()
  i=Float()
@equations begin
  v = p.v - n.v
  0 = p.i + n.i
  i = p.i
  end
end

@model VoltageSource begin
    # Interface for voltage sources
  @extends OnePort()
  @inherits v
  offset=0 # Voltage offset
  startTime=0 # Time offset
  signalSource=SignalSource(offset = offset, startTime=startTime)
@equations begin
  v = signalSource.y
  end
end

@model SineVoltage1 begin
    # Sine voltage source
  V=Parameter() # Amplitude of sine wave
  phase=0 # Phase of sine wave
  freqHz=Parameter() # Frequency of sine wave
  @extends VoltageSource(signalSource = Sine(amplitude=V, freqHz=freqHz, phase=phase))
end

@model SineVoltage begin
    # Sine voltage source
  V=Parameter() # Amplitude of sine wave
  phase=0 # Phase of sine wave
  freqHz=Parameter() # Frequency of sine wave
#  @extends VoltageSource(signalSource = Sine())
  @extends VoltageSource(signalSource = Sine(amplitude=1,freqHz=1, offset=0, startTime=0))
#  @extends VoltageSource(signalSource = Sine(amplitude=V, freqHz=freqHz, phase=phase,
#    offset=offset, startTime=startTime))
#  @inherits offset, startTime
end

@model IdealOpAmp3Pin begin
    # Ideal operational amplifier (norator-nullator pair), but 3 pins
  in_p=Pin()
  in_n=Pin()
  out=Pin()
@equations begin
  in_p.v = in_n.v
  in_p.i = 0
  in_n.i = 0
  end
end

@model IdealDiode begin # Ideal diode
  @extends OnePort()
  @inherits v, i
  s = Float(start=0.0) # Auxiliary variable for actual position on the ideal diode characteristic
#=   s = 0: knee point
     s < 0: below knee point, diode conducting
     s > 0: above knee point, diode locking
=#
@equations begin
  v = if positive(s); 0 else s end
  i = if positive(s); s else 0 end
  end
end

@model IdealDiode2 begin # Ideal diode
  @extends OnePort()
  @inherits v, i
  Ron = 0 # 1.0E-5 # Forward state-on differential resistance (closed diode resistance)
  Goff = 0 # 1.0E-5 # Backward state-off conductance (opened diode conductance)
  Vknee = 0 # Forward threshold voltage
#  off = Variable(start=true) # Switching state
  s = Float(start=0.0) # Auxiliary variable for actual position on the ideal diode characteristic
#=  s = 0: knee point
     s < 0: below knee point, diode conducting
     s > 0: above knee point, diode locking
=#
@equations begin
#  off := s < 0
  v = s*if ! positive(s); 1 else Ron end + Vknee
  i = s*if ! positive(s); Goff else 1 end + Goff*Vknee
  end
end

@model IdealDiode1 begin # Ideal diode
  @extends OnePort()
  @inherits v, i
  Ron = 0 # 1.0E-5 # Forward state-on differential resistance (closed diode resistance)
  Goff = 0 # 1.0E-5 # Backward state-off conductance (opened diode conductance)
  Vknee = 0 # Forward threshold voltage
  off = Variable(start=true) # Switching state
  s = Float(start=0) # Auxiliary variable for actual position on the ideal diode characteristic
#=  s = 0: knee point
     s < 0: below knee point, diode conducting
     s > 0: above knee point, diode locking
=#
@equations begin
  off := positive(-s*1.0)
  v = s*if off; 1 else Ron end + Vknee
  i = s*if off; Goff else 1 end + Goff*Vknee
  end
end

end


text:
# ---
# title: 321. Create Maximum Number
# id: problem321
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/create-maximum-number/description/>
# hidden: true
# ---
#
# Given two arrays of length `m` and `n` with digits `0-9` representing two
# numbers. Create the maximum number of length `k <= m + n` from digits of the
# two. The relative order of the digits from the same array must be preserved.
# Return an array of the `k` digits.
#
# **Note:** You should try to optimize your time and space complexity.
#
# **Example 1:**
#
#
#
#     Input:
#     nums1 = [3, 4, 6, 5]
#     nums2 = [9, 1, 2, 5, 8, 3]
#     k = 5
#     Output:
#     [9, 8, 6, 5, 3]
#
# **Example 2:**
#
#
#
#     Input:
#     nums1 = [6, 7]
#     nums2 = [6, 0, 4]
#     k = 5
#     Output:
#     [6, 7, 6, 0, 4]
#
# **Example 3:**
#
#
#
#     Input:
#     nums1 = [3, 9]
#     nums2 = [8, 9]
#     k = 3
#     Output:
#     [9, 8, 9]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module VectorValuedMultiFEOperatorsTests

# Testing Stokes problem with Dirichlet BCs for the velocity
# and pressure prescribed in a single point

using Test
using Gridap
import Gridap: ∇

const T = VectorValue{2,Float64}

# Define manufactured functions
u1(x) = VectorValue(x[1], x[2])
u2(x) = x[1] - x[2]

∇u1(x) = TensorValue(1.0,0.0,0.0,1.0)

∇(::typeof(u1)) = ∇u1

b1(x) = VectorValue(1.0,-1.0)
b2(x) = 2.0

# Construct the discrete model
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))

# Construct the FEspace 1
order = 2
fespace1 = FESpace(
  reffe = :QLagrangian,
  conformity = :H1,
  valuetype = T,
  model = model,
  order = order,
  diritags = "boundary")

# Construct the FEspace 2
fespace2 = FESpace(
  reffe = :PLagrangian,
  conformity = :L2,
  valuetype = Float64,
  model = model,
  order = order-1,
  constraint = :zeromean)

# Define test and trial
V1 = TestFESpace(fespace1)
V2 = TestFESpace(fespace2)
V = [V1, V2]

U1 = TrialFESpace(fespace1,u1)
U2 = TrialFESpace(fespace2)
U = [U1, U2]

# Define integration mesh and quadrature for volume
trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

a(v,u) =
  inner(∇(v[1]),∇(u[1])) - inner(div(v[1]),u[2]) + inner(v[2],div(u[1]))
b(v) = inner(v[1],b1) + inner(v[2],b2)
t_Ω = AffineFETerm(a,b,trian,quad)

# Define the FEOperator
op = LinearFEOperator(V,U,t_Ω)

# Solve!
uh = solve(op)

# Define exact solution and error
e1 = u1 - uh[1]

e2 = u2 - uh[2]

#writevtk(trian,"trian",cellfields=["uh2"=>uh[2],"p"=>p])

# Define norms to measure the error
l2(u) = inner(u,u)
h1(u) = inner(∇(u),∇(u)) + l2(u)

# Compute errors
e1l2 = sqrt(sum( integrate(l2(e1),trian,quad) ))
e1h1 = sqrt(sum( integrate(h1(e1),trian,quad) ))

e2l2 = sqrt(sum( integrate(l2(e2),trian,quad) ))

@test e1l2 < 1.e-8
@test e1h1 < 1.e-8

@test e2l2 < 1.e-8

end # module


text:
@testset "log" begin

    V1 = [1; -2; 3; 4; 0; -6]
    V2 = -V1
    M = [1 1 0.5; 0.0 1.0 2.0]

    @testset "LogTransform" begin

        logV1 = [log(2); -log(3); log(4); log(5); log(1); -log(7)]
        logV2 = -logV1
        logM = [log(2) log(2) log(1.5); log(1) log(2) log(3)]

        @testset "simple" for x in (V1, V2, M)
            transform = LogTransform()
            @test cardinality(transform) == OneToOne()
            @test transform isa Transform
        end

        @testset "Apply" for (x, y) in ((V1, logV1), (V2, logV2), (M, logM))
            f = LogTransform()
            transformed = FeatureTransforms.apply(x, f)
            @test transformed ≈ y atol=1e-5
            @test FeatureTransforms.apply(transformed, f; inverse=true) ≈ x atol=1e-5
        end
    end

    @testset "InverseHyperbolicSine" begin

        logV1 = asinh.(V1)
        logV2 = asinh.(V2)
        logM = asinh.(M)

        @testset "simple" for x in (V1, V2, M)
            transform = InverseHyperbolicSine()
            @test cardinality(transform) == OneToOne()
            @test transform isa Transform
        end

        @testset "Apply" for (x, y) in ((V1, logV1), (V2, logV2), (M, logM))
            f = InverseHyperbolicSine()
            transformed = FeatureTransforms.apply(x, f)
            @test transformed ≈ y atol=1e-5
            @test FeatureTransforms.apply(transformed, f; inverse=true) ≈ x atol=1e-5
        end
    end
end


text:
using Parameters

struct DSGE_Model{data <: AbstractArray, alfa <: Float64}
    "Data"
    data::data
    "Alpha"
    alfa::alfa
end

function (problem::DSGE_Model)(pars)
    @unpack alfa, data = problem
    @unpack bet, epsilon, theta, sig, s2, phi, phi_pi, phi_y, rho_v = pars
    llh = log_like_dsge([alfa,bet,epsilon,theta,sig,s2,phi,phi_pi,phi_y,rho_v],data)
        return llh+logpdf(prior_bet,bet)+logpdf(prior_epsilon,epsilon)+logpdf(prior_theta,theta)+logpdf(prior_sig,sig)+logpdf(prior_s2,s2)+logpdf(prior_phi,phi)+logpdf(prior_phi_pi,phi_pi)+logpdf(prior_phi_y,phi_y)+logpdf(prior_rho_v,rho_v)
end

function problem_transform(p::DSGE_Model)
    as((bet = as_unit_interval, epsilon=as_positive_real,theta = as_unit_interval, sig = as_positive_real, s2 = as_positive_real, phi = as_positive_real, phi_pi = as_positive_real, phi_y = as_positive_real, rho_v = as_unit_interval))
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
module Arena

# due to issue with ImageMagick and Pkg importing, the order is very sensitive here!
# see https://github.com/JuliaIO/ImageMagick.jl/issues/142
using ImageMagick
using PlotUtils
using Caesar, ImageView, Images, MeshIO, MeshCat
using LinearAlgebra

using Rotations, CoordinateTransformations
using TransformUtils
using Graphs, NLsolve
using GeometryTypes, ColorTypes
using DocStringExtensions, ProgressMeter
# using CaesarLCMTypes
using Requires
using FileIO
using JSON
using Base64

const CTs = CoordinateTransformations
const TUs = TransformUtils

export
  meshgrid,
  DepthCamera,
  buildmesh!,
  reconstruct,
  VisualizationContainer,
  startDefaultVisualization,
  newtriad!,
  visualize,
  stopVis!,
  setGlobalDrawTransform!,
  visualizetriads,
  visualizeallposes!,
  visualizeDensityMesh!,
  updaterealtime!,
  visualizerealtime,
  # new tree interface
  drawpose!,
  drawposepoints!,
  drawLine!,
  drawLineBetween!,
  drawAllOdometryEdges!,
  pointToColor,
  findAllBinaryFactors,
  drawAllBinaryFactorEdges!,
  loadmodel,
  DrawModel,
  DrawROV,
  DrawScene,
  #deleting functions
  deletemeshes!,

  # more drawing utils
  ArcPointsRangeSolve,
  findaxiscenter!,
  parameterizeArcAffineMap,
  animatearc,

  # ImageUtils
  # image_tToRgb,
  rgbUint8ToRgb,
  rgbToJpeg,
  rgbToPng,
  imshowhackpng,
  cloudimshow,
  imshowhack,
  roi,

  # BotVis
  CameraModel,
  initBotVis2,
  drawPoses2!,
  cloudFromDepthImage,

  # point clouds
  visPointCloudOnPose!,
  drawPointCloudonPose!,
  cloudFromDepthImage,

  # colour gradients
  # re-exports
  cgrad,
  clibraries,
  cgradients

const NothingUnion{T} = Union{Nothing, T}

# Minimal globals
global loopvis = true
global drawtransform = Translation(0.0,0.0,0.0) ∘ LinearMap(Quat(1.0,0.0,0.0,0.0))

# types and models
include("Common/CameraModel.jl")
include("VisualizationTypes.jl")
include("RobotSceneModels.jl")

# Common
include("Common/DepthImages.jl")
# utils
include("GeneralUtils.jl")
include("ColorUtils.jl")
include("VisualizeLines.jl")
include("MeshUtils.jl")
include("BigDataUtils.jl")
include("ImageUtils.jl")
include("AnimationUtils.jl")
include("VisualizationUtils.jl")
include("VisualizePosesPoints.jl")
include("ModelVisualizationUtils.jl")
include("deprecated/Deprecated.jl")

# user interaction
include("HighLevelAPI.jl")
include("Amphitheatre/Amphitheatre.jl")
using .Amphitheatre

# plugins

# Developer tools
include("BotVis.jl")

# Used by Requires.jl to check if packages are imported. Much cleaner than janky isdefined().
function __init__()
  @info "Conditionally importing RoMEPlotting, GraffSDK, and Director..."
  # Checking what to import from the calling module
  @require GraffSDK="d47733cc-d211-5467-9efc-951b5b83f246" begin
    @info "--- GraffSDK is defined in the calling namespace, importing Graff functions..."
    include("plugins/GraffVisualizationService.jl")
    include("deprecated/DeprecatedGraff.jl")
    # Graff exports
    export visualizeSession
  end
  @require RoMEPlotting="238d586b-a4bf-555c-9891-eda6fc5e55a2" begin
    @info "--- RoMEPlotting is defined in the calling namespace, importing RoMEPlotting functions..."
  end
  @require DrakeVisualizer="49c7015b-b8db-5bc5-841b-fcb31c578176" begin
    @info "--- DrakeVisualizer is defined in the calling namespace, importing DrakeVisualizer functions..."
    include("DirectorVisService.jl")
    # DrakeVisualizer exports
    export drawdbdirector
    end
  end

# load the internal plugins that may or may not depend on the @requires above
include("plugins/VisualizationDefault.jl")
include("plugins/GetRobotConfiguration.jl")
include("plugins/ReprojectBearingRange.jl")
include("plugins/PointClouds.jl")

end # module


text:
## Currently, if/else will create its own local scope
## Assignments will not register in the loop's main scope
## although stores and return values will.
negateop!(ls::LoopSet, condop::Operation, elementbytes::Int) = add_compute!(ls, gensym!(ls, "negated#mask"), :~, [condop], elementbytes)

function add_if!(ls::LoopSet, LHS::Symbol, RHS::Expr, elementbytes::Int, position::Int, mpref::Union{Nothing,ArrayReferenceMetaPosition} = nothing)
    # for now, just simple 1-liners
    @assert length(RHS.args) == 3 "if statements without an else cannot be assigned to a variable."
    condition = first(RHS.args)
    condop = if condition isa Symbol
        getop(ls, condition, elementbytes)
    elseif isnothing(mpref)
        add_operation!(ls, gensym!(ls, "mask"), condition, elementbytes, position)
    else
        add_operation!(ls, gensym!(ls, "mask"), condition, mpref, elementbytes, position)
    end
    iftrue = RHS.args[2]
    if iftrue isa Expr
        trueop = add_operation!(ls, gensym!(ls, "iftrue"), iftrue, elementbytes, position)
        if iftrue.head === :ref && all(ld -> ld ∈ loopdependencies(trueop), loopdependencies(condop)) && !search_tree(parents(condop), trueop)
            trueop.instruction = Instruction(:conditionalload)
            push!(parents(trueop), condop)
        end
    else
        trueop = getop(ls, iftrue, elementbytes)
    end
    iffalse = RHS.args[3]
    if iffalse isa Expr
        falseop = add_operation!(ls, gensym!(ls, "iffalse"), iffalse, elementbytes, position)
        if iffalse.head === :ref && all(ld -> ld ∈ loopdependencies(falseop), loopdependencies(condop)) && !search_tree(parents(condop), falseop)
            falseop.instruction = Instruction(:conditionalload)
            push!(parents(falseop), negateop!(ls, condop, elementbytes))
        end
    else
        falseop = getop(ls, iffalse, elementbytes)
    end
    add_compute_ifelse!(ls, LHS, condop, trueop, falseop, elementbytes)
end

function add_andblock!(ls::LoopSet, condop::Operation, LHS, rhsop::Operation, elementbytes::Int, position::Int)
    if LHS isa Symbol
        altop = getop(ls, LHS, elementbytes)
        return add_compute_ifelse!(ls, LHS, condop, rhsop, altop, elementbytes)
    elseif LHS isa Expr && LHS.head === :ref
        return add_conditional_store!(ls, LHS, condop, rhsop, elementbytes)
    else
        throw("Don't know how to assign onto $LHS.")
    end
end
function add_andblock!(ls::LoopSet, condop::Operation, LHS, RHS::Expr, elementbytes::Int, position::Int)
    rhsop = add_compute!(ls, gensym!(ls, "iftruerhs"), RHS, elementbytes, position)
    add_andblock!(ls, condop, LHS, rhsop, elementbytes, position)
end
function add_andblock!(ls::LoopSet, condop::Operation, LHS, RHS, elementbytes::Int, position::Int)
    rhsop = getop(ls, RHS, elementbytes)
    add_andblock!(ls, condop, LHS, rhsop, elementbytes, position)
end
function add_andblock!(ls::LoopSet, condexpr::Expr, condeval::Expr, elementbytes::Int, position::Int)
    condop = add_operation!(ls, gensym!(ls, "mask"), condexpr, elementbytes, position)
    if condeval.head === :call
        @assert first(condeval.args) === :setindex!
        array, raw_indices = ref_from_setindex!(ls, condeval)
        ref = Expr(:ref, array); append!(ref.args, raw_indices)
        return add_andblock!(ls, condop, ref, condeval.args[3], elementbytes, position)
    end
    @assert condeval.head === :(=)
    @assert length(condeval.args) == 2
    LHS = condeval.args[1]
    RHS = condeval.args[2]
    add_andblock!(ls, condop, LHS, RHS, elementbytes, position)
end
function add_andblock!(ls::LoopSet, ex::Expr, elementbytes::Int, position::Int)
    add_andblock!(ls, first(ex.args)::Expr, last(ex.args)::Expr, elementbytes, position)
end

function add_orblock!(ls::LoopSet, condop::Operation, LHS, rhsop::Operation, elementbytes::Int, position::Int)
    negatedcondop = negateop!(ls, condop, elementbytes)
    if LHS isa Symbol
        altop = getop(ls, LHS, elementbytes)
        # return add_compute!(ls, LHS, :ifelse, [condop, altop, rhsop], elementbytes)
        # Placing altop second seems to let LLVM fuse operations; but as of LLVM 9.0.1 it will not if altop is first
        # therefore, we negate the condition and switch order so that the altop is second.
        return add_compute_ifelse!(ls, LHS, negatedcondop, rhsop, altop, elementbytes)
    elseif LHS isa Expr && LHS.head === :ref
        # negatedcondop = add_compute!(ls, gensym(:negated_mask), :~, [condop], elementbytes)
        return add_conditional_store!(ls, LHS, negatedcondop, rhsop, elementbytes)
    else
        throw("Don't know how to assign onto $LHS.")
    end
end
function add_orblock!(ls::LoopSet, condop::Operation, LHS, RHS::Expr, elementbytes::Int, position::Int)
    rhsop = add_compute!(ls, gensym!(ls, "iffalserhs"), RHS, elementbytes, position)
    add_orblock!(ls, condop, LHS, rhsop, elementbytes, position)
end
function add_orblock!(ls::LoopSet, condop::Operation, LHS, RHS, elementbytes::Int, position::Int)
    rhsop = getop(ls, RHS, elementbytes)
    add_orblock!(ls, condop, LHS, rhsop, elementbytes, position)
end
function add_orblock!(ls::LoopSet, condexpr::Expr, condeval::Expr, elementbytes::Int, position::Int)
    condop = add_operation!(ls, gensym!(ls, "mask"), condexpr, elementbytes, position)
    if condeval.head === :call
        @assert first(condeval.args) === :setindex!
        array, raw_indices = ref_from_setindex!(ls, condeval)
        ref = Expr(:ref, array); append!(ref.args, raw_indices)
        return add_orblock!(ls, condop, ref, condeval.args[3], elementbytes, position)
    end
    @assert condeval.head === :(=)
    @assert length(condeval.args) == 2
    LHS = condeval.args[1]
    RHS = condeval.args[2]
    add_orblock!(ls, condop, LHS, RHS, elementbytes, position)
end
function add_orblock!(ls::LoopSet, ex::Expr, elementbytes::Int, position::Int)
    add_orblock!(ls, first(ex.args)::Expr, last(ex.args)::Expr, elementbytes, position)
end