text:
using Definitions
using DataFrames
using CSVFiles
using Model_Household

function map_person( model_person :: DataFrameRow )

    income = Incomes_Dict()
    for i in instances(Incomes_Type)
        ikey = make_sym_for_frame("income", i)
        if ! ismissing(model_person[ikey])
            if model_person[ikey] != 0.0
                v = model_person[ikey] # this is a hack because in j 1.3 csv is parsing many cols as strings & I don't understand why
                tv = typeof(v)
                if tv == String # FIXME delete not needed
                    v = parse( Float64, v )
                end
                # println( "setting ikey = $ikey to $v of type $tv")
                income[i] = v
            end
        end
    end

    assets = Asset_Dict()
    for i in instances(Asset_Type)
        if i != Missing_Asset_Type
            ikey = make_sym_for_asset( i )
            if ! ismissing(model_person[ikey])
                v = model_person[ikey]
                if typeof(v) == String
                    v = parse( Float64, v )
                end
                if model_person[ikey] != 0.0
                    assets[i] = model_person[ikey]
                end
            end
        end
    end

    disabilities = Disability_Dict()
    for i in instances(Disability_Type)
        ikey = make_sym_for_frame("disability", i)
        if ! ismissing(model_person[ikey])
            if model_person[ikey] == 1
                disabilities[i] = Bool(model_person[ikey])
            end
        end
    end

    relationships = Relationship_Dict()
    for i in 1:14
        relmod = Symbol( "relationship_$(i)") # :relationship_10 or :relationship_2
        irel = model_person[relmod]
        if (! ismissing( irel )) & ( irel >= 0 )
            pid = get_pid(
                FRS,
                model_person.data_year,
                model_person.hid,
                i )
            relationships[pid] = Relationship( irel )
        end
    end

    Person(

        BigInt(model_person.hid),  # BigInt# == sernum
        BigInt(model_person.pid),  # BigInt# == unique id (year * 100000)+
        model_person.pno,  # Integer# person number in household
        model_person.default_benefit_unit,  # Integer
        model_person.age,  # Integer

        Sex(model_person.sex),  # Sex
        Ethnic_Group(safe_assign(model_person.ethnic_group)),  # Ethnic_Group
        Marital_Status(safe_assign(model_person.marital_status)),  # Marital_Status
        Qualification_Type(safe_assign(model_person.highest_qualification)),  # Qualification_Type

        SIC_2007(safe_assign(model_person.sic)),  # SIC_2007
        Standard_Occupational_Classification(safe_assign(model_person.occupational_classification)),  # Standard_Occupational_Classification
        Employment_Sector(safe_assign(model_person.public_or_private)),  #  Employment_Sector
        Employment_Type(safe_assign(model_person.principal_employment_type)),  #  Employment_Type

        Socio_Economic_Group(safe_assign(model_person.socio_economic_grouping)),  # Socio_Economic_Group
        m2z(model_person.age_completed_full_time_education),  # Integer
        m2z(model_person.years_in_full_time_work),  # Integer
        ILO_Employment(safe_assign(model_person.employment_status)),  # ILO_Employment
        m2z(model_person.actual_hours_worked),  # Real
        m2z(model_person.usual_hours_worked),  # Real

        income,
        assets,

        safe_to_bool(model_person.registered_blind),
        safe_to_bool(model_person.registered_partially_sighted),
        safe_to_bool(model_person.registered_deaf),

        disabilities,
        Health_Status(safe_assign(model_person.health_status)),
        relationships,
        Relationship(model_person.relationship_to_hoh),
        safe_to_bool(model_person.is_informal_carer),
        safe_to_bool(model_person.receives_informal_care_from_non_householder),
        m2z(model_person.hours_of_care_received),
        m2z(model_person.hours_of_care_given),
        m2z(model_person.hours_of_childcare),
        m2z(model_person.cost_of_childcare),
        Child_Care_Type(safe_assign(model_person.childcare_type )),
        safe_to_bool( model_person.employer_provides_child_care ),

        Fuel_Type( m2z(model_person.company_car_fuel_type )),
        m2z(model_person.company_car_value),
        m2z(model_person.company_car_contribution),
        m2z(model_person.fuel_supplied)

    )

end

function map_hhld( hno::Integer, frs_hh :: DataFrameRow )

    people = People_Dict()

    Household(
        hno,
        frs_hh.hid,
        frs_hh.interview_year,
        frs_hh.interview_month,
        frs_hh.quarter,
        Tenure_Type(frs_hh.tenure),
        Standard_Region(frs_hh.region),
        CT_Band(frs_hh.ct_band),
        m2z(frs_hh.council_tax),
        m2z(frs_hh.water_and_sewerage ),
        m2z(frs_hh.mortgage_payment),
        m2z(frs_hh.mortgage_interest),
        m2z(frs_hh.years_outstanding_on_mortgage),
        m2z(frs_hh.mortgage_outstanding),
        m2z(frs_hh.year_house_bought),
        m2z(frs_hh.gross_rent),
        safe_to_bool(frs_hh.rent_includes_water_and_sewerage),
        m2z(frs_hh.other_housing_charges),
        m2z(frs_hh.gross_housing_costs),
        m2z(frs_hh.total_income),
        m2z(frs_hh.total_wealth),
        m2z(frs_hh.house_value),
        m2z(frs_hh.weight),
        people )
end

function load_hhld_from_frame( year :: Integer, hid :: Integer; hhld_fr :: DataFrame, pers_fr :: DataFrame ) :: Household
     frs_hh = hhld_fr[((hhld_fr.data_year .== year).& (hhld_fr.hid .== hid)),:]
     nhh = size( frs_hh )[1]
     @assert nhh in [0,1]
     if nhh == 1
         return load_hhld_from_frs( 1, frs_hh[1,:], pers_fr )
     else
        return missing
     end
end

function load_hhld_from_frame( hseq::Integer, hhld_fr :: DataFrameRow, pers_fr :: DataFrame ) :: Household
     hh = map_hhld( hseq, hhld_fr )
     pers_fr_in_this_hh = pers_fr[((pers_fr.data_year .== hhld_fr.data_year).&(pers_fr.hid .== hh.hid)),:]
     npers = size( pers_fr_in_this_hh )[1]
     @assert npers in 1:19
     for p in 1:npers
         pers = map_person( pers_fr_in_this_hh[p,:])
         hh.people[pers.pid] = pers
     end
     hh
 end


text:
using IterTools

cd(dirname(@__FILE__))
const lines = open("four.txt", "r") do f
    readlines(f)
end

parse_line(l::String)::Array{String} =
    [split(s,":")[1] for s in split(l)]

join_by_space(l) = join(l, " ")

valid_set = Set(["hcl", "iyr", "eyr", "ecl", "pid", "byr", "hgt"])
is_valid(s) = issubset(valid_set, s)

n = (groupby(!isempty, lines)
    |> l -> map(join_by_space, l)
    |> l -> filter(!isempty, l)
    |> l -> map(parse_line, l)
    |> s -> map(is_valid, s)
    |> sum)

println(n)

ext: jl


text:
"""
    disableCustomHandler()

Disables the user custom error handler.

"""
function disableCustomHandler()
    ccall((:calceph_seterrorhandler, libcalceph), Cvoid,
                  (Cint, Ptr{Cvoid}), 1, C_NULL)
end

mutable struct UserHandlerContainer
    f::Function
end

const userHandlerContainerInstance = UserHandlerContainer(s::String->Nothing)

function userHandlerWrapper(msg::Cstring)::Cvoid
    global userHandlerContainerInstance
    s = unsafe_string(msg)
    userHandlerContainerInstance.f(s)
    return
end

# see https://discourse.julialang.org/t/cfunction-error-handler/20678
#userHandlerCWrapper = @cfunction(userHandlerWrapper, Cvoid, (Cstring,))
userHandlerCWrapper = Nothing

"""
    setCustomHandler(f::Function)

Sets the user custom error handler.

# Arguments
- `f`: function taking a single argument of type String which will contain the CALCEPH error message. f should return Nothing.

Use setCustomHandler(s->Nothing) to disable CALCEPH error messages printout to the console.

"""
function setCustomHandler(f::Function)
    global userHandlerContainerInstance
    global userHandlerCWrapper
    userHandlerCWrapper = @cfunction(userHandlerWrapper, Cvoid, (Cstring,))
    userHandlerContainerInstance.f = f
    ccall((:calceph_seterrorhandler, libcalceph), Cvoid,
                  (Cint, Ptr{Cvoid}), 3, userHandlerCWrapper)
end


text:
module FoldsPreviewTests

using Distributed
using LoadAllPackages
using Test

include("utils.jl")

function include_tests(m = @__MODULE__, dir = @__DIR__)
    for file in readdir(dir)
        if match(r"^test_.*\.jl$", file) !== nothing
            Base.include(m, joinpath(dir, file))
        end
    end
end

include_tests()

function collect_modules(root::Module)
    modules = Module[]
    for n in names(root, all = true)
        m = getproperty(root, n)
        m isa Module || continue
        m === root && continue
        startswith(string(nameof(m)), "Test") || continue
        push!(modules, m)
    end
    return modules
end

collect_modules() = collect_modules(@__MODULE__)

this_project() = joinpath(dirname(@__DIR__), "Project.toml")

function is_in_path()
    project = this_project()
    paths = Base.load_path()
    project in paths && return true
    realproject = realpath(project)
    realproject in paths && return true
    matches(path) = path == project || path == realproject
    return any(paths) do path
        matches(path) || matches(realpath(path))
    end
end

function with_project(f)
    is_in_path() && return f()
    load_path = copy(LOAD_PATH)
    push!(LOAD_PATH, this_project())
    try
        f()
    finally
        append!(empty!(LOAD_PATH), load_path)
    end
end

function load_me_everywhere()
    prepare_impl()
    @everywhere append!(empty!(LOAD_PATH), $(copy(LOAD_PATH)))
    pkgid = Base.PkgId(@__MODULE__)
    @everywhere Base.require($pkgid)
    @everywhere $prepare_impl()
end

function prepare_impl()
    LoadAllPackages.loadall(this_project())
end

function runtests(modules = collect_modules())
    # `with_project` and `LoadAllPackages.loadall` are for doctests
    with_project() do
        if get(ENV, "CI", "false") == "true"
            if nprocs() < 4
                addprocs(4 - nprocs())
            end
        end
        @info "Testing with:" nprocs()

        load_me_everywhere()
        runtests_impl(modules)
    end
end

function runtests_impl(modules)
    @testset "$(nameof(m))" for m in modules
        tests = map(names(m, all = true)) do n
            n == :test || startswith(string(n), "test_") || return nothing
            f = getproperty(m, n)
            f !== m || return nothing
            parentmodule(f) === m || return nothing
            applicable(f) || return nothing  # removed by Revise?
            return f
        end
        filter!(!isnothing, tests)
        @testset "$f" for f in tests
            @debug "Testing $m.$f"
            f()
        end
    end
end

end # module


text:
# # Box-Muller method to Generate normal distribution
using NiLang

# In this tutorial, we introduce using Box-Muller method to transform a uniform distribution to a normal distribution.
# The transformation and inverse transformation of `Box-Muller` method could be found in
# [this blog](https://mathworld.wolfram.com/Box-MullerTransformation.html)
@i function boxmuller(x::T, y::T) where T
    @routine @invcheckoff begin
        @zeros T θ logx _2logx
        θ += 2π * y
        logx += log(x)
        _2logx += -2 * logx
    end

    ## store results
    z1 ← zero(T)
    z2 ← zero(T)
    z1 += _2logx ^ 0.5
    ROT(z1, z2, θ)
    ~@routine

    SWAP(x, z1)
    SWAP(y, z2)

    ## arithmetic uncomputing: recomputing the original values of `x` and `y` to deallocate z1 and z2
    @routine @invcheckoff begin
        @zeros T at sq _halfsq
        at += atan(y, x)
        if (y < 0, ~)
            at += T(2π)
        end
        sq += x ^ 2
        sq += y ^ 2
        _halfsq -= sq / 2
    end
    z1 -= exp(_halfsq)
    z2 -= at / (2π)
    @invcheckoff z1 → zero(T)
    @invcheckoff z2 → zero(T)
    ~@routine
end

# One may wonder why this implementation is so long,
# should't NiLang generate the inverse for user?
# The fact is, although Box-Muller is arithmetically reversible.
# It is not finite precision reversible.
# Hence we need to "uncompute" it manually,
# this trick may introduce reversibility error.

using Plots
N = 5000
x = rand(2*N)

Plots.histogram(x, bins = -3:0.1:3, label="uniform",
    legendfontsize=16, xtickfontsize=16, ytickfontsize=16)

# forward
@instr boxmuller.(x[1:N], x[N+1:end])
Plots.histogram(x, bins = -3:0.1:3, label="normal",
    legendfontsize=16, xtickfontsize=16, ytickfontsize=16)

# backward
@instr (~boxmuller).(x[1:N], x[N+1:end])
Plots.histogram(x, bins = -3:0.1:3, label="uniform",
    legendfontsize=16, xtickfontsize=16, ytickfontsize=16)

# ## Check the probability distribution function
using LinearAlgebra, Test

normalpdf(x) = sqrt(1/2π)*exp(-x^2/2)

# obtain `log(abs(det(jacobians)))`
@i function f(x::Vector)
    boxmuller(x[1], x[2])
end
jac = NiLang.AD.jacobian(f, [0.5, 0.5], iin=1)
ladj = log(abs(det(jac)))

# check if it matches the `log(p/q)`.
z1, z2 = boxmuller(0.5, 0.5)
@test ladj ≈ log(1.0 / (normalpdf(z1) * normalpdf(z2)))

# ## To obtaining Jacobian - a simpler approach
# We can define a function that exactly reversible from the instruction level,
# but costs more space for storing output.
@i function boxmuller2(x1::T, x2::T, z1::T, z2::T) where T
    @routine @invcheckoff begin
        @zeros T θ logx _2logx

        θ += 2π * x2
        logx += log(x1)
        _2logx += -2 * logx
    end

    ## store results
    z1 += _2logx ^ 0.5
    ROT(z1, z2, θ)
    ~@routine
end

# However, this is not a bijector from that maps `x` to `z`,
# because computing the backward just erases the content in `z`.
# However, this function can be used to obtain `log(abs(det(jacobians)))`
@i function f2(x::Vector, z::Vector)
    boxmuller2(x[1], x[2], z[1], z[2])
end
jac = NiLang.AD.jacobian(f2, [0.5, 0.5], [0.0, 0.0], iin=1, iout=2)
ladj = log(abs(det(jac)))

# check if it matches the `log(p/q)`.
_, _, z1, z2 = boxmuller2(0.5, 0.5, 0.0, 0.0)
@test ladj ≈ log(1.0 / (normalpdf(z1) * normalpdf(z2)))


text:
using RNASplicing
Lib = RNASplicing.Library
using FASTX
using BioSequences

# Generate 1000-length RNA with 4% chance of N, 24% for A, C, G, or U
sp = SamplerWeighted(rna"ACGU", fill(0.25, 3))
seq = randseq(RNAAlphabet{4}(), sp, 50)
splice(seq, Lib.majorsos, Lib.coremachinery)


text:
using Test
using Rectangle
using Rectangle: isnil

const VISUALIZATION=false

function xy_predicate(v, sdir)
    dir  = sdir == 1 ? 2 : 1
    odir = dir == 1 ? 2 : 1

    sz = size(v)

    vloc = dir == 1 ? div(sz[dir], 4) : div(sz[dir], 2)

    vloc <= 8 && return nothing, nothing, XYData(0, 0, (0, 0)), 0

    prange = parentindices(v)
    loc = vloc + prange[dir][1] - 1
    range = (prange[odir][1] + firstindex(v, odir) -1,
             prange[odir][1] + lastindex(v, odir)  -1)

    nvl, nvr = dir == 1 ? ((@view v[1:vloc, :]), (@view v[vloc+1:end, :])) :
        ((@view v[:, 1:vloc]), (@view v[:, vloc+1:end]))

    return (()->xy_predicate(nvl, dir),
            ()->xy_predicate(nvr, dir),
            XYData(dir, loc, range), 0)
end

@testset "XY Tree Tests" begin
    m = falses(64, 64)

    t = XYTree{Int, Int}() do
        xy_predicate(m, 0)
    end

    ls, vs = get_values(t)

    @test ls == [Line(32, 1, 32, 16),
                 Line(1, 16, 64, 16),
                 Line(16, 17, 16, 28),
                 Line(1, 28, 32, 28),
                 Line(1, 37, 16, 37),
                 Line(16, 29, 16, 64),
                 Line(17, 37, 32, 37),
                 Line(32, 17, 32, 64),
                 Line(48, 17, 48, 28),
                 Line(33, 28, 64, 28),
                 Line(33, 37, 48, 37),
                 Line(48, 29, 48, 64),
                 Line(49, 37, 64, 37)]
    @test vs == zeros(Int, 13)
end

@static if VISUALIZATION

    import PyPlot
    using PyCall

    function plot_map(bimg, t)
        fig = PyPlot.gcf()
        PyPlot.clf()
        fig[:clf](true)
        fig[:set_dpi](72f0)

        sz = size(bimg)
        hpi = sz[1]/16f0
        wpi = sz[2]/16f0
        println(hpi, " ", wpi)
        fig[:set_figwidth](wpi, forward=true)
        fig[:set_figheight](hpi, forward=true)

        ax = fig[:gca]()

        ax[:set_xlim](0f0, sz[2])
        ax[:set_ylim](0f0, sz[1])

        img = map(x->x ? 0xff : 0x00, bimg)
        ls, vs = get_values(t)
        for l in ls
            img[ly(l):ry(l), lx(l):rx(l)] .= 0x7f
        end
        imgplot = PyPlot.imshow(img)
        PyPlot.show()
    end
end


text:
"""
    Trim{N} <: AbstractModel

Trimmed aircraft controller model with state variables corresponding to control
inputs which are used to trim the aircraft and inputs corresponding to the total
forces and moments on the aircraft ``F_x, F_y, F_z, M_x, M_y, M_z``.
"""
struct Trim{N} <: AbstractModel
    state_indices::NTuple{N,Int}
    force_indices::NTuple{N,Int}
end

"""
    Trim(state_indices::NTuple{N}, force_indices::NTuple{N})

Initialize an object of type [`Trim`](@ref). Use the state variables
corresponding to `state_indices` to trim the forces/moments corresponding to
`force_indices`.
"""
Trim(state_indices, force_indices)

Trim() = Trim((1,2,3,4,5,6), (1,2,3,4,5,6))

# --- Traits --- #

number_of_states(::Type{Trim{N}}) where N = N
number_of_inputs(::Type{Trim{N}}) where N = N
number_of_parameters(::Type{<:Trim}) = 0

inplaceness(::Type{<:Trim}) = OutOfPlace()

rate_jacobian_type(::Type{<:Trim}) = Zeros()
state_jacobian_type(::Type{<:Trim}) = Zeros()
input_jacobian_type(::Type{<:Trim}) = Identity()
parameter_jacobian_type(::Type{<:Trim}) = Empty()
time_gradient_type(::Type{<:Trim}) = Zeros()

# --- Methods --- #

get_residual(::Trim, dx, x, y, p, t) = y

# --- Performance Overloads --- #

# NOTE: No performance overloads are needed since the jacobians are fully
# defined by their provided types

# --- Convenience Methods --- #

set_states!(x, model::Trim; delta) = x .= delta

function set_inputs!(y, model::Trim; Fx, Fy, Fz, Mx, My, Mz)

    y[1] = Fx
    y[2] = Fy
    y[3] = Fz
    y[4] = Mx
    y[5] = My
    y[6] = Mz

    return y
end

set_parameters!(p, model::Trim) = p

separate_states(model::Trim, x) = (delta = x,)

function separate_inputs(model::Trim, y)

    return (Fx = y[1], Fy = y[2], Fz = y[3], Mx = y[4], My = y[5], Mz = y[6])
end

separate_parameters(model::Trim, p) = ()

# --- Internal Methods for Couplings --- #


text:
export BrickletVoltageCurrentV2CurrentCallbackConfiguration
struct BrickletVoltageCurrentV2CurrentCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletVoltageCurrentV2VoltageCallbackConfiguration
struct BrickletVoltageCurrentV2VoltageCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletVoltageCurrentV2PowerCallbackConfiguration
struct BrickletVoltageCurrentV2PowerCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletVoltageCurrentV2Configuration
struct BrickletVoltageCurrentV2Configuration
    averaging::Integer
    voltage_conversion_time::Integer
    current_conversion_time::Integer
end

export BrickletVoltageCurrentV2Calibration
struct BrickletVoltageCurrentV2Calibration
    voltage_multiplier::Integer
    voltage_divisor::Integer
    current_multiplier::Integer
    current_divisor::Integer
end

export BrickletVoltageCurrentV2SPITFPErrorCount
struct BrickletVoltageCurrentV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletVoltageCurrentV2Identity
struct BrickletVoltageCurrentV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletVoltageCurrentV2
"""
Measures power, DC voltage and DC current up to 720W/36V/20A
"""
mutable struct BrickletVoltageCurrentV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletVoltageCurrentV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_voltage_current_v2")
        deviceInternal = package.BrickletVoltageCurrentV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_current
"""
    $(SIGNATURES)

Returns the current.

If you want to get the value periodically, it is recommended to use the
:cb:`Current` callback. You can set the callback configuration
with :func:`Set Current Callback Configuration`.
"""
function get_current(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_current()
end

export set_current_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Current` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Current` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_current_callback_configuration(device::BrickletVoltageCurrentV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_current_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_current_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Current Callback Configuration`.
"""
function get_current_callback_configuration(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_current_callback_configuration()
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage.

If you want to get the value periodically, it is recommended to use the
:cb:`Voltage` callback. You can set the callback configuration
with :func:`Set Voltage Callback Configuration`.
"""
function get_voltage(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_voltage()
end

export set_voltage_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Voltage` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Voltage` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_voltage_callback_configuration(device::BrickletVoltageCurrentV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_voltage_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_voltage_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Voltage Callback Configuration`.
"""
function get_voltage_callback_configuration(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_voltage_callback_configuration()
end

export get_power
"""
    $(SIGNATURES)

Returns the power.

If you want to get the value periodically, it is recommended to use the
:cb:`Power` callback. You can set the callback configuration
with :func:`Set Power Callback Configuration`.
"""
function get_power(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_power()
end

export set_power_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Power` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Power` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_power_callback_configuration(device::BrickletVoltageCurrentV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_power_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_power_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Power Callback Configuration`.
"""
function get_power_callback_configuration(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_power_callback_configuration()
end

export set_configuration
"""
    $(SIGNATURES)

Sets the configuration of the Voltage/Current Bricklet 2.0. It is
possible to configure number of averages as well as
voltage and current conversion time.
"""
function set_configuration(device::BrickletVoltageCurrentV2, averaging, voltage_conversion_time, current_conversion_time)
    device.deviceInternal.set_configuration(averaging, voltage_conversion_time, current_conversion_time)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_configuration()
end

export set_calibration
"""
    $(SIGNATURES)

Since the ADC and the shunt resistor used for the measurements
are not perfect they need to be calibrated by a multiplier and
a divisor if a very precise reading is needed.

For example, if you are expecting a current of 1000mA and you
are measuring 1023mA, you can calibrate the Voltage/Current Bricklet
by setting the current multiplier to 1000 and the divisor to 1023.
The same applies for the voltage.

The calibration will be saved on the EEPROM of the Voltage/Current
Bricklet and only needs to be done once.
"""
function set_calibration(device::BrickletVoltageCurrentV2, voltage_multiplier, voltage_divisor, current_multiplier, current_divisor)
    device.deviceInternal.set_calibration(voltage_multiplier, voltage_divisor, current_multiplier, current_divisor)
end

export get_calibration
"""
    $(SIGNATURES)

Returns the calibration as set by :func:`Set Calibration`.
"""
function get_calibration(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_calibration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletVoltageCurrentV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletVoltageCurrentV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletVoltageCurrentV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletVoltageCurrentV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletVoltageCurrentV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletVoltageCurrentV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletVoltageCurrentV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletVoltageCurrentV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""

Linear-Stability-Calculators
============================
Repo:      Linear-Stability-Calculators
Code:      ShallowWater/Julia/linear_stability_shallow_water.jl
Model:     Rotating Shallow Water
Geometry:  Cartesian/Spherical
Structure: Bickley Jet

Nondimensional Parameters
=========================
Fr = Froude Number
Ro = Rossby Number

Created Feb 26, 2021
By Francis J. Poulin

"""

abstract type AbstractGeometry end

struct Cartesian <: AbstractGeometry end
struct Spherical <: AbstractGeometry end

include("Grid.jl")
include("Physics.jl")
include("Basic_States.jl")
include("Plotting_Functions.jl")
include("Compute_Eigenvalues.jl")


text:
# # Linear regression

import MaximumLikelihoodProblems

import Distributions
import ForwardDiff
import LogDensityProblems
import TransformVariables

struct LinearRegression{Ty, TX}
    y::Ty
    X::TX
end

function (problem::LinearRegression)(θ)
    y = problem.y
    X = problem.X

    β = θ.β
    σ = θ.σ

    η = X * β
    μ = η
    ε = y - μ

    ## these two lines are equivalent:
    ## log_likelihood = Distributions.loglikelihood(Distributions.Normal(0, σ), ε)
    ## log_likelihood = sum(Distributions.logpdf.(Distributions.Normal(0, σ), ε))

    log_likelihood = sum(Distributions.logpdf.(Distributions.Normal(0, σ), ε))
    return log_likelihood
end

N = 10_000

## the first column (the column of all ones) is the intercept
X = hcat(ones(N), randn(N), randn(N))

size_β = (3,)
β_true = [1.0, 2.0, -1.0]
σ_true = 0.5
η_true = X * β_true
μ_true = η_true
ε_true = randn(N) .* σ_true
y = μ_true + ε_true

function generate_problem_transformation(p::LinearRegression)
    return TransformVariables.as((β = TransformVariables.as(Array, size(p.X, 2)),
                                  σ = TransformVariables.asℝ₊))
end

problem = LinearRegression(y, X)
transformation = generate_problem_transformation(problem)
transformed_problem = LogDensityProblems.TransformedLogDensity(transformation,
                                                               problem)
transformed_gradient_problem = LogDensityProblems.ADgradient(:ForwardDiff,
                                                             transformed_problem)

β_hat_initial_guess = zeros(size_β)
σ_hat_initial_guess = 1.0

θ_hat_initial = (; β = β_hat_initial_guess,
                   σ = σ_hat_initial_guess)

# θ_hat:

θ_hat = MaximumLikelihoodProblems.fit(transformed_gradient_problem,
                                      θ_hat_initial;
                                      learning_rate = 1e-6)

# β_hat:

β_hat = θ_hat[:β]

# σ_hat:

σ_hat = θ_hat[:σ]

# Value of the log likelihood function evaluated at θ_hat:

MaximumLikelihoodProblems.loglikelihood(transformed_gradient_problem, θ_hat)


text:
using LinearSolve
using Test

@testset "LinearSolve.jl" begin
    A = rand(5, 5)
    b = rand(5)
    prob = LinearProblem(A, b)

    # Factorization
    @test A * solve(prob, LUFactorization(); alias_A = false, alias_b = false) ≈
          b
    @test A * solve(prob, QRFactorization(); alias_A = false, alias_b = false) ≈
          b
    @test A *
          solve(prob, SVDFactorization(); alias_A = false, alias_b = false) ≈ b

    # Krylov
    @test A * solve(prob, KrylovJL(A, b)) ≈ b
end


text:
# This file was generated, do not modify it. # hide
using MLJ, RDatasets, PrettyPrinting
MLJ.color_off() # hide
import Distributions
const D = Distributions

@load LinearRegressor pkg=MLJLinearModels
@load RidgeRegressor pkg=MLJLinearModels
@load LassoRegressor pkg=MLJLinearModels

hitters = dataset("ISLR", "Hitters")
@show size(hitters)
names(hitters) |> pprint

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
ENV["GKSwstype"] = "100"
using Documenter, GPGPU

DocMeta.setdocmeta!(GPGPU, :DocTestSetup,
                   :(using GPGPU); recursive=true)

# generate Literate documentation
#include("generate.jl")

makedocs(
    format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS"),
                             collapselevel = 1,
                             assets = ["assets/juliareach-light.css"]),
    sitename = "GPGPU",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "Assignments" => Any["Affine cipher" => "man/cipher.md"]]
)

deploydocs(
    repo = "github.com/mforets/GPGPU",
    push_preview = true,
)

max_stars_repo_name: mforets/GPGPU

max_issues_repo_name: mforets/GPGPU

max_forks_repo_name: mforets/GPGPU


text:
# Run with
#   include("learn_struct.jl")
# not with
#   Execute File
module Allnix
mutable struct Well
    Well() = (
        me = new(); # new Well object
        # Assign default values
        #me.well_name = nothing;

        me # return the new object
    )

    well_name::String
end

mutable struct Simulation
    Simulation() = (
        me = new(); # new Simulation object
        # Assign default values
        me.well_list = Vector{Well}();

        me # return the new object
    )
    well_list::Vector{Well}
end

export Well, Simulation

end

sim = Allnix.Simulation()
push!(sim.well_list, Allnix.Well())


text:
using CmeSIM
using Documenter

DocMeta.setdocmeta!(CmeSIM, :DocTestSetup, :(using CmeSIM); recursive=true)

makedocs(;
    modules=[CmeSIM],
    authors="Martin Kunz <martinkunz@email.cz> and contributors",
    repo="https://github.com/kunzaatko/CmeSIM.jl/blob/{commit}{path}#{line}",
    sitename="CmeSIM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://kunzaatko.github.io/CmeSIM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/kunzaatko/CmeSIM.jl.git",
    devbranch="trunk"
)


text:
# # Evolution via method of lines

using FFTW
using Plots
using OrdinaryDiffEq
using LinearAlgebra: norm, mul!
using Test
using BenchmarkTools
using Random

@info "Threads: $(FFTW.nthreads())"

# ## Operators

function get_operators(N, κ₀)
    ## Differentiation in spectral space
    Dx_hat = im * κ₀ * [
        ifelse(k1 ≤ div(N, 2) + 1, k1 - 1, k1 - 1 - N) for k2 in 1:div(N, 2)+1, k1 in 1:N
    ]
    Dy_hat = im * κ₀ * [k2 - 1 for k2 in 1:div(N, 2)+1, k1 in 1:N]

    Delta_hat = - κ₀^2 * [
        ifelse(k1 ≤ div(N, 2) + 1, (k1 - 1)^2 + (k2 - 1)^2, (k1 - 1 - N)^2 + (k2 - 1)^2)
        for k2 in 1:div(N, 2)+1, k1 in 1:N
    ]

    ## For the Basdevant formulation
    DxsqDysq_hat = Dx_hat.^2 .- Dy_hat.^2
    Dxy_hat = Dx_hat .* Dy_hat

    ## Recovering of the velocity field from the vorticity
    Hu_hat = - Dy_hat ./ Delta_hat
    Hu_hat[1, 1] = 0.0
    Hv_hat = Dx_hat ./ Delta_hat
    Hv_hat[1, 1] = 0.0

    return Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat
end

# Methods to generate a scalar field from a list of wavenumbers and amplitudes

function field_from_modes(L, N, modes::Matrix{<:Integer}, amps::Matrix{<:Real})
    κ₀ = 2π/L
    x = y = (L/N):(L/N):L
    field = sum(
        [
            2κ₀^2 * (k[1]^2 + k[2]^2) * (
                a[1] * cos.(κ₀ * (k[1] * one.(y) * x' + k[2] * y * one.(x)'))
                - a[2] * sin.(κ₀ * (k[1] * one.(y) * x' + k[2] * y * one.(x)'))
            )
            for (k, a) in zip(eachrow(modes), eachrow(amps))
        ]
    )
    return field
end

function field_from_modes(rng::AbstractRNG, L, N, num_modes::Int)

    modes = rand(rng, 1:div(N,10), num_modes, 2)
    amps = rand(rng, num_modes, 2)

    field = field_from_modes(L, N, modes, amps)

    return field
end

field_from_modes(L, N, num_modes::Int) = field_from_modes(Xoshiro(), L, N, num_modes)

# Differential equations

function nsepervort_hat_rhs!(dvorhattdt, vort_hat, params, t)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    ν, g_hat, N, Nsub = vars
    u_hat, v_hat, u, v, uv, v2u2, uv_hat, v2u2_hat = auxs
    plan, plan_inv = plans

    vort_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    vort_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im

    u_hat .= Hu_hat .* vort_hat
    v_hat .= Hv_hat .* vort_hat

    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)

    uv .= u .* v
    v2u2 .= v.^2 .- u.^2
    mul!(uv_hat, plan, uv)
    mul!(v2u2_hat, plan, v2u2)

    dvorhattdt .= g_hat .+ ν .* Delta_hat .* vort_hat .- DxsqDysq_hat .* uv_hat .- Dxy_hat .* v2u2_hat

    ## dealiasing
    dvorhattdt[div(Nsub,2) + 1:end, :] .= 0.0im
    dvorhattdt[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im

    return dvorhattdt
end

# ## The spatial domain and its discretization

L = 2π
κ₀ = 2π/L
N = 128
Nsub = 84
x = y = (L/N):(L/N):L

# ## Test convergence to one-mode steady state

ν = 1.0e-0 # viscosity

Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat = get_operators(N, κ₀)

vort_init = field_from_modes(L, N, 4)
vort_init_hat = rfft(vort_init)

g_steady = field_from_modes(L, N, 1)
g_steady_hat = rfft(g_steady)

vort_steady_hat = - g_steady_hat ./ Delta_hat
vort_steady_hat[1, 1] = 0.0im
vort_steady = irfft(vort_steady_hat, N)

vort_hat = copy(vort_init_hat)
vort = irfft(vort_hat, N)

u_hat = similar(vort_hat)
v_hat = similar(vort_hat)
u = similar(vort)
v = similar(vort)
uv = similar(vort)
v2u2 = similar(vort)
uv_hat = similar(vort_hat)
v2u2_hat = similar(vort_hat)
plan = plan_rfft(vort, flags=FFTW.MEASURE)
plan_inv = plan_irfft(vort_hat, N, flags=FFTW.MEASURE)

operators = Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat
vars = ν, g_steady_hat, N, Nsub
auxs = u_hat, v_hat, u, v, uv, v2u2, uv_hat, v2u2_hat
plans = plan, plan_inv

params = (
    operators,
    vars,
    auxs,
    plans
)

nothing
#

tspan = (0.0, 10.0)

prob = ODEProblem(nsepervort_hat_rhs!, vort_init_hat, tspan, params)
nothing

#

OrdinaryDiffEq.ForwardDiff.can_dual(::Type{ComplexF64}) = true

sol = solve(prob, Vern9()) # QNDF() etc.

sol.retcode

# Distância entre a solução estacionária e a solução ao final do intervalo de tempo.

norm(sol.u[end] - vort_steady_hat) * L/N /N

#

vort = irfft(sol.u[end], N)
nothing

# No espaço físico:

norm(vort .- vort_steady) * L / N

# Na norma do máximo:

maximum(abs, vort .- vort_steady)

#

heatmap(x, y, vort, xlabel="x", ylabel="y", title="vorticity", titlefont=12)

#

surface(x, y, vort, xlabel="x", ylabel="y", title="vorticity", titlefont=12)

#

surface(x, y, vort .- vort_steady, xlabel="x", ylabel="y", title="vorticity", titlefont=12)

# We can check how the $L^2$ distance of the vorticity to the steady state evolves.

error = [norm(sol.u[n] .- vort_steady_hat) * L/N / N for n in eachindex(sol.u)]
nothing

#

plt1 = plot(sol.t, error, title="Enstrophy convergence", titlefont = 10, xaxis="time", yaxis = "enstrophy", label=false)

plt2 = plot(sol.t[1:div(end, 100)], error[1:div(end, 100)], title="Enstrophy convergence", titlefont = 10, xaxis="time", yaxis = "enstrophy", label=false)

plt3 = plot(sol.t[div(end,100):div(end, 50)], error[div(end,100):div(end, 50)], title="Enstrophy convergence", titlefont = 10, xaxis="time", yaxis = "enstrophy", label=false)

plt4 = plot(sol.t[end-1000:end], error[end-1000:end], title="Enstrophy convergence", titlefont = 10, xaxis="time", yaxis = "enstrophy", label=false)

plot(plt1, plt2, plt3, plt4, layout = 4)
#


text:
function randint(start, end_)
    r = rand() * (end_ - start)
    r = floor(r) + start
    r = Int(r)
end

function swap(lst, i, j)
    t = lst[i]
    lst[i] = lst[j]
    lst[j] = t
end

function partition(lst, p, r, random=false)
    if random == true
        r_ = randint(p, r+1)
        swap(lst, r, r_)
    end
    x = lst[r]
    i = p
    for j in p:r-1
        if lst[j] <= x
            swap(lst, j, i)
            i += 1
        end
    end
    swap(lst, r, i)
    return i
end

function quicksort(lst, p, r, random=false)
    if p < r
        q = partition(lst, p, r, random)
        quicksort(lst, p, q-1)
        quicksort(lst, q+1, r)
    end
end

lst = [8, 11, 2, 1, 4, 19, 9, 14, 13, 5, 6, 22, ]
len = length(lst)
println(lst)
quicksort(lst, 1, len, true)
println(lst)

max_stars_repo_name: Bingwen-Hu/hackaway

max_issues_repo_name: Bingwen-Hu/hackaway

max_forks_repo_name: Bingwen-Hu/hackaway


text:
print("Last updated: $(read(`date`, String))")

# Load environment
import Pkg; Pkg.activate("../../../")

# Import Libraries
using Turing
using Turing: Variational
using Distributions
using Distances
using PyPlot
using StatsFuns
import Random
using Flux
import LinearAlgebra
using JSON3

# For getting quantiles along array dimensions
function quantiles(X, q; dims, drop=false)
  Q = mapslices(x -> quantile(x, q), X, dims=dims)
  out = drop ? dropdims(Q, dims=dims) : Q
  return out
end

# To extract parameters from ADVI model.
function make_extractor(m, q, nsamples=1000)
    qsamples = rand(q, nsamples)
    _, sym2range = Variational.bijector(m; sym_to_ranges = Val(true));
    return sym -> qsamples[collect(sym2range[sym][1]), :]
end
;

# Squared-exponential covariance function
function sqexp_cov_fn(D, alpha, rho; d_is_squared=false)
    if d_is_squared
        return alpha ^ 2 * exp.(-0.5 * D/(rho^2))
    else
        return alpha ^ 2 * exp.(-0.5 * (D/rho) .^ 2)
    end
end

@model function GP(y, X, m_alpha=0.0, s_alpha=1.0, m_rho=0.0, s_rho=1.0, m_sigma=0.0, s_sigma=1.0)
    # Distance matrix.
    D2 = pairwise(Distances.SqEuclidean(), X, dims=1)

    # Priors.
    alpha ~ LogNormal(m_alpha, s_alpha)
    rho ~ LogNormal(m_rho, s_rho)
    sigma ~ LogNormal(m_sigma, s_sigma)

    # Realized covariance function
    K = sqexp_cov_fn(D2, alpha, rho, d_is_squared=true)

    # Sampling Distribution.
    y ~ MvNormal(K + LinearAlgebra.I * sigma^2)  # mean=0, covariance=K.
end;

# Read data.

# Directory where all simulation data are stored.
data_path = joinpath(@__DIR__, "../data/gp-data-N30.json")

# Load data in JSON format.
data = let
    x = open(f -> read(f, String), data_path)
    JSON3.read(x)
end

# Reshape data if needed.
y = Float64.(data[:y])
X = Float64.(reshape(data[:x], length(y), 1))

f = Float64.(data[:f])
x_grid = Float64.(data[:x_grid])

N = size(X, 1);

# Plot data
plt.scatter(vec(X), y, label="Data")
plt.plot(x_grid, f, c="grey", ls=":", label="True f(x)")
plt.xlabel("x")
plt.ylabel("y")
plt.legend();

# Create model.
m = begin
    m_alpha = 0.0
    s_alpha = 0.1
    m_rho = 0.0
    s_rho = 1.0
    m_sigma = 0.0
    s_sigma = 1.0
    GP(y, X, m_alpha, s_alpha, m_rho, s_rho, m_sigma, s_sigma)
end;

# Fit via ADVI. You can also use HMC.
Random.seed!(0)

# initialize variational distribution (optional)
q0 = Variational.meanfield(m)

# NOTE: ADVI(num_elbo_samples, max_iters)

# Compile
@time q = vi(m, ADVI(1, 1), q0, optimizer=Flux.ADAM(1e-1));

# RUN
@time q = vi(m, ADVI(1, 2000), q0, optimizer=Flux.ADAM(1e-1));

# Get posterior samples
extract_gp = make_extractor(m, q)
alpha = vec(extract_gp(:alpha));
rho = vec(extract_gp(:rho));
sigma = vec(extract_gp(:sigma));
advi_samples = Dict(:alpha => alpha, :rho => rho, :sigma => sigma);

# Fit via HMC.
Random.seed!(0)
burn = 1000
nsamples = 1000

# Compile
@time _ = sample(m, HMC(0.01, 1), 1)

# Run
@time chain = sample(m, HMC(0.01, 100), burn + nsamples)

# Get posterior samples
alpha = vec(group(chain, :alpha).value.data[end-nsamples:end, :, 1]);
rho = vec(group(chain, :rho).value.data[end-nsamples:end, :, 1]);
sigma = vec(group(chain, :sigma).value.data[end-nsamples:end, :, 1]);
hmc_samples = Dict(:alpha => alpha, :rho => rho, :sigma => sigma);

# Fit via NUTS.
Random.seed!(7)

# Compile
@time _ = sample(m, NUTS(4, 0.8), 9);

# Run
@time chain = begin
    nsamples = 1000  # number of MCMC samples
    nadapt = 1000  # number of iterations to adapt tuning parameters in NUTS
    iterations = nsamples + nadapt
    target_accept_ratio = 0.8

    sample(m, NUTS(nadapt, target_accept_ratio, max_depth=10), iterations);
end

# Get posterior samples
alpha = vec(group(chain, :alpha).value.data[:, :, 1]);
rho = vec(group(chain, :rho).value.data[:, :, 1]);
sigma = vec(group(chain, :sigma).value.data[:, :, 1]);
nuts_samples = Dict(:alpha => alpha, :rho => rho, :sigma => sigma);

# This funciton returns a function for predicting at new points given parameter values.
function make_gp_predict_fn(Xnew, y, X)
    N = size(X, 1)
    M = size(Xnew, 1)
    Q = N + M
    Z = [Xnew; X]
    D2 = pairwise(SqEuclidean(), Z, dims=1)

    return (alpha, rho, sigma) -> let
        K = sqexp_cov_fn(D2, alpha, rho, d_is_squared=true) + LinearAlgebra.I * sigma^2
        Koo_inv = inv(K[(M+1):end, (M+1):end])
        Knn = K[1:M, 1:M]
        Kno = K[1:M, (M+1):end]
        C = Kno * Koo_inv
        m = C * y
        S = Matrix(LinearAlgebra.Hermitian(Knn - C * Kno'))
        mvn = MvNormal(m, S)
        rand(mvn)
    end
end

# Function for plotting parameter posterior.
function plot_post(samples, name; bins=nothing, suffix="")
    param = samples[name]
    plt.hist(param, density=true, bins=bins)
    plt.xlabel(name)
    plt.ylabel("density")
    if suffix !== ""
        suffix = "($(suffix))"
    end
    plt.title("Histogram of $(name) $(suffix)");
end

function plot_fn_posterior(samples; figsize=(12, 4), figsize_f=figsize, suffix="")
    # Get parameters
    alpha = samples[:alpha]
    rho = samples[:rho]
    sigma = samples[:sigma]
    nsamps = length(alpha)

    # Plot parameters posterior.
    plt.figure(figsize=figsize)
    plt.subplot(1, 3, 1)
    plot_post(samples, :alpha, bins=30, suffix=suffix)
    plt.subplot(1, 3, 2)
    plot_post(samples, :rho, bins=30, suffix=suffix)
    plt.subplot(1, 3, 3)
    plot_post(samples, :sigma, bins=30, suffix=suffix)
    plt.axvline(data[:sigma], ls=":", label="truth", color="red")
    plt.legend()
    plt.tight_layout()

    # Make predictions at new locations.
    X_new = reshape(collect(range(-3.5, 3.5, length=100)), 100, 1)
    gp_predict = make_gp_predict_fn(X_new, y, X)

    ynew = [gp_predict(alpha[m], rho[m], sigma[m]) for m in 1:nsamps]
    ynew = hcat(ynew...)

    # Summary statistics of posterior function.
    pred_mean = mean(ynew, dims=2)
    pred_lower = quantiles(ynew, 0.025, dims=2)
    pred_upper = quantiles(ynew, 0.975, dims=2)

    if suffix !== ""
        suffix = "($(suffix))"
    end

    # Plot predictions.
    plt.figure(figsize=figsize_f)
    plt.plot(vec(X_new), vec(pred_mean), color="blue", label="Posterior predictive mean")
    plt.plot(x_grid, f, color="red", ls=":", label="True f(x)")

    plt.scatter(vec(X), vec(y), color="black", label="Data")
    plt.fill_between(vec(X_new), vec(pred_upper), vec(pred_lower), color="blue",
                     alpha=0.2, label="95% Credible Interval")
    plt.legend(loc="upper left", fontsize=8)
    plt.title("GP Posterior predictive with 95% credible interval $(suffix)");

    plt.xlabel("x")
    plt.ylabel("y")

    plt.ylim(-2, 2)
    plt.tight_layout()
end

plot_fn_posterior(advi_samples, figsize=(9, 3), suffix="ADVI")

plot_fn_posterior(hmc_samples, figsize=(9, 3), suffix="HMC");

plot_fn_posterior(nuts_samples, figsize=(9, 3), suffix="NUTS");

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
include("coal.jl")
include("coal_sort.jl")


text:
""" Type for location at the cell centres. """
struct Centre <: AbstractLocation end

# US spelling alias
Center = Centre

""" Type for location at the cell faces. """
struct Face <: AbstractLocation end

"""
    struct Field1D{LX<:AbstractLocation, G, D} <: AbstractField

A field datatype for 1D objects containing data of type `D`.

$(TYPEDFIELDS)
"""
struct Field1D{LX<:AbstractLocation, G, D} <: AbstractField
    "Array with the values of the field."
    data :: D
    "The grid on which the field lives."
    grid :: G

    Field1D(LX, data::D, grid::G) where {G, D} = new{LX, G, D}(data, grid)
end

function Field1D(LX, data::Array, grid::Grid1D)
    nx, hx = grid.nx, grid.hx

    data_with_halos = OffsetArray(zeros(nx + 2hx), -hx)

    @. data_with_halos[1:nx] = data

    field = Field1D(LX, data_with_halos, grid)

    fill_halos!(field)

    return field
end

"""
    struct Field2D{LX<:AbstractLocation, LY<:AbstractLocation, G, D} <: AbstractField

A field datatype for 2D objects containing data of type `D`.

$(TYPEDFIELDS)
"""
struct Field2D{LX<:AbstractLocation, LY<:AbstractLocation, G, D} <: AbstractField
    "Array with the values of the field."
    data :: D
    "The grid on which the field lives."
    grid :: G

    Field2D(LX, LY, data::D, grid::G) where {G, D} = new{LX, LY, G, D}(data, grid)
end

function Field2D(LX, LY, data::Array, grid::Grid2D)
    nx, hx = grid.nx, grid.hx
    ny, hy = grid.ny, grid.hy

    data_with_halos = OffsetArray(zeros(nx + 2hx, ny + 2hy), -hx, -hy)

    @. data_with_halos[1:nx, 1:ny] = data

    field = Field2D(LX, LY, data_with_halos, grid)

    fill_halos!(field)

    return field
end

"""
    Field(LX, data, grid::Grid1D)

Constructs a 1D field of `data` at location `LX` on `grid`.
"""
Field(LX, data, grid::Grid1D) = Field1D(LX, data, grid)

"""
    Field(LX, LY, data, grid::Grid2D)

Constructs a 2D field of `data` at location `(LX, LY)` on `grid`.
"""
Field(LX, LY, data, grid::Grid2D) = Field2D(LX, LY, data, grid)

#####
##### Intepolations
#####

𝐼xᶜ(i, f::Field1D{Face})   = (f.data[i] + f.data[i+1]) / 2
𝐼xᶠ(i, f::Field1D{Centre}) = (f.data[i-1] + f.data[i]) / 2

𝐼xᶜᶜ(i, j, f::Field2D{Face, Centre})   = (f.data[i, j] + f.data[i+1, j]) / 2
𝐼xᶠᶜ(i, j, f::Field2D{Centre, Centre}) = (f.data[i-1, j] + f.data[i, j]) / 2
𝐼yᶜᶜ(i, j, f::Field2D{Centre, Face})   = (f.data[i, j] + f.data[i, j+1]) / 2
𝐼yᶜᶠ(i, j, f::Field2D{Centre, Centre}) = (f.data[i, j-1] + f.data[i, j]) / 2

"""
    𝐼x!(output::Field1D, input::Field1D{<:Any, Grid1D{Periodic}})

Interpolates a 1D `input` field to the location where the `output` field lives for 1D grids
with periodic boundary conditions.
"""
function 𝐼x!(output::Field1D{Centre}, input::Field1D{Face, Grid1D{Periodic}})
    nx = input.grid.nx

    for i in 1:nx
        output.data[i] = 𝐼xᶜ(i, input)
    end

    fill_halos!(output)

    return nothing
end

function 𝐼x!(output::Field1D{Face}, input::Field1D{Centre, Grid1D{Periodic}})
    nx = input.grid.nx

    for i in 1:nx
        output.data[i] = 𝐼xᶠ(i, input)
    end

    fill_halos!(output)

    return nothing
end

"""
    𝐼x!(output::Field2D, input::Field2D{<:Any, <:Any, Grid2D{Periodic, Periodic}})

Interpolates a 2D `input` field to the location where the `output` field lives for 2D grids
with periodic boundary conditions.
"""
function 𝐼x!(output::Field2D{Face, Centre}, input::Field2D{Centre, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny

    for j in 1:ny, i in 1:nx
        output.data[i, j] = 𝐼xᶠᶜ(i, j, input)
    end

    fill_halos!(output)

    return nothing
end

function 𝐼x!(output::Field2D{Centre, Centre}, input::Field2D{Face, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny

    for j in 1:ny, i in 1:nx
        output.data[i, j] = 𝐼xᶜᶜ(i, j, input)
    end

    fill_halos!(output)

    return nothing
end

"""
    𝐼y!(output::Field2D, input::Field2D{<:Any, <:Any, Grid1D{Periodic, Periodic}})

Interpolates a 2D `input` field to the location where the `output` field lives for 2D grids
with periodic boundary conditions.
"""
function 𝐼y!(output::Field2D{Centre, Face}, input::Field2D{Centre, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny

    for j in 1:ny, i in 1:nx
        output.data[i, j] = 𝐼yᶜᶠ(i, j, input)
    end

    fill_halos!(output)

    return nothing
end

function 𝐼y!(output::Field2D{Centre, Centre}, input::Field2D{Centre, Face, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny

    for j in 1:ny, i in 1:nx
        output.data[i, j] = 𝐼yᶜᶜ(i, j, input)
    end

    fill_halos!(output)

    return nothing
end

#####
##### Derivatives
#####

δxᶜ(i, f::Field1D{Face})   = f.data[i+1] - f.data[i]
δxᶜ(i, f::Field1D{Centre}) = f.data[i+1] - f.data[i-1]
δxᶠ(i, f::Field1D{Centre}) = f.data[i]   - f.data[i-1]
δxᶠ(i, f::Field1D{Face})   = f.data[i+1] - f.data[i-1]

δxᶜᶜ(i, j, f::Field2D{Face, Centre})   = f.data[i+1, j] - f.data[i, j]
δxᶜᶜ(i, j, f::Field2D{Centre, Centre}) = f.data[i+1, j] - f.data[i-1, j]
δxᶠᶜ(i, j, f::Field2D{Centre, Centre}) = f.data[i, j]   - f.data[i-1, j]
δxᶠᶜ(i, j, f::Field2D{Face, Centre})   = f.data[i+1, j] - f.data[i-1, j]
δyᶜᶜ(i, j, f::Field2D{Centre, Face})   = f.data[i, j+1] - f.data[i, j]
δyᶜᶜ(i, j, f::Field2D{Centre, Centre}) = f.data[i, j+1] - f.data[i, j-1]
δyᶜᶠ(i, j, f::Field2D{Centre, Centre}) = f.data[i, j]   - f.data[i, j-1]
δyᶜᶠ(i, j, f::Field2D{Centre, Face})   = f.data[i, j+1] - f.data[i, j-1]

"""
    ∂x!(output::Field1D, input::Field1D{<:Any, Grid1D{Periodic}})

Compute the derivative of a 1D `input` field onto the location where the `output` field lives
for 1D grids with periodic boundary conditions.
"""
function ∂x!(output::Field1D{Centre}, input::Field1D{Face, Grid1D{Periodic}})
    nx, dx = input.grid.nx, input.grid.dx

    for i in 1:nx
        output.data[i] = δxᶜ(i, input) / dx
    end

    fill_halos!(output)

    return nothing
end

function ∂x!(output::Field1D{Face}, input::Field1D{Centre, Grid1D{Periodic}})
    nx, dx = input.grid.nx, input.grid.dx

    for i in 1:nx
        output.data[i] = δxᶠ(i, input) / dx
    end

    fill_halos!(output)

    return nothing
end

function ∂x!(output::Field1D{Centre}, input::Field1D{Centre, Grid1D{Periodic}})
    nx, dx = input.grid.nx, input.grid.dx

    for i in 1:nx
        output.data[i] = δxᶜ(i, input) / (2*dx)
    end

    fill_halos!(output)

    return nothing
end

function ∂x!(output::Field1D{Face}, input::Field1D{Face, Grid1D{Periodic}})
    nx, dx = input.grid.nx, input.grid.dx

    for i in 1:nx
        output.data[i] = δxᶠ(i, input) / (2*dx)
    end

    fill_halos!(output)

    return nothing
end

"""
    ∂x!(output::Field2D, input::Field2D{<:Any, <:Any, Grid2D{Periodic, Periodic}})

Compute the ``x`` derivative of a 2D `input` field onto the location where the `output` field lives
for 2D grids with periodic boundary conditions.
"""
function ∂x!(output::Field2D{Centre, Centre}, input::Field2D{Face, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dx = input.grid.dx

    for j in 1:ny, i = 1:nx
        output.data[i, j] = δxᶜᶜ(i, j, input) / dx
    end

    fill_halos!(output)

    return nothing
end

function ∂x!(output::Field2D{Face, Centre}, input::Field2D{Face, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dx = input.grid.dx

    for j in 1:ny, i = 1:nx
        output.data[i, j] = δxᶠᶜ(i, j, input) / (2*dx)
    end

    fill_halos!(output)

    return nothing
end

function ∂x!(output::Field2D{Face, Centre}, input::Field2D{Centre, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dx = input.grid.dx

    for j in 1:ny, i in 1:nx
        output.data[i, j] = δxᶠᶜ(i, j, input)/dx
    end

    fill_halos!(output)

    return nothing
end

function ∂x!(output::Field2D{Centre, Centre}, input::Field2D{Centre, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dx = input.grid.dx

    for j in 1:ny, i in 1:nx
        output.data[i, j] = δxᶜᶜ(i, j, input) / (2*dx)
    end

    fill_halos!(output)

    return nothing
end

"""
    ∂y!(output::Field2D, input::Field2D{<:Any, <:Any, Grid2D{Periodic, Periodic}})

Compute the ``y`` derivative of a 2D `input` field onto the location where the `output` field lives
for 2D grids with periodic boundary conditions.
"""
function ∂y!(output::Field2D{Centre, Centre}, input::Field2D{Centre, Face, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dy = input.grid.dy

    for j in 1:ny, i in 1:nx
        output.data[i, j] = δyᶜᶜ(i, j, input) / dy
    end

    fill_halos!(output)

    return nothing
end

function ∂y!(output::Field2D{Centre, Face}, input::Field2D{Centre, Face, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dy = input.grid.dy

    for j in 1:ny, i in 1:nx
        output.data[i, j] = δyᶜᶠ(i, j, input) / (2*dy)
    end

    fill_halos!(output)

    return nothing
end

function ∂y!(output::Field2D{Centre, Face}, input::Field2D{Centre, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dy = input.grid.dy

    for j in 1:ny, i in 1:nx
        output.data[i, j] = δyᶜᶠ(i, j, input) / dy
    end

    fill_halos!(output)

    return nothing
end

function ∂y!(output::Field2D{Centre, Centre}, input::Field2D{Centre, Centre, Grid2D{Periodic, Periodic}})
    nx, ny = input.grid.nx, input.grid.ny
    dy = input.grid.dy

    for j in 1:ny, i in 1:nx
        output.data[i, j] = δyᶜᶜ(i, j, input) / (2*dy)
    end

    fill_halos!(output)

    return nothing
end

#####
##### Filling halos
#####

"""
    fill_halos!(field::Field1D{<:Any, Grid1D{Periodic}})

Fill halos for a 1D `field` that lives on a grid with periodic boundary conditions.
"""
function fill_halos!(field::Field1D{<:Any, Grid1D{Periodic}})
    nx, hx = field.grid.nx, field.grid.hx

    for i in 1:hx
        field.data[nx+i] = field.data[i]
        field.data[-i+1] = field.data[nx-i+1]
    end

    return nothing
 end

"""
    fill_halos!(field::Field2D{<:Any, <:Any, Grid2D{Periodic, Periodic}})

Fill halos for a 2D `field` that lives on a grid with periodic boundary conditions in both
directions.
"""
function fill_halos!(field::Field2D{<:Any, <:Any, Grid2D{Periodic, Periodic}})
    nx, hx = field.grid.nx, field.grid.hx
    ny, hy = field.grid.ny, field.grid.hy

    for j in 1:ny, i in 1:hx
        field.data[nx+i, j] = field.data[i, j]
        field.data[-i+1, j] = field.data[nx-i+1, j]
    end

    for j in 1:hy, i in 1:nx
        field.data[i, ny+j] = field.data[i, j]
        field.data[i, -j+1] = field.data[i, ny-j+1]
    end

    return nothing
 end


text:
# Sum of the positive components

export SumPositive

"""
**Sum of the positive coefficients**

    SumPositive()

Returns the function
```math
f(x) = ∑_i \\max\\{0, x_i\\}.
```
"""
struct SumPositive <: ProximableFunction end

is_separable(f::SumPositive) = true
is_convex(f::SumPositive) = true

function (f::SumPositive)(x::AbstractArray{T}) where T <: Real
    return sum(xi -> max(xi, 0), x)
end

function prox!(y::AbstractArray{R}, f::SumPositive, x::AbstractArray{R}, gamma::R=R(1)) where R <: Real
    fsum = R(0)
    for i in eachindex(x)
        y[i] = x[i] < gamma ? (x[i] > 0 ? R(0) : x[i]) : x[i]-gamma
        fsum += y[i] > 0 ? y[i] : R(0)
    end
    return fsum
end

function gradient!(y::AbstractArray{R}, f::SumPositive, x::AbstractArray{R}) where R <: Real
    y .= max.(0, sign.(x))
    return sum(xi -> max(xi, 0), x)
end

fun_name(f::SumPositive) = "Sum of the positive coefficients"
fun_dom(f::SumPositive) = "AbstractArray{Real}"
fun_expr(f::SumPositive) = "x ↦ sum(max(0, x))"

function prox_naive(f::SumPositive, x::AbstractArray{R}, gamma::R=R(1)) where R <: Real
    y = copy(x)
    indpos = x .> 0
    y[indpos] = max.(R(0), x[indpos] .- gamma)
    return y, sum(max.(R(0), y))
end

# ######################### #
# Prox with multiple gammas #
# ######################### #

function prox!(y::AbstractArray{R}, f::SumPositive, x::AbstractArray{R}, gamma::AbstractArray{R}) where R <: Real
    fsum = R(0)
    for i in eachindex(x)
        y[i] = x[i] < gamma[i] ? (x[i] > 0 ? R(0) : x[i]) : x[i]-gamma[i]
        fsum += y[i] > 0 ? y[i] : R(0)
    end
    return fsum
end

function prox_naive(f::SumPositive, x::AbstractArray{R}, gamma::AbstractArray{R}) where R <: Real
    y = copy(x)
    indpos = x .> 0
    y[indpos] = max.(R(0), x[indpos] .- gamma[indpos])
    return y, sum(max.(R(0), y))
end


text:
export iauS2c
"""
**  Convert spherical coordinates to Cartesian.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  vector/matrix support function.
**
**  Given:
**     theta    double       longitude angle (radians)
**     phi      double       latitude angle (radians)
**
**  Returned:
**     c        double[3]    direction cosines
**
**  This revision:  2013 June 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauS2c(double theta, double phi, double c[3])
function iauS2c(theta::Real, phi::Real)
    c = zeros(Float64, 3)

    ccall((:iauS2c, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Ptr{Cdouble}),
       convert(Float64, theta),
       convert(Float64, phi),
       c)

   return SVector{3}(c)
end


text:
s = "rosetta code phrase reversal"

println("The original phrase.")
println("   ", s)

println("Reverse the string.")
t = reverse(s)
println("   ", t)

println("Reverse each individual word in the string.")
t = join(map(reverse, split(s, " ")), " ")
println("   ", t)

println("Reverse the order of each word of the phrase.")
t = join(reverse(split(s, " ")), " ")
println("   ", t)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here
move_to(cr, 128.0, 25.6);
line_to(cr, 230.4, 230.4);
rel_line_to(cr, -102.4, 0.0);
curve_to(cr, 51.2, 230.4, 51.2, 128.0, 128.0, 128.0);
close_path(cr);

move_to(cr, 64.0, 25.6);
rel_line_to(cr, 51.2, 51.2);
rel_line_to(cr, -51.2, 51.2);
rel_line_to(cr, -51.2, -51.2);
close_path(cr);

set_line_width(cr, 10.0);
set_source_rgb(cr, 0, 0, 1);
fill_preserve(cr);
set_source_rgb(cr, 0, 0, 0);
stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_fill_and_stroke2.png");


text:
# SucVector
# ---------
#
# Almost the same as CompactBitVector, but the layout of data is different.
# Four chunks of bits are stored in a block and blocks are stored in a vector.
# As a result, all data are interleaved and thus the CPU cache may be used
# efficiently. Note that in v0.3 the blocks will not be packed in a vector.

immutable Block
    # bit chunks (64bits × 4 = 256bits)
    chunks::NTuple{4,Uint64}
    # large block
    large::Uint32
    # 8-bit extension of the large block (32 + 8 = 40 bits in total)
    ext::Uint8
    # small blocks
    smalls::NTuple{3,Uint8}
end

const bits_per_chunk =  64
const bits_per_block = 256

function Block(chunks::NTuple{4,Uint64}, offset::Int)
    @assert offset ≤ 2^40 - 1
    a =     convert(Uint8, count_ones(chunks[1]))
    b = a + convert(Uint8, count_ones(chunks[2]))
    c = b + convert(Uint8, count_ones(chunks[3]))
    Block(chunks, offset, offset >>> 32, (a, b, c))
end

type SucVector <: AbstractIndexableBitVector
    blocks::Vector{Block}
    len::Int
end

SucVector() = SucVector(Block[], 0)

function SucVector(src::Union(BitVector,Vector{Bool}))
    len = length(src)
    n_blocks = div(len, bits_per_block) + 1
    blocks = Array(Block, n_blocks)
    offset = 0
    for i in 1:n_blocks
        chunks = read_4chunks(src, (i - 1) * bits_per_block + 1)
        blocks[i] = Block(chunks, offset)
        for j in 1:4
            offset += count_ones(chunks[j])
        end
    end
    return SucVector(blocks, len)
end

function read_chunk(src, from::Int)
    # read a 64-bit chunk from a bitvector
    chunk = zero(Uint64)
    for i in from:from+63
        chunk <<= 1
        if i ≤ endof(src)
            chunk += src[i] ≠ 0
        end
    end
    return chunk
end

function read_4chunks(src, from::Int)
    # read four 64-bit chunks from a bitvector at once
    a = read_chunk(src, from)
    b = read_chunk(src, from + 64 * 1)
    c = read_chunk(src, from + 64 * 2)
    d = read_chunk(src, from + 64 * 3)
    (a, b, c, d)
end

function convert(::Type{SucVector}, v::Union(BitVector,Vector{Bool}))
    return SucVector(v)
end

Base.length(v::SucVector) = v.len

function getindex(v::SucVector, i::Integer)
    if !(1 ≤ i ≤ endof(v))
        throw(BoundsError())
    end
    return unsafe_getindex(v, i)
end

function unsafe_getindex(v::SucVector, i::Integer)
    q, r = divrem(i - 1, bits_per_block)
    block = v.blocks[q+1]
    q, r = divrem(r, bits_per_chunk)
    chunk = block.chunks[q+1]
    return bitat(chunk, r + 1)
end

function rank1(v::SucVector, i::Integer)
    if !(0 ≤ i ≤ endof(v))
        throw(BoundsError())
    end
    return unsafe_rank1(v, i)
end

function unsafe_rank1(v::SucVector, i::Integer)
    q, r = divrem(i - 1, bits_per_block)
    @inbounds block = v.blocks[q+1]
    ret = 0
    # large block
    ret += block.ext << 32 + convert(Int64, block.large)
    # small block
    q, r = divrem(r, bits_per_chunk)
    if q > 0
        @inbounds ret += block.smalls[q]
    end
    # remaining bits
    @inbounds chunk = block.chunks[q+1]
    ret += count_ones(chunk & lmask(Uint64, r + 1))
    return ret
end


text:
module Jus

using Match
using JSON3
using JSON3: StructTypes
using Base.Filesystem
using HTTP
using Pkg
using DefaultApplication
using Sockets
using Sockets: TCPServer
using Mustache

import Base.@kwdef

include("types.jl")
include("vars.jl")
include("server.jl")
include("gen.jl")
include("local.jl")

const CONVERT_ENUM = r"^enum:(.*)$"

verbose = false

strdefault(str, default) = str == "" ? default : str

function usage()
    print("""Usage: jus NAMESPACE [-v] [-s ADDR | -c ADDR [-x SECRET]] [CMD [ARGS...]]

NAMESPACE is this jus instance's namespace

-i          interactive -- don't start Jus, just run the Julia repl
-s ADDR     serve HTTP requests on ADDR if ADDR starts with '/', use a UNIX domain socket
-c ADDR     connect to ADDR using NAMESPACE
-x SECRET   use secret to prove ownership of NAMESPACE. If NAMESPACE does not yet exist,
            the server creates it and associates SECRET with it.
-d DIR      add dir to the http server file path
-e EXPR     evaluate Julia expression
-b TYPE     open browser on TYPE
-v          verbose

COMMANDS
set [[-c] [PATH VALUE] | [-m PATH]]...
  Set a variable, optionally creating it if it does not exist. NAMESPACE defaults to the
  current namespace.

  -c    create a variable
  -m    set only metadata for a variable (PATH should contain META -- see below)

  The format of PATH (without spaces) is
    NAME ['.' NAME]... [: META]
    or
    NAMESPACE '/' ID ['.' NAME]... [: META]

  NAMESPACE '/' ID defaults to ROOT/0.
  The name '?' means to create a new numbered child of NAMESPACE '/' ID.
  The name '?N' refers to the Nth variable created in this command (to help create trees)

  The format of META is
    NAME [',' META]...
    or
    NAME '=' VALUE [',' META]...

  Set returns a JSON list of ids the variables that were created:
    {"result": [id1, id2, ...]}

  COMMON META VALUE NAMES
  app       create an instance of the named application's object model
  path      the path from the parent's value to the corresponding value in the object model
            an empty path refers to the parent
  monitor   whether to monitor the value in the object model value is yes/no/true/false/on/off
  call      call a function on the value, optionally with a context object

get PATH...
  See set command for PATH format.

  returns values of requested variables:
    {"result": [value1, value2, ...]}

delete ID
  Remove a variable and all of its children

observe ID...
  receive updates whenever variables or their children change. Update format is
    {"update": [id1, name1, id2, name2, ...]}

  returns the ids and values for the given paths plus all of their descendants
    {"result": [id1, value1, id2, value2, ...]}
""")
    exit(1)
end

log(args...) = @info join(args)

function shutdown(config)
    println("SHUTTING DOWN")
end

function parseAddr(config::Config, s)
    m = match(r"^(([^:]+):)?([^:]+)$", s)
    config.host == m[2] ? "127.0.0.1" : m[2]
    config.port = parse(UInt16, m[3])
end

function abort(msg...)
    println(Base.stderr, msg...)
    exit(1)
end

output(ws; args...) = output(ws, args)
function output(ws, data)
    write(ws, JSON3.write(data))
    flush(ws)
    @debug("WROTE: $(JSON3.write(data))")
end

input(ws) = JSON3.read(readavailable(ws))

# no need for JSON with FakeWS
output(ws::FakeWS, data) = put!(ws.output, data)
input(ws::FakeWS) = take!(ws.input)

function client(config::Config)
    if config.namespace.secret === "" abort("Secret required") end
    @debug("CLIENT $(config.namespace) connecting to ws//$(config.host):$(config.port)")
    HTTP.WebSockets.open("ws://$(config.host):$(config.port)") do ws
        output(ws, (namespace = config.namespace, secret = config.namespace.secret))
        output(ws, config.args)
        result = input(ws) # read one message
        @debug("RESULT: $(result)")
        if config.args[1] == "observe"
            @debug("READING UPDATES")
            while true
                result = input(ws)
                @debug("RESULT: $(result)")
            end
        end
    end
end

function exec(serverfunc, args::Vector{String}; config = Config())
    local browse = ""

    interactive = "-i" in args
    interactive && (args = filter(x-> x != "-i", args))
    if interactive && (isempty(args) || match(r"^-.*$", args[1]) !== nothing)
        i = 1
    else
        # name required -- only one instance per name allowed
        (length(args) === 0 || match(r"^-.*$", args[1]) !== nothing) && !("-i" in args) && usage()
        config.serverfunc = serverfunc
        config.namespace.name = args[1]
        i = 2
    end
    while i <= length(args)
        @match args[i] begin
            "-i" => (interactive = true)
            "-e" => Main.eval(Meta.parse(args[i += 1]))
            "-d" => add_file_dir(args[i += 1])
            "-v" => (config.verbose = true)
            "-x" => (config.namespace.secret = args[i += 1])
            "-c" => parseAddr(config, args[i += 1])
            "-b" => (browse = args[i += 1])
            "-s" => begin
                parseAddr(config, args[i += 1])
                config.server = true
            end
            unknown => begin
                @debug("MATCHED DEFAULT: $(args[i:end])")
                push!(config.args, args[i:end]...)
                i = length(args)
                break
            end
        end
        i += 1
    end
    interactive && return
    atexit(()-> shutdown(config))
    config.server && browse != "" && present(config, browse)
    (config.server ? server : client)(config)
    config
end

"""
    start(data | type;
        port = 0, host = "0.0.0.0", dirs = [], async = true, metadata = (;),
        create_dir = "", output_dir = create_dir)
    start(config, server, data; dirs = [], async = true, metadata = (;), create_dir = "")

Starts Jus in a new task on host and port (port 0 means pick a random port)
and, if given a type, opens a browser page editing the type.

output_dir specifies a directory (if any) in which to place generated viewdefs
create_dir specifies a directory to create and also supplies the default value for output_dir

Returns the config.
"""
function start(data = nothing;
               port = 0, host = nothing, dirs = [], async = true, metadata = (;),
               create_output = "", output_dir = create_output)
    port, socket = host === nothing ? Sockets.listenany(port) : Sockets.listenany(IPv4(host), port)
    confhost = getsockname(socket)[1]
    hostname = host !== nothing ? host : confhost == ip"0.0.0.0" ? "localhost" : string(confhost)
    println("HOST: $(confhost)")
    start(Config(; output_dir, port, host = confhost, hostname), socket, data; dirs, async, metadata, create_output)
end

function start(config::Config, socket::TCPServer, data; dirs = [], async = true, metadata = (;),
               create_output = "")
    create_output != "" && mkpath(create_output)
    data !== nothing && data !== Nothing && present(config, data, metadata)
    for dir in dirs
        add_file_dir(dir)
    end
    if async
        @async server(config, socket)
    else
        server(config, socket)
    end
    config
end

function refresh(config::Config)
    for con in values(config.connections)
        if !con.refresh_queued && !haskey(con.pending_result, :result)
            finish_command(JusCmd(config, con, "refresh"))
        end
    end
end

function present(config::Config, data, metadata = (;))
    var = addvar(config, data)
    #make the connection track the variable so it cleans up on disconnect
    config.init_connection = con-> addvar(con, var)
    @async begin
        try
            sleep(0.5)
            DefaultApplication.open("http://$(config.hostname):$(config.port)/shell.html?$(var.id)")
        catch err
            @error "Error opening browser" exception=err,catch_backtrace()
            exit(1)
        end
    end
end
present(config::Config, type::Type, metadata = (;)) =
    present(config, "$(Base.parentmodule(type)).$type", metadata)
function present(config::Config, type::String, metadata = (;))
    md = join(["$k:$v" for (k, v) in pairs(metadata)], ",")
    !isempty(md) && (md = ",$md")
    @async begin
        try
            sleep(0.5)
            DefaultApplication.open("http://$(config.hostname):$(config.port)/shell.html?@/0:create=$type$md")
        catch err
            @error "Error opening browser" exception=err,catch_backtrace()
            exit(1)
        end
    end
end

changed(cfg::Config, var::Var) = get!(()-> Dict(), cfg.changes, var.id)[:set] = true
function changed(cfg::Config, var::Var, metaproperties...)
    push!(get!(()-> Set(), get!(()-> Dict(), cfg.changes, var.id), :metadata), metaproperties...)
end

set_var(cmd::VarCommand, value) = set_var(cmd, cmd.var, value)
set_var(cmd::VarCommand, name::Symbol, value) = set_var(cmd, cmd.var[name], value)
function set_var(cmd::VarCommand, var::Var, value)
    var.value = value
    changed(cmd.config, var)
end

delete_var(cmd::VarCommand) = delete_var(cmd, cmd.var)
delete_var(cmd::VarCommand, name::Symbol) = delete_var(cmd, cmd.var[name])
function delete_var(cmd::VarCommand, var::Var)
    if var.parent != EMPTYID
        parent = cmd.config[var.parent]
        if var.name != Symbol("")
            delete!(parent.namedchildren, var.name)
        else
            parent.indexedchildren = filter(c-> c !== var, parent.indexedchildren)
        end
        delete!(cmd.config.vars, var.id)
    end
    changes = get!(()-> Dict(), cmd.config.changes, var.id)
    delete!(changes, :set)
    delete!(changes, :metadata)
    changes[:delete] = true
end

set_metadata(cmd::VarCommand, name::Symbol, value) = set_metadata(cmd, cmd.var, name, value)
set_metadata(cmd::VarCommand, varname::Symbol, name::Symbol, value) =
    set_metadata(cmd, cmd.var[varname], name, value)
function set_metadata(cmd::VarCommand, var::Var, name::Symbol, value::AbstractString)
    if !haskey(var.metadata, name) || get(var.metadata, name, "") != value
        var.metadata[name] = value
        changed(cmd.config, var, name)
    end
end

"""
    route(parent_value, cmd)

Route a command, calling handle_child for each ancestor of the variable, starting with the root.
Returns the command.

See handle() for details on commands.
"""
function route(parent_value, cmd::VarCommand)
    path = cmd.path = []
    cur = cmd.var.parent
    while cur !== EMPTYID
        var = cmd.config[cur]
        push!(path, var)
        cur = var.parent
    end
    handle_route(parent_value, cmd)
end

# implemented recursively so developers can overrride it
# using the parent_value's type
function handle_route(parent_value, cmd::VarCommand)
    cmd.cancel && return cmd
    #println("@@@ HANDLE ROUTE $(parent_value), $(cmd)")
    if isempty(cmd.path)
        #println("@@@ FINAL HANDLE $(parent_value), $(cmd)")
        override(cmd, handle(parent_value, cmd))
    else
        var = cmd.path[1]
        cmd.path = @view cmd.path[2:end]
        cmd = override(cmd, handle_child(var, var.value, parent_value, cmd))
        cmd = override(cmd, handle_route(parent_value, cmd))
        override(cmd, finish_handle_child(var, var.value, parent_value, cmd))
    end
end

override(cmd::VarCommand, result) = result isa VarCommand ? result : cmd

"""
    handle_child(var, value, cmd, path)

Allows ancestor variables to alter or cancel commands.
`handle_child` can replace the given VarCommand by return a different one
See handle() for details on commands.
"""
handle_child(ancestor_var, ancestor, value, cmd) = nothing

"""
    finish_handle_child(var, value, cmd, path)

Allows ancestor variables to alter or cancel commands after they have been processed.
See handle() for details on commands.
"""
finish_handle_child(ancestor_var, ancestor, value, cmd) = nothing

"""
    handle(value, cmd::VarCommand{Command, Arg})

Base-level command routing.
Developers can specialize on value and cmd.

COMMANDS

  :metadata Set metadata for a variable. Arg will be a tuple with the metadata's name.
            initially sent before :set and :create
            NOTE: value is the PARENT VARIABLE'S value

  :set      Change the value of a variable. Initially called before :create.
            NOTE: value is the PARENT VARIABLE'S value

  :create   The variable has just been created (called after :metadata and :set)
            NOTE: value is the VARIABLE'S value

  :get      Determine the new value for a variable. By default, variables retain their values
            but handlers can change this behavior.
            NOTE: value is the PARENT VARIABLE'S value
"""
handle(value, cmd::VarCommand) = default_handle(value, cmd)

default_handle(value, cmd::VarCommand{:metadata}) = nothing

function default_handle(value, cmd::VarCommand{:metadata, (:create,)})
    cmd.creating && create_from_metadata(cmd)
end

function default_handle(value, cmd::VarCommand{:metadata, (:path,)})
    @debug("@@@ PATH METADATA: $(repr(cmd))")
    set_path_from_metadata(cmd)
end

function default_handle(value, cmd::VarCommand{:metadata, (:access,)})
    println("@@@ ACCESS METADATA: $(repr(cmd))")
    set_access_from_metadata(cmd.var)
end

default_handle(value, cmd::VarCommand{:metadata, (:observe,)}) =
    push!(cmd.connection.observing, cmd.var.id)

default_handle(value, cmd::VarCommand{:metadata, (:genview,)}) = generate_viewdef(cmd)

function default_handle(value, cmd::VarCommand{:metadata, (:convert,)})
    m = match(CONVERT_ENUM, cmd.var.metadata[:convert])
    if m !== nothing && match(JPATH_COMPONENT, m[1]) !== nothing && !endswith(m[1], "()")
        # safe to eval because m[1] is just "ident(.ident)*"
        values = Dict(lowercase(string(inst)) => inst for inst in instances(Main.eval(Meta.parse(m[1]))))
        cmd.var.value_conversion = v-> values[lowercase(v)]
        println("###\n### CREATED CONVERTER $(cmd.var.metadata[:convert]) FOR $(cmd.var)")
    else
        throw("Could not create converter $(cmd.var.metadata[:convert]) for $(cmd.var)")
    end
end

function default_handle(value, cmd::VarCommand{:set})
    #println("@@@ BASIC SET $(cmd.var.id): $(repr(cmd))")
    #println("@@@@@@ VAR: $(cmd.var)")
    if isempty(cmd.var.path)
        cmd.var.internal_value = cmd.var.value = cmd.arg
        cmd.creating && changed(cmd.config, cmd.var)
        set_type(cmd)
    elseif cmd.var.writeable
        println("@@@@@@@@ SET PATH: $(cmd.var.path) OF $(cmd.var)")
        set_path(cmd)
    end
    !cmd.cancel && !cmd.var.action && changed(cmd.config, cmd.var)
end

function default_handle(value, cmd::VarCommand{:get})
    if has_path(cmd.var)
        get_path(cmd)
        if haskey(cmd.var.metadata, :transformerId)
            cur = cmd.var
            while cur.parent != EMPTYID
                cur = cmd.config[cur.parent]
                if haskey(cmd.var.metadata, :transformer)
                    cmd.var.value = transform(cur, cur.internal_value, cmd)
                    break
                end
            end
        end
    end
end

default_handle(value, cmd::VarCommand{:observe}) = nothing

default_handle(value, cmd::VarCommand{:create}) = nothing

transform(parent_var::Var, parent, cmd::VarCommand) = cmd.var.internal_value

end


text:
using Revise
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

code = @dice begin
    function uniform(b::Int, w::Int) # uniform [0, b)
        num_b = 1 * (b == 1) + (ndigits(b, base = 2) - 1) * (b != 1)
        x = Vector(undef, num_b)
        for i = num_b:-1:1
            x[i] = if b != 1 flip(0.5) else flip(0) end
        end

        ans = DistInt(x)
        return add_bits(ans, w)
    end

    function uniform2(i::Int)
        num_b = ndigits(i, base = 2)
        bits = Vector{Bool}(undef, num_b)
        for bit_idx = 1:num_b
            b::Bool = i & 1
            @inbounds bits[bit_idx] = b
            i = i >> 1
        end
        sum = 0
        sum_v = Vector(undef, num_b)
        pos = 1
        println(bits)
        for i = 1:num_b
            println(bits[i])
            sum+=2^(i-1) * bits[i]
            sum_v[pos] = sum
            pos += 1
        end
        println(sum_v)
        ans = (DistInt([flip(0)]))
        println(num_b)
        for i = 2:num_b
            println(i)
            ans = if flip((2^(i-1) * bits[i])/sum_v[i])
                (uniform(2^(i-1), 3) + (sum_v[i] - 2^(i-1)*bits[i]))[1] else (ans) end
        end
        return ans
    end
    uniform2(63)
end

# BDD analysis
bdd = compile(code)
num_flips(bdd)
num_nodes(bdd)
infer(code, :bdd)

# IR analysis
# to_dice_ir(code)
# has_dice_binary() && rundice(code)
# has_dice_binary() && infer(code, :ocaml)


text:
"""
    function sigma(para, extK = DiagTree.getK(para.totalLoopNum, 1), subdiagram = false; name = :Σ, resetuid = false)

    Build sigma diagram.
    When sigma is created as a subdiagram, then no Fock diagram is generated if para.filter contains NoFock, and no sigma diagram is generated if para.filter contains Girreducible

# Arguments
- `para`            : parameters. It should provide internalLoopNum, interactionTauNum, firstTauIdx
- `extK`            : basis of external loop.
- `subdiagram`      : a sub-vertex or not
- `name`            : name of the diagram
- `resetuid`        : restart uid count from 1

# Output
- A DataFrame with fields `:type`, `:extT`, `:diagram`, `:hash`
- All sigma share the same incoming Tau index, but not the outgoing one
"""
function sigma(para, extK=DiagTree.getK(para.totalLoopNum, 1), subdiagram=false; name=:Σ, resetuid=false)
    resetuid && uidreset()
    @assert para.diagType == SigmaDiag
    @assert para.innerLoopNum >= 1
    # @assert length(extK) == para.totalLoopNum
    # @assert (para.innerLoopNum <= 1) || ((NoBubble in para.filter) == false) "The many-body correction sigma only accounts for half of the bubble counterterm right now."
    if (para.innerLoopNum > 1) && (NoBubble in para.filter)
        @warn "Sigma with two or more loop orders still contain bubble subdiagram even if NoBubble is turned on in para.filter!"
    end

    @assert length(extK) >= para.totalLoopNum "expect dim of extK>=$(para.totalLoopNum), got $(length(extK))"
    extK = extK[1:para.totalLoopNum]

    if isValidSigma(para.filter, para.innerLoopNum, subdiagram) == false
        return DataFrame(type=[], extT=[], diagram=[])
    end

    if (para.extra isa ParquetBlocks) == false
        parquetblocks = ParquetBlocks(phi=[PPr, PHEr], ppi=[PHr, PHEr], Γ4=[PPr, PHr, PHEr])
        para = reconstruct(para, extra=parquetblocks)
    end

    K = zero(extK)
    LoopIdx = para.firstLoopIdx
    K[LoopIdx] = 1.0
    @assert (K ≈ extK) == false
    legK = [extK, K, K, extK]

    function GWwithGivenExTtoΣ(group, oW, paraG)
        # println(group)
        # @assert length(group[:, :diagram]) == 1
        # allsame(group, [:response, :type, :GT])
        @assert group[:response] == UpUp || group[:response] == UpDown
        #type: Instant or Dynamic
        response, type = group[:response], group[:type]
        sid = SigmaId(para, type, k=extK, t=group[:extT])
        g = green(paraG, K, group[:GT], true; name=oW == 0 ? :Gfock : :G_Σ) #there is only one G diagram for a extT
        @assert g isa Diagram
        # Sigma = G*(2 W↑↑ - W↑↓)
        # ! The sign of ↑↓ is from the spin symmetry, not from the fermionic statistics!
        spinfactor = (response == UpUp) ? 2 : -1
        # spinfactor = (response == UpUp) ? 0 : 1
        if oW > 0 # oW are composte Sigma, there is a symmetry factor 1/2
            spinfactor *= 0.5
            # elseif oW == 0 # the Fock diagram requires an additional minus sign, because the interaction is currently a direct one, but is expected to be exchange.
            #     spinfactor *= paraG.isFermi ? -1.0 : 1.0
        end
        # plot_tree(mergeby(DataFrame(group)), maxdepth = 7)
        sigmadiag = Diagram(sid, Prod(), [g, group[:diagram]], factor=spinfactor, name=name)
        # plot_tree(sigmadiag, maxdepth = 7)
        return (type=type, extT=group[:extT], diagram=sigmadiag)
    end

    compositeSigma = DataFrame()

    for (oG, oW) in orderedPartition(para.innerLoopNum - 1, 2, 0)

        idx, maxLoop = findFirstLoopIdx([oW, oG], LoopIdx + 1)
        @assert maxLoop <= para.totalLoopNum "maxLoop = $maxLoop > $(para.totalLoopNum)"
        WfirstLoopIdx, GfirstLoopIdx = idx

        # it is important to do W first, because the left in of W is also the incoming leg of sigma, they have the same Tidx
        idx, maxTau = findFirstTauIdx([oW, oG], [Ver4Diag, GreenDiag], para.firstTauIdx, para.interactionTauNum)
        @assert maxTau <= para.totalTauNum
        WfirstTauIdx, GfirstTauIdx = idx

        paraG = reconstruct(para, diagType=GreenDiag, innerLoopNum=oG,
            firstLoopIdx=GfirstLoopIdx, firstTauIdx=GfirstTauIdx)
        paraW = reconstruct(para, diagType=Ver4Diag, innerLoopNum=oW,
            firstLoopIdx=WfirstLoopIdx, firstTauIdx=WfirstTauIdx)

        #TODO: add validation for paraW
        if isValidG(paraG)
            if oW == 0 # Fock-type Σ
                paraW0 = reconstruct(paraW, filter=union(paraW.filter, Proper), transferLoop=zero(K))
                ver4 = vertex4(paraW0, legK, [], true)
            else # composite Σ
                # paraW0 = reconstruct(paraW, filter=union(paraW.filter, Proper), transferLoop=extK-K)
                ver4 = vertex4(paraW, legK, [PHr,], true, phi_toplevel=[], Γ4_toplevel=[PHr, PHEr, PPr,])
                # plot_tree(mergeby(ver4).diagram[1])
            end
            #transform extT coloum intwo extT for Σ and extT for G
            # plot_tree(ver4)
            df = transform(ver4, :extT => ByRow(x -> [(x[INL], x[OUTR]), (x[OUTL], x[INR])]) => [:extT, :GT])

            groups = mergeby(df, [:response, :type, :GT, :extT], operator=Sum())
            for mergedVer4 in eachrow(groups)
                push!(compositeSigma, GWwithGivenExTtoΣ(mergedVer4, oW, paraG))
            end
        end
    end

    if isempty(compositeSigma)
        return DataFrame(type=[], extT=[], diagram=[])
    end

    Factor = 1 / (2π)^para.loopDim
    sigmadf = mergeby(compositeSigma, [:type, :extT], name=name, factor=Factor,
        getid=g -> SigmaId(para, g[1, :type], k=extK, t=g[1, :extT]))

    @assert all(x -> x[1] == para.firstTauIdx, sigmadf.extT) "all sigma should share the same in Tidx\n$sigmadf"
    # println(sigmadf)
    # plot_tree(sigmadf)
    return sigmadf
end


text:
# Medici et al. Western blot data
# BACKGROUND SUBTRACTED VALUES FROM RAG

# Figure 1B # Ojb1

Fig1B = [
12,207,18, #Snail # 1
12,18,207, #Slug  # 4
1,182,191, #LEF-1 # 7
212,209,224 #a-tubulin #11
]
#Snail,-,+,-,,
#Slug,-,-,+,,

# Figure 1F # Ojb2

Fig1F = [
208,1,1,184,    #E-cadherin # 1
1,197,203,32,   #Vimentin  # 5
1,180,176,30,   #Fibronectin # 9
190,203,202,199 #a-tubulin # 13
]
#Snail,-,'+,-,'+,
#Slug,-,-,'+,-,
#DN LEF-1,-,-,-,'+,

# Figure 2A # Ojb3
Fig2A = [
11,225,221,30, #TGF-b3
209,216,212 #a-tubulin
]
#Snail,-,+,-,fig3F(+ECAD),
#Slug,-,-,+,fig3F(+ECAD),

#Figure 2E # Ojb4
Fig2E = [
1,184,184,18,21, #LEF-1 # 1
169,1,1,173,170, #E-cadherin # 6
1,182,193,36,39, #Vimentin   # 11
1,178,178,38,39, #Fibronectin # 16
201,201,201,201,201, #a-tubulin
]
#Snail,-,'+,-,'+,-
#Slug,-,-,'+,-,'+
#TGF-b3 Ab,-,-,-,'+,'+

#Figure 3D # Ojb5
Fig3D = [
1,141,164,1,1, #b-catenin #1
149,149,149,149,149 #TCF-4 #6
]
#Snail,-,'+,-,'+,-
#Slug,-,-,'+,-,'+
#E-cadherin,-,-,-,'+,'+

#Figure 3F # Ojb6
Fig3F = [
17,188,190,30,32, #TGF-b3 # 1
172,1,1,185,189, #E-cadherin # 6
1,154,159,1,1, #Vimentin #  11
1,160,158,1,1, #Fibronectin #
184,184,184,184,184 #a-tubulin #
]
#Snail,-,'+,-,'+,-
#Slug,-,-,'+,-,'+
#E-cadherin,-,-,-,'+,'+

#Figure 4D # Ojb7

Fig4D = [
10,25, #Snail
5,10, #Slug
217,213 #a-tubulin
]
#TGF-b3,-,'+,,,

#Figure 4E # Ojb8
Fig4E = [
1,186,11, #LEF-1
216,229,223 #a-tubulin
]
#TGF-b3,-,'+,'+,,
#DN Smad4,-,-,'+,,

#Figure 4G # Ojb9
Fig4G = [
192,1,199,183, #E-cadherin # 1
1,182,27,26, #Vimentin     # 5
1,183,9,6, #Fibronectin    # 9
222,228,198,206 #a-tubulin #
]
#TGF-b3,-,'+,'+,'+,
#DN Smad4,-,-,'+,-,
#DN LEF-1,-,-,-,'+,

#Figure 5A # Ojb10
Fig5A = [
48,211,22,209,35, #Snail
40,213,20,211,21, #Slug
204,220,213,224,209 #a-tubulin
]
#TGF-b1,-,'+,'+,-,-
#TGF-b2,-,'+,-,'+,'+
#U0126,-,-,'+,-,'+

#Figure 5C # # Ojb11 # NOT FIGURE 5 C
Fig5C = [
1,217,20,200, #LEF-1
196,199,210,214 #a-tubulin
]

#TGF-b1,-,'+,'+,-,
#TGF-b2,-,'+,-,'+,
#MEKINH,-,-,'+,-,
# WARNING according to paper, 3rd condition is actually TGFB3_AB
# different from RAG excel sheet
# THE OBJECTIVE IS FORMULATED CORRECTLY BUT IT DOES NOT RELATE TO FIGURE 5C WHICH IS NOT TRAIED TO

Fig4F = [
56,224,31 #GSK3_P
]
#TGF-b3,-,+,-
#PI3Kin,-,-,+


text:
function EllipsoidTessellator(arg0::OneAxisEllipsoid, arg1::TileAiming, arg2::jint)
    return EllipsoidTessellator((OneAxisEllipsoid, TileAiming, jint), arg0, arg1, arg2)
end

function build_simple_zone(::Type{EllipsoidTessellator}, arg0::jdouble, arg1::Vector{GeodeticPoint})
    return jcall(EllipsoidTessellator, "buildSimpleZone", SphericalPolygonsSet, (jdouble, Vector{GeodeticPoint}), arg0, arg1)
end

function build_simple_zone(::Type{EllipsoidTessellator}, arg0::jdouble, arg1::Vector{Vector{jdouble}})
    return jcall(EllipsoidTessellator, "buildSimpleZone", SphericalPolygonsSet, (jdouble, Vector{Vector{jdouble}}), arg0, arg1)
end

function sample(obj::EllipsoidTessellator, arg0::SphericalPolygonsSet, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "sample", List, (SphericalPolygonsSet, jdouble, jdouble), arg0, arg1, arg2)
end

function tessellate(obj::EllipsoidTessellator, arg0::SphericalPolygonsSet, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jboolean, arg6::jboolean)
    return jcall(obj, "tessellate", List, (SphericalPolygonsSet, jdouble, jdouble, jdouble, jdouble, jboolean, jboolean), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end