text:
function fillblockconfig(stockdata,n)
    i = findall(stockdata.nUPRNs_undersameroof.==n)

    # Find UPRN groups for N-dwellings per footprint
    UPRN_groups = unique(stockdata.UPRNs_undersameroof[i,:])

    # Find arrays of indices representing each UPRN groups
    findall_arr(arr,val) = findall(arr.==val)
    j = findall_arr.((stockdata.UPRNs_undersameroof,),UPRN_groups)
    # Find dwelling types for each group (if available)
    index_arr(arr,j) = arr[j]
    all_dwellingtypes = index_arr.((stockdata.EPC_Dwelling_Type,),j)

    # Find if dwellings are either flats or missings
    isflat(dwlgtypes) = isempty(findall((occursin.("flat",skipmissing(dwlgtypes)).+occursin.("maisonette",skipmissing(dwlgtypes))).==false))
    simpleflats = DataFrame(:i_block => findall(isflat.(all_dwellingtypes)))

    # Check area tolerence
    simpleflats[!,:arearatio] = index_arr.((stockdata.arearatio,),j)[simpleflats.i_block]
    simpleflats[!,:EPC_area] = index_arr.((stockdata.EPC_Total_floor_area_m2_,),j)[simpleflats.i_block]
    # Keep flats with areas in tolerence
    function areatol(simpleflats)
        i_keep = intersect(
            findall(minimum.(simpleflats.arearatio) .> 0.5),
            findall(maximum.(simpleflats.arearatio) .< 1.2),
        )
        simpleflats = simpleflats[i_keep,:]
    end
    simpleflats = areatol(simpleflats)

    # Add EPC_Dwelling_Type
    simpleflats[!,:EPC_Dwelling_Type] = index_arr.((stockdata.EPC_Dwelling_Type,),j)[simpleflats.i_block]

    # Estimate Dwelling_Type where EPCs are missing
    function eval_dwellingtypes(block,n)
        # Prep block composition "Top" comes before "Ground" to ensure every block
        # ... has a roof (mid-floor can exist in 2-UPRN blocks, possible mixed-use)
        eval_blocktemplate(n) = vcat([["Top-floor ","Ground-floor "],repeat(["Mid-floor "],maximum([0,n-2]))]...)
        block_template = eval_blocktemplate(n)

        # Find any flats in block with missing EPC
        i_missing = findall(ismissing.(block))

        # Remove already occurring flat type(s) from template
        occursin_arr(val,arr) = occursin.(val,arr)
        if isempty(i_missing)
            block_template = []
        elseif !isempty(block[Not(i_missing)])
            i_rm = vcat(collect(1:sum(occursin.("Top-floor ",block[Not(i_missing)])).+0),
                collect(2:sum(occursin.("Ground-floor ",block[Not(i_missing)])).+1),
                collect(3:sum(occursin.("Mid-floor ",block[Not(i_missing)])).+2),
            )
            # `i_rm[findall(i_rm.<=n)]` necessary, as mid-floor can exist in 2-UPRN blocks
            block_template = block_template[Not(i_rm[findall(i_rm.<=n)])]
        end
        [block[i_missing[i]]="$(block_template[i])flat (ParaDwell:ASSUMED)" for i in 1:length(i_missing)]

        # Apply assumptions from block_template for missing dwelling types
        return block
    end
    assumed_dwellingtypes = eval_dwellingtypes.(simpleflats.EPC_Dwelling_Type,n)

    # Map back to master dataframe
    i_master = vcat(j[simpleflats.i_block]...)
    stockdata.EPC_Dwelling_Type[i_master] .= vcat(assumed_dwellingtypes...)
    return stockdata
end


text:
#using Lint

function readSMAT(filename::AbstractString)
    (rows,header) = readdlm(filename;header=true)
    A = sparse(
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),1])+1,
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),2])+1,
               convert(Array{Int64,1},rows[1:parse(Int,header[3]),3]),
               parse(Int,header[1]),
               parse(Int,header[2])
               )
    return A
end

function load_matrix_network(name::AbstractString)
    pathname = joinpath(Pkg.dir("MatrixNetworks"),"data")
    smatfile = joinpath(pathname,"$(name).smat")
    if isfile(smatfile)
        return readSMAT(smatfile)
    else
        error(@sprintf "The example datafile '%s' does not seem to exist where it should\n" name)
    end
end

function load_matrix_network_metadata(name::AbstractString)
    pathname = joinpath(Pkg.dir("MatrixNetworks"),"data")
    smatfile = joinpath(pathname,"$(name).smat")
    meta_xy = joinpath(pathname,"$(name).xy")
    meta_labels = joinpath(pathname,"$(name).labels")
    if isfile(smatfile)
        if isfile(meta_xy) && isfile(meta_labels)
            xy = readdlm(meta_xy)
            labels = readdlm(meta_labels)
            return (readSMAT(smatfile),xy,labels)
        end
    else
        error(@sprintf "The example datafile '%s' does not seem to exist where it should\n" name)
    end
end

function matrix_network_datasets()
    datasets_location = joinpath(Pkg.dir("MatrixNetworks"),"data")
    content = readdir(datasets_location)
    smat_files = filter(x->contains(x,".smat"),content)
    for i = 1:length(smat_files)
        # no need for this lintpragma anymore
        # @lintpragma( "Ignore use of undeclared variable end")
        smat_files[i] = smat_files[i][1:end-5]
    end
    return smat_files
end

ext: jl


text:
#!/usr/bin/env julia

using JCheck
using Test

import Aqua

Aqua.test_all(JCheck, ambiguities = false)

include("types_with_generate.jl")

@time begin
    @testset "JCheck's tests" begin
        include("Quickcheck.jl")
        include("shrink.jl")
        include("generate.jl")
    end
end


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

"""
    LeaveBallOut(ball)

Leave-`ball`-out (a.k.a. spatial leave-one-out) validation.

    LeaveBallOut(radius)

By default, use Euclidean ball of given `radius`.

## References

* Le Rest et al. 2014. Spatial leave-one-out
  cross-validation for variable selection in
  the presence of spatial autocorrelation.
"""
struct LeaveBallOut{B<:BallNeighborhood} <: AbstractErrorEstimator
  ball::B
end

LeaveBallOut(radius::Real) = LeaveBallOut(BallNeighborhood(radius))

function Base.error(solver::AbstractLearningSolver,
                    problem::LearningProblem,
                    eestimator::LeaveBallOut)
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  ball  = eestimator.ball

  # efficient neighborhood search
  searcher = NeighborhoodSearcher(sdata, ball)

  # pre-allocate memory for coordinates
  coords = MVector{ndims(sdata),coordtype(sdata)}(undef)

  solutions = pmap(1:npoints(sdata)) do i
    coordinates!(coords, sdata, i)

    # points inside and outside ball
    inside  = search(coords, searcher)
    outside = [j for j in 1:npoints(sdata) if j ∉ inside]

    # setup and solve learning sub-problem
    subproblem = LearningProblem(view(sdata, outside),
                                 view(sdata, [i]),
                                 task(problem))
    solve(subproblem, solver)
  end

  result = pmap(ovars) do var
    𝔏 = defaultloss(sdata[1,var])
    ŷ = [solutions[i][1,var] for i in 1:npoints(sdata)]
    y = [sdata[i,var] for i in 1:npoints(sdata)]
    var => 𝔏(ŷ, y)
  end

  Dict(result)
end


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

mutable struct PublicIPAddressDnsSettings2 <: SwaggerModel
    domainNameLabel::Any # spec type: Union{ Nothing, String } # spec name: domainNameLabel
    fqdn::Any # spec type: Union{ Nothing, String } # spec name: fqdn
    reverseFqdn::Any # spec type: Union{ Nothing, String } # spec name: reverseFqdn

    function PublicIPAddressDnsSettings2(;domainNameLabel=nothing, fqdn=nothing, reverseFqdn=nothing)
        o = new()
        validate_property(PublicIPAddressDnsSettings2, Symbol("domainNameLabel"), domainNameLabel)
        setfield!(o, Symbol("domainNameLabel"), domainNameLabel)
        validate_property(PublicIPAddressDnsSettings2, Symbol("fqdn"), fqdn)
        setfield!(o, Symbol("fqdn"), fqdn)
        validate_property(PublicIPAddressDnsSettings2, Symbol("reverseFqdn"), reverseFqdn)
        setfield!(o, Symbol("reverseFqdn"), reverseFqdn)
        o
    end
end # type PublicIPAddressDnsSettings2

const _property_map_PublicIPAddressDnsSettings2 = Dict{Symbol,Symbol}(Symbol("domainNameLabel")=>Symbol("domainNameLabel"), Symbol("fqdn")=>Symbol("fqdn"), Symbol("reverseFqdn")=>Symbol("reverseFqdn"))
const _property_types_PublicIPAddressDnsSettings2 = Dict{Symbol,String}(Symbol("domainNameLabel")=>"String", Symbol("fqdn")=>"String", Symbol("reverseFqdn")=>"String")
Base.propertynames(::Type{ PublicIPAddressDnsSettings2 }) = collect(keys(_property_map_PublicIPAddressDnsSettings2))
Swagger.property_type(::Type{ PublicIPAddressDnsSettings2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PublicIPAddressDnsSettings2[name]))}
Swagger.field_name(::Type{ PublicIPAddressDnsSettings2 }, property_name::Symbol) =  _property_map_PublicIPAddressDnsSettings2[property_name]

function check_required(o::PublicIPAddressDnsSettings2)
    true
end

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


text:
# # Comparison with MathOptInterface on a Probability Simplex

# In this example, we project a random point onto a probability simplex with the Frank-Wolfe algorithm using
# either the specialized LMO defined in the package or a generic LP formulation using `MathOptInterface.jl` (MOI) and
# `GLPK` as underlying LP solver.
# It can be found as Example 4.4 [in the paper](https://arxiv.org/abs/2104.06675).

using FrankWolfe

using LinearAlgebra
using LaTeXStrings

using Plots

using JuMP
const MOI = JuMP.MOI

import GLPK

n = Int(1e3)
k = 10000

xpi = rand(n);
total = sum(xpi);
const xp = xpi ./ total;

f(x) = norm(x - xp)^2
function grad!(storage, x)
    @. storage = 2 * (x - xp)
    return nothing
end

lmo_radius = 2.5
lmo = FrankWolfe.FrankWolfe.ProbabilitySimplexOracle(lmo_radius)

x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n))
gradient = collect(x00)

x_lmo, v, primal, dual_gap, trajectory_lmo = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    collect(copy(x00)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=false,
    trajectory=true,
);

# Create a MathOptInterface Optimizer and build the same linear constraints:

o = GLPK.Optimizer()
x = MOI.add_variables(o, n)

for xi in x
    MOI.add_constraint(o, xi, MOI.GreaterThan(0.0))
end

MOI.add_constraint(
    o,
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
    MOI.EqualTo(lmo_radius),
)

lmo_moi = FrankWolfe.MathOptLMO(o)

x, v, primal, dual_gap, trajectory_moi = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo_moi,
    collect(copy(x00)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=false,
    trajectory=true,
);

# Alternatively, we can use one of the modelling interfaces based on `MOI` to formulate the LP.
# The following example builds the same set of constraints using `JuMP`:

m = JuMP.Model(GLPK.Optimizer)
@variable(m, y[1:n] ≥ 0)

@constraint(m, sum(y) == lmo_radius)

lmo_jump = FrankWolfe.MathOptLMO(m.moi_backend)

x, v, primal, dual_gap, trajectory_jump = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo_jump,
    collect(copy(x00)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=false,
    trajectory=true,
);

x_lmo, v, primal, dual_gap, trajectory_lmo_blas = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x00,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.blas,
    verbose=false,
    trajectory=true,
);

x, v, primal, dual_gap, trajectory_jump_blas = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo_jump,
    x00,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.blas,
    verbose=false,
    trajectory=true,
);

iteration_list = [[x[1] + 1 for x in trajectory_lmo], [x[1] + 1 for x in trajectory_moi]]
time_list = [[x[5] for x in trajectory_lmo], [x[5] for x in trajectory_moi]]
primal_gap_list = [[x[2] for x in trajectory_lmo], [x[2] for x in trajectory_moi]]
dual_gap_list = [[x[4] for x in trajectory_lmo], [x[4] for x in trajectory_moi]]

label = [L"\textrm{Closed-form LMO}", L"\textrm{MOI LMO}"]

FrankWolfe.plot_results(
    [primal_gap_list, primal_gap_list, dual_gap_list, dual_gap_list],
    [iteration_list, time_list, iteration_list, time_list],
    label,
    ["", "", L"\textrm{Iteration}", L"\textrm{Time}"],
    [L"\textrm{Primal Gap}", "", L"\textrm{Dual Gap}", ""],
    xscalelog=[:log, :identity, :log, :identity],
    yscalelog=[:log, :log, :log, :log],
    legend_position=[:bottomleft, nothing, nothing, nothing]
)


text:
"""Julia program to check if a number is automorphic or not.
An Automorphic Number is such a number whose square ends in the same digits as  the number itself."""

function check_automophic(num)
    temp = num * num

    while(num > 0)
        # Extract the last digits of the given number and its square.
        lastNum = num % 10
        lastSquare = temp % 10

        # Check if they are equal
        if(lastNum != lastSquare)
            return "The given number $n is not an Automorphic Number."
        end

        num = num ÷ 10
        temp = temp ÷ 10
    end

    return "The given number $n is an Automorphic Number."
end

print("Enter the number: ")
n = readline()
n = parse(Int, n)
res = check_automophic(abs(n))
print(res)

"""
Time Complexity: O(log(n)), where 'n' is the given number
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE 1

Enter the number: 24
The given number 24 is not an Automorphic Number.

SAMPLE 2

Enter the number: 25
The given number 25 is an Automorphic Number.
"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
module CATKEVerticalDiffusivities

using Adapt
using KernelAbstractions: @kernel, @index

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Utils
using Oceananigans.Fields
using Oceananigans.Fields: ZeroField
using Oceananigans.BoundaryConditions: default_prognostic_bc, DefaultBoundaryCondition
using Oceananigans.BoundaryConditions: BoundaryCondition, FieldBoundaryConditions, DiscreteBoundaryFunction, FluxBoundaryCondition
using Oceananigans.BuoyancyModels: ∂z_b, top_buoyancy_flux
using Oceananigans.Operators: ℑzᵃᵃᶜ

using Oceananigans.TurbulenceClosures:
    getclosure,
    time_discretization,
    AbstractTurbulenceClosure,
    AbstractScalarDiffusivity,
    ConvectiveAdjustmentVerticalDiffusivity,
    ExplicitTimeDiscretization,
    VerticallyImplicitTimeDiscretization,
    ThreeDimensionalFormulation,
    VerticalFormulation

import Oceananigans.BoundaryConditions: getbc
import Oceananigans.Utils: with_tracers
import Oceananigans.TurbulenceClosures:
    validate_closure,
    add_closure_specific_boundary_conditions,
    calculate_diffusivities!,
    DiffusivityFields,
    implicit_linear_coefficient,
    viscosity,
    diffusivity,
    diffusive_flux_x,
    diffusive_flux_y,
    diffusive_flux_z

function hydrostatic_turbulent_kinetic_energy_tendency end

struct CATKEVerticalDiffusivity{TD, CD, CL, CQ} <: AbstractScalarDiffusivity{TD, VerticalFormulation}
    Cᴰ :: CD
    mixing_length :: CL
    surface_TKE_flux :: CQ
end

function CATKEVerticalDiffusivity{TD}(Cᴰ:: CD,
                                      mixing_length :: CL,
                                      surface_TKE_flux :: CQ) where {TD, CD, CL, CQ}

    return CATKEVerticalDiffusivity{TD, CD, CL, CQ}(Cᴰ, mixing_length, surface_TKE_flux)
end

const CATKEVD{TD} = CATKEVerticalDiffusivity{TD} where TD

# Support for "ManyIndependentColumnMode"
const CATKEVDArray{TD} = AbstractArray{<:CATKEVD{TD}} where TD
const FlavorOfCATKE{TD} = Union{CATKEVD{TD}, CATKEVDArray{TD}} where TD

function with_tracers(tracer_names, closure::FlavorOfCATKE)
    :e ∈ tracer_names ||
        throw(ArgumentError("Tracers must contain :e to represent turbulent kinetic energy " *
                            "for `CATKEVerticalDiffusivity`."))

    return closure
end

# For tuples of closures, we need to know _which_ closure is CATKE.
# Here we take a "simple" approach that sorts the tuple so CATKE is first.
# This is not sustainable though if multiple closures require this.
# The two other possibilities are:
# 1. Recursion to find which closure is CATKE in a compiler-inferrable way
# 2. Store the "CATKE index" inside CATKE via validate_closure.
validate_closure(closure_tuple::Tuple) = Tuple(sort(collect(closure_tuple), lt=catke_first))

catke_first(closure1, catke::FlavorOfCATKE) = false
catke_first(catke::FlavorOfCATKE, closure2) = true
catke_first(closure1, closure2) = false

# Two CATKEs?!
catke_first(catke1::FlavorOfCATKE, catke2::FlavorOfCATKE) = error("Can't have two CATKEs in one closure tuple.")

include("mixing_length.jl")
include("surface_TKE_flux.jl")
include("turbulent_kinetic_energy_equation.jl")

for S in (:MixingLength, :SurfaceTKEFlux)
    @eval @inline convert_eltype(::Type{FT}, s::$S) where FT = $S{FT}(; Dict(p => getproperty(s, p) for p in propertynames(s))...)
    @eval @inline convert_eltype(::Type{FT}, s::$S{FT}) where FT = s
end

"""
    CATKEVerticalDiffusivity(time_discretization = VerticallyImplicitTimeDiscretization, FT=Float64;
                             Cᴰ = 0.215,
                             mixing_length = MixingLength{FT}(),
                             surface_TKE_flux = SurfaceTKEFlux{FT}(),
                             warning = true)

Returns the `CATKEVerticalDiffusivity` turbulence closure for vertical mixing by
small-scale ocean turbulence based on the prognostic evolution of subgrid
Turbulent Kinetic Energy (TKE).
"""
CATKEVerticalDiffusivity(FT::DataType; kw...) = CATKEVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), FT; kw...)

function CATKEVerticalDiffusivity(time_discretization = VerticallyImplicitTimeDiscretization(), FT=Float64;
                                  Cᴰ = 0.81,
                                  mixing_length = MixingLength{FT}(),
                                  surface_TKE_flux = SurfaceTKEFlux{FT}(),
                                  warning = true)

    if warning
        @warn "CATKEVerticalDiffusivity is an experimental turbulence closure that \n" *
              "is unvalidated and whose default parameters are not calibrated for \n" *
              "realistic ocean conditions or for use in a three-dimensional \n" *
              "simulation. Use with caution and report bugs and problems with physics \n" *
              "to https://github.com/CliMA/Oceananigans.jl/issues."
    end

    Cᴰ = convert(FT, Cᴰ)
    mixing_length = convert_eltype(FT, mixing_length)
    surface_TKE_flux = convert_eltype(FT, surface_TKE_flux)

    return CATKEVerticalDiffusivity{typeof(time_discretization)}(Cᴰ, mixing_length, surface_TKE_flux)
end

#####
##### Diffusivities and diffusivity fields utilities
#####

function DiffusivityFields(grid, tracer_names, bcs, closure::FlavorOfCATKE)

    default_diffusivity_bcs = (Kᵘ = FieldBoundaryConditions(grid, (Center, Center, Face)),
                               Kᶜ = FieldBoundaryConditions(grid, (Center, Center, Face)),
                               Kᵉ = FieldBoundaryConditions(grid, (Center, Center, Face)))

    bcs = merge(default_diffusivity_bcs, bcs)

    Kᵘ = CenterField(grid, boundary_conditions=bcs.Kᵘ)
    Kᶜ = CenterField(grid, boundary_conditions=bcs.Kᶜ)
    Kᵉ = CenterField(grid, boundary_conditions=bcs.Kᵉ)
    Lᵉ = CenterField(grid) #, boundary_conditions=nothing)

    # Secret tuple for getting tracer diffusivities with tuple[tracer_index]
    _tupled_tracer_diffusivities         = NamedTuple(name => name === :e ? Kᵉ : Kᶜ          for name in tracer_names)
    _tupled_implicit_linear_coefficients = NamedTuple(name => name === :e ? Lᵉ : ZeroField() for name in tracer_names)

    return (; Kᵘ, Kᶜ, Kᵉ, Lᵉ, _tupled_tracer_diffusivities, _tupled_implicit_linear_coefficients)
end

@inline viscosity_location(::FlavorOfCATKE) = (Center(), Center(), Face())
@inline diffusivity_location(::FlavorOfCATKE) = (Center(), Center(), Face())

function calculate_diffusivities!(diffusivities, closure::FlavorOfCATKE, model)

    arch = model.architecture
    grid = model.grid
    velocities = model.velocities
    tracers = model.tracers
    buoyancy = model.buoyancy
    clock = model.clock
    top_tracer_bcs = NamedTuple(c => tracers[c].boundary_conditions.top for c in propertynames(tracers))

    event = launch!(arch, grid, :xyz,
                    calculate_CATKE_diffusivities!,
                    diffusivities, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    return nothing
end

@kernel function calculate_CATKE_diffusivities!(diffusivities, grid, closure::FlavorOfCATKE, args...)
    i, j, k, = @index(Global, NTuple)

    # Ensure this works with "ensembles" of closures, in addition to ordinary single closures
    closure_ij = getclosure(i, j, closure)

    @inbounds begin
        diffusivities.Kᵘ[i, j, k] = Kuᶜᶜᶠ(i, j, k, grid, closure_ij, args...)
        diffusivities.Kᶜ[i, j, k] = Kcᶜᶜᶠ(i, j, k, grid, closure_ij, args...)
        diffusivities.Kᵉ[i, j, k] = Keᶜᶜᶠ(i, j, k, grid, closure_ij, args...)
        diffusivities.Lᵉ[i, j, k] = implicit_dissipation_coefficient(i, j, k, grid, closure_ij, args...)
    end
end

@inline function implicit_linear_coefficient(i, j, k, grid, closure::FlavorOfCATKE{<:VITD}, K, ::Val{id}, args...) where id
    L = K._tupled_implicit_linear_coefficients[id]
    return @inbounds L[i, j, k]
end

@inline turbulent_velocity(i, j, k, grid, e) = @inbounds sqrt(max(zero(eltype(grid)), e[i, j, k]))

@inline function Kuᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    u★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocity, tracers.e)
    ℓu = momentum_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    return ℓu * u★
end

@inline function Kcᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    u★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocity, tracers.e)
    ℓc = tracer_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    return ℓc * u★
end

@inline function Keᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    u★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocity, tracers.e)
    ℓe = TKE_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    return ℓe * u★
end

@inline viscosity(::FlavorOfCATKE, diffusivities) = diffusivities.Kᵘ
@inline diffusivity(::FlavorOfCATKE, diffusivities, ::Val{id}) where id = diffusivities._tupled_tracer_diffusivities[id]

#####
##### Show
#####

function Base.summary(closure::CATKEVD)
    TD = nameof(typeof(time_discretization(closure)))
    return string("CATKEVerticalDiffusivity{$TD}")
end

Base.show(io::IO, closure::FlavorOfCATKE) = print(io, summary(closure))

end # module


text:
using Test, MathematicalPredicates

@time @testset "Basic functionality" begin include("basic.jl") end
@time @testset "LazySets support" begin include("LazySets.jl") end


text:
handle("updateeditor") do data
    @destruct [
        text || "",
        mod || "Main",
        path || nothing,
        updateSymbols || true
    ] = data

    return try
        updateeditor(text, mod, path, updateSymbols)
    catch err
        OutlineItem[]
    end
end

# NOTE: update outline and symbols cache all in one go
function updateeditor(text, mod = "Main", path = nothing, updateSymbols = true)
    # update symbols cache
    # ref: https://github.com/JunoLab/Juno.jl/issues/407
    updateSymbols && updatesymbols(mod, path, text)

    # return outline
    outline(toplevelitems(text))
end

struct OutlineItem
    name::String
    type::String
    icon::String
    start::Int
    stop::Int
end
OutlineItem(name, type, icon, item::ToplevelItem) =
    OutlineItem(name, type, icon, first(item.lines), last(item.lines))

outline(items)::Vector{OutlineItem} = filter!(item -> item !== nothing, outlineitem.(items))

outlineitem(item::ToplevelItem) = nothing # fallback case
outlineitem(binding::ToplevelBinding) = begin
    expr = binding.expr
    bind = binding.bind
    name = CSTParser.has_sig(expr) ? str_value(CSTParser.get_sig(expr)) : bind.name
    OutlineItem(name, static_type(bind), static_icon(bind), binding)
end
outlineitem(call::ToplevelCall) = begin
    # show includes
    return if isinclude(call.expr)
        OutlineItem(call.verbatim, "module", "icon-file-code", call)
    elseif isprecompile(call.expr)
        OutlineItem(call.verbatim, "module", "icon-file-binary", call)
    else
        nothing
    end
end
outlineitem(macrocall::ToplevelMacroCall) = begin
    # don't show doc strings
    isdoc(macrocall.expr) && return nothing
    # show first line verbatim of macro calls
    firstline = strlimit(first(split(macrocall.verbatim, '\n')), 100)
    OutlineItem(firstline, "snippet", "icon-mention", macrocall)
end
outlineitem(usage::ToplevelModuleUsage) = begin
    useline = replace(str_value(usage.expr), ":" => ": ")
    OutlineItem(useline, "mixin", "icon-package", usage)
end


text:
include(joinpath("Thesis model", "Mos model.jl"))
import .Mos_Model

res_mo = []
for var in [0.01:0.05:5;]
    push!(res_mo, Mos_Model.mos_model(var))
end

################################################################

using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Thesis model", "Figures", "EMP_Model")
resultspath = joinpath("Thesis model", "Results", "EMP_Model")

res_df = Array{String,Float64}(res_emp)
df = DataFrame[]
for i = 1:100
    push!(df, res_emp[i])
end
CSV.write(joinpath(resultspath, "Results_EMP.csv"), res_emp)

x = [results_emp["mu"] for results_emp in res_emp]
a = [results_emp["g6p"] for results_emp in res_emp]
b = [results_emp["glu"] for results_emp in res_emp]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)

ax.xlabel = "Biomass [1/h]"
ax.ylabel = "Metabolite concentration"

fig[1, 1] = Legend(
    fig,
    [[scatter1, line1], [scatter2, line2]],
    ["Glucose-6-Phosphate", "Glutamate"],
    tellheight = false,
    tellwidth = false,
    halign = :right,
    valign = :center,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Glu_EMP.pdf"), fig)


text:
using Butterfly
using CSV
using DataFrames
using Statistics
using BenchmarkTools
using NPZ
using ProgressMeter

simdata = CSV.read(joinpath(@__DIR__, "benchmark_data.csv"))

datadir = joinpath(@__DIR__, "data")
mkpath(datadir)

times = DataFrame(time = Float64[], nallocs = Int[], datafile = String[])
@showprogress for (i, row) in enumerate(eachrow(simdata))
    b = @benchmarkable simulate($row)
    tune!(b)
    t, df = BenchmarkTools.run_result(b)
    dfile = joinpath(datadir, "julia_$i.npy")
    npzwrite(dfile, df)
    push!(times, [mean(t).time * 1e-9, t.allocs, dfile])
end

CSV.write(joinpath(@__DIR__, "julia_times.csv"), times)

max_stars_repo_name: zclaytor/butterpy

max_issues_repo_name: zclaytor/butterpy

max_forks_repo_name: zclaytor/butterpy


text:
# workspace()

# include("./types.jl")
#  Global variables are bad style.
# g_num_features = 0 # Int64
# g_num_examples = 0 # Int64
# g_the_normalizer = 0.0  # Float64
# g_is_initialized = false
#
# # Initialize some global variables.
# function InitializeTreeData(examples::Vector{Example}, normalizer::Float64)
# # leave it here, but avoid using it
#   global g_num_examples = length(examples)
#   # length of feature, num_examples should >= 1
#   global g_num_features = length(examples[1].values)
#   global g_the_normalizer = normalizer # the global normalizer
#   global g_is_initialized = true # globally initialized
#   return (g_num_examples, g_num_features, g_the_normalizer, g_is_initialized)
# end

# Return root node for a tree.
function MakeRootNode(examples::Vector{Example})
  # return Node root
  root = Node()
  root.examples = examples
  root.positive_weight = root.negative_weight = 0
  for example in examples
    if example.label == 1
      root.positive_weight += example.weight
    else   # label == -1
      root.negative_weight += example.weight
    end
  end
  root.leaf = true
  root.depth = 0
  return root
end

# Return a tree trained on examples.
function TrainTree(examples::Vector{Example}, beta::Float64, lambda::Float64,
                  tree_depth::Int64, normalizer::Float64)
  tree = Tree() # vector of Nodes
  push!(tree, MakeRootNode(examples))
  node_id = Int64(1) # julia starts from 1
  # println("length(tree) is ", length(tree))
  # for node_id in 1:length(tree)
  while node_id <= length(tree) # by yy: for the previous range is static
    # reference: any problems?
    # fine: need to be updated in the end, by function return
    node = tree[node_id]
    best_split_feature = Int64(0)
    best_split_value = 0.
    best_delta_gradient = 0.

    # following three are global in cpp, we need to calculate here
    num_features = length(examples[1].values) # length of feature, num_examples should >= 1
    num_examples = length(examples) # only used in ComplexityPenalty
    the_normalizer = normalizer # the normalizer, only used in ComplexityPenalty

    # search for the best split feature and the best split value
    for split_feature in 1:num_features
      value_to_weights = MakeValueToWeightsMap(node, split_feature)
      # this is a little different from the cpp code,
      # since cpp achieves changes of multiple by pointers in the input params/args
      (split_value, delta_gradient) = BestSplitValue(value_to_weights, node,
                                          length(tree), beta, lambda,
                                          the_normalizer,
                                          num_examples, num_features)
      if delta_gradient > best_delta_gradient + kTolerance
        best_delta_gradient = delta_gradient
        best_split_feature = split_feature
        best_split_value = split_value
      end
    end
    # add nodes to the tree if the depth doesn't exceed the limit
    # and the best delta gradient is not zero
    if (node.depth < tree_depth) && (best_delta_gradient > kTolerance)
      node, tree = MakeChildNodes(best_split_feature, best_split_value, node, tree)
    end
    # println("Run for ", node_id, " iterations to add the nodes.")
    node_id += 1
  end # end of train a new tree
  return tree
end

# Given a value-to-weights map for a feature (constructed by MakeValueToWeightsMap()),
# determine the best split value for the feature and
# the improvement in the gradient of the objective if we split on that value.
# Note that delta_gradient <= 0 indicates that we should not split on this feature.
function BestSplitValue(value_to_weights::Dict{Float64, Vector{Float64}},
                    node::Node, tree_size::Int64, beta::Float64, lambda::Float64,
                    the_normalizer::Float64, num_examples::Int64, num_features::Int64)
# used in TrainTree
# split_value: float pointer, no need to be in the input params
# delta_gradient: float pointer, no need to be in input params
# in julia, there is no pointer, so we return them to store the changed values
  split_value = 0. # by yy: split_value should be initialized, not sure if 0 is
                   # reasonable
  delta_gradient = 0.
  left_positive_weight = 0.
  left_negative_weight = 0.
  right_positive_weight = node.positive_weight
  right_negative_weight = node.negative_weight

  old_error = min(left_positive_weight + right_positive_weight,
                  left_negative_weight + right_negative_weight)
  old_gradient = Gradient(old_error, tree_size, 0., -1, beta, lambda,
                          the_normalizer, num_examples, num_features)
  for elem in sort(collect(value_to_weights), by=x->x[1]) # elem now is (key, value) Tuple
    left_positive_weight  += elem[2][1]
    right_positive_weight -= elem[2][1]
    left_negative_weight  += elem[2][2]
    right_negative_weight -= elem[2][2]

    new_error = min(left_positive_weight, left_negative_weight) +
                min(right_positive_weight, right_negative_weight)
    new_gradient = Gradient(new_error, tree_size + 2, 0., -1, beta, lambda,
                            the_normalizer, num_examples, num_features)
    if abs(new_gradient) - abs(old_gradient) > delta_gradient + kTolerance
      delta_gradient = abs(new_gradient) - abs(old_gradient)
      split_value = elem[1]
    end
  end
  return (split_value, delta_gradient)
end

# Return the (sub)gradient of the objective with respect to a tree.
function Gradient(wgtd_error::Float64, tree_size::Int64, alpha::Float64,
                  sign_edge::Int64, beta::Float64, lambda::Float64,
                  the_normalizer::Float64, num_examples::Int64,
                  num_features::Int64)
  # TODO(usyed): Can we make some mild assumptions and get rid of sign_edge?
  # const float
  complexity_penalty = ComplexityPenalty(tree_size, beta, lambda,
                                         the_normalizer,
                                         num_examples, num_features)
  edge = wgtd_error - 0.5 # const float
  sign_alpha = (alpha >= 0) ? 1 : -1 # const int
  if abs(alpha) > kTolerance
    return edge + sign_alpha * complexity_penalty
  elseif abs(edge) <= complexity_penalty
    return 0
  else
    return edge - sign_edge * complexity_penalty
  end
end

# Return complexity penalty.
function ComplexityPenalty(tree_size::Int64, beta::Float64, lambda::Float64,
                  the_normalizer::Float64, num_examples::Int64, num_features::Int64)
# we need to handle these global variables
# the_normalizer
# num_examples
# num_features
  rademacher = sqrt((((2. * tree_size+1.) * (log(num_features+2.)/log(2.))*log(num_examples)) / num_examples))
  return ((lambda*rademacher+beta)*num_examples) / (2*the_normalizer)
end

# Return a map from each value of feature to a pair of weights. The first
# weight in the pair is the total weight of positive examples at node that have
# that value for feature, and the second weight in the pair is the total weight
# of negative examples at node that have that value for feature. This map is
# used to determine the best split feature/value.
# map<Value::Float64, pair<Weight::Float64, Weight::Float64>>
function MakeValueToWeightsMap(node::Node, feature::Int64)
  # map<Value, pair<Weight, Weight>>
  # modified here by yy: for tuples are immutable
  value_to_weights = Dict{Float64, Vector{Float64}}()
  for example in node.examples
    if ~haskey(value_to_weights, example.values[feature])
      value_to_weights[example.values[feature]] = Vector{Float64}([0.,0.])
    end
    if example.label == 1
      value_to_weights[example.values[feature]][1] += example.weight
    else   # label = -1
      value_to_weights[example.values[feature]][2] += example.weight
    end
  end
  return value_to_weights
end

# Make child nodes using split feature/value and add them to the tree. Also
# update info in the parent node, like child pointers.
function MakeChildNodes(split_feature::Int64, split_value::Float64, parent::Node,
                    tree::Tree)
# parent: pointer Node
# tree: pointer Tree
  parent.split_feature = split_feature
  parent.split_value = split_value
  parent.leaf = false
  # println(tree[1].leaf) # To check whether the tree has been modified
  # Node left_child, right_child
  left_child = Node()
  right_child = Node()
  # Initialize the variables in child nodes
  left_child.examples = Array{Example, 1}()
  right_child.examples = Array{Example, 1}()
  left_child.depth = right_child.depth = parent.depth + 1
  left_child.leaf = right_child.leaf = true
  left_child.positive_weight = 0.
  left_child.negative_weight = 0.
  right_child.positive_weight = 0.
  right_child.negative_weight = 0.

  for example in parent.examples
    # Node* child, how to handle this?
    child = Node()
    if example.values[split_feature] <= split_value
      child = left_child
    else
      child = right_child
    end
    # TODO(usyed): Moving examples around is inefficient.
    push!(child.examples, example)
    if example.label == 1
      child.positive_weight += example.weight
    else #   label == -1
      child.negative_weight += example.weight
    end
  end
  # update the info of parent
  # by yy: modified to adapt to julia style
  parent.left_child_id = length(tree) + 1
  parent.right_child_id = length(tree) + 2
  push!(tree, left_child)
  push!(tree, right_child)
  return (parent, tree)
end

# Given an example and a tree, classify the example with the tree.
# NB: This function assumes that if an example has a feature value that is
# _less than or equal to_ a node's split value then the example should be sent
# to the left child, and otherwise sent to the right child.
function ClassifyExampleTree(example::Example, tree::Tree)
# return label of an example
# check the tree
  # for i in 1:length(tree)
  #   println("child info of node ", i, " is: left -> ", tree[i].left_child_id,
  #           " right -> ", tree[i].right_child_id,
  #           " leaf ? -> ", tree[i].leaf )
  # end
  node = tree[1] # 0-->1; const Node pointer, fine to do in this way in julia
  while node.leaf == false
    if example.values[node.split_feature] <= node.split_value
      node = tree[node.left_child_id]
    else
      node = tree[node.right_child_id]
    end
    # println(node.leaf)
  end
  # println("Complete reaching the leaf")
  if node.positive_weight >= node.negative_weight
    return 1
  else
    return -1
  end
end

# Given a set of examples and a tree,
# return the weighted error of tree on the examples.
function EvaluateTreeWgtd(examples::Vector{Example}, tree::Tree)
  wgtd_error = 0.
  for example in examples
    # println("...")
    if ClassifyExampleTree(example, tree) != example.label
      wgtd_error += example.weight
    end
  end
  return wgtd_error
end


text:
"""
    cosine_similarity(u::AbstractVector, v::AbstractVector) -> Float

Compute the cosine similarity of two vectors.

# Examples

```jldoctest
julia> cosine_similarity([1, 0], [1, 0])
1.0

julia> round(cosine_similarity([1,1], [1,0]), digits=3)
0.707

julia> cosine_similarity([1, 0], [0, 1])
0.0

julia> cosine_similarity([-1, 0], [1, 0])
-1.0

julia> round(cosine_similarity([1,0,0], [1,1,1]), digits=3)
0.577
```

"""
function cosine_similarity(u::AbstractVector, v::AbstractVector)
    u ⋅ v / (norm(u) * norm(v))
end

"""
    angle_between(u::AbstractVector, v::AbstractVector) -> Float

Compute the angle between two vectors. The angle is returned in radians.

# Examples

```jldoctest
julia> angle_between([1, 0], [1, 0])
0.0

julia> round(angle_between([1,1], [1,0]), digits=3)
0.785

julia> round(angle_between([1, 0], [0, 1]), digits=3)
1.571

julia> round(angle_between([-1, 0], [1, 0]), digits=3)
3.142

julia> round(angle_between([1,0,0], [1,1,1]), digits=3)
0.955
```

"""
function angle_between(u::AbstractVector, v::AbstractVector)
    acos(cosine_similarity(u, v))
end

"""
    distance(point_a::AbstractVector, point_b::AbstractVector) -> Float

Compute the distance between two points.

# Examples

```jldoctest
julia> distance([0, 0], [0, 0])
0.0

julia> distance([1, 0], [0, 0])
1.0

julia> round(distance([1, 1], [2, 2]), digits=3)
1.414

julia> round(distance([0, 0, 0], [-1, -1, -1]), digits=3)
1.732
```

"""
function distance(point_a::AbstractVector, point_b::AbstractVector)
    return norm(Vector(point_a, point_b))
end

"""
    distance(point_a::AbstractVector, line::AbstractLine) -> Float

Compute the distance from a point to a line.

This is the distance from the point to its projection on the line.

# Examples

```jldoctest
julia> distance([0, 0], Line([0, 0], [1, 0]))
0.0

julia> round(distance([1, 0], Line([0, 0], [1, 1])), digits=3)
0.707

julia> round(distance([1, 2, 3], Line([-1, 3, 2], [7, 4, 2])), digits=3)
1.978
```

"""
function distance(point::AbstractVector, line::AbstractLine)
    point_projected = project(point, line)
    return distance(point, point_projected)
end


text:
module JacobianTests

using DiffTests, ForwardDiff, ReverseDiff, Test

include(joinpath(dirname(@__FILE__), "../utils.jl"))

function test_unary_jacobian(f, x)
    test_val = f(x)
    test = ForwardDiff.jacobian!(DiffResults.JacobianResult(test_val, x), f, x, ForwardDiff.JacobianConfig(f, x))

    # without JacobianConfig

    test_approx(ReverseDiff.jacobian(f, x), DiffResults.jacobian(test))

    out = similar(DiffResults.jacobian(test))
    ReverseDiff.jacobian!(out, f, x)
    test_approx(out, DiffResults.jacobian(test))

    result = DiffResults.JacobianResult(test_val, x)
    ReverseDiff.jacobian!(result, f, x)
    test_approx(DiffResults.value(result), DiffResults.value(test))
    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))

    # with JacobianConfig

    cfg = ReverseDiff.JacobianConfig(x)

    test_approx(ReverseDiff.jacobian(f, x, cfg), DiffResults.jacobian(test))

    out = similar(DiffResults.jacobian(test))
    ReverseDiff.jacobian!(out, f, x, cfg)
    test_approx(out, DiffResults.jacobian(test))

    result = DiffResults.JacobianResult(test_val, x)
    ReverseDiff.jacobian!(result, f, x, cfg)
    test_approx(DiffResults.value(result), DiffResults.value(test))
    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))

    # with JacobianTape

    tp = ReverseDiff.JacobianTape(f, rand(eltype(x), size(x)))

    test_approx(ReverseDiff.jacobian!(tp, x), DiffResults.jacobian(test))

    out = similar(DiffResults.jacobian(test))
    ReverseDiff.jacobian!(out, tp, x)
    test_approx(out, DiffResults.jacobian(test))

    result = DiffResults.JacobianResult(test_val, x)
    ReverseDiff.jacobian!(result, tp, x)
    test_approx(DiffResults.value(result), DiffResults.value(test))
    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))

    # with compiled JacobianTape

    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
        ctp = ReverseDiff.compile(tp)

        test_approx(ReverseDiff.jacobian!(ctp, x), DiffResults.jacobian(test))

        out = similar(DiffResults.jacobian(test))
        ReverseDiff.jacobian!(out, ctp, x)
        test_approx(out, DiffResults.jacobian(test))

        result = DiffResults.JacobianResult(test_val, x)
        ReverseDiff.jacobian!(result, ctp, x)
        test_approx(DiffResults.value(result), DiffResults.value(test))
        test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
    end
end

function test_unary_jacobian(f!, y, x)
    y_original = copy(y)
    y_copy = copy(y)
    test = ForwardDiff.jacobian!(DiffResults.JacobianResult(y_copy, x), f!, y_copy, x)

    # without JacobianConfig

    out = ReverseDiff.jacobian(f!, y, x)
    test_approx(y, DiffResults.value(test))
    test_approx(out, DiffResults.jacobian(test))
    copyto!(y, y_original)

    out = similar(DiffResults.jacobian(test))
    ReverseDiff.jacobian!(out, f!, y, x)
    test_approx(y,   DiffResults.value(test))
    test_approx(out, DiffResults.jacobian(test))
    copyto!(y, y_original)

    result = DiffResults.JacobianResult(y, x)
    ReverseDiff.jacobian!(result, f!, y, x)
    @test DiffResults.value(result) == y
    test_approx(y, DiffResults.value(test))
    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
    copyto!(y, y_original)

    # with JacobianConfig

    cfg = ReverseDiff.JacobianConfig(y, x)

    out = ReverseDiff.jacobian(f!, y, x, cfg)
    test_approx(y,   DiffResults.value(test))
    test_approx(out, DiffResults.jacobian(test))
    copyto!(y, y_original)

    out = similar(DiffResults.jacobian(test))
    ReverseDiff.jacobian!(out, f!, y, x, cfg)
    test_approx(y,   DiffResults.value(test))
    test_approx(out, DiffResults.jacobian(test))
    copyto!(y, y_original)

    result = DiffResults.JacobianResult(y, x)
    ReverseDiff.jacobian!(result, f!, y, x, cfg)
    @test DiffResults.value(result) == y
    test_approx(y, DiffResults.value(test))
    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
    copyto!(y, y_original)

    # with JacobianTape

    tp = ReverseDiff.JacobianTape(f!, y, rand(eltype(x), size(x)))

    out = ReverseDiff.jacobian!(tp, x)
    test_approx(out, DiffResults.jacobian(test))

    out = similar(DiffResults.jacobian(test))
    ReverseDiff.jacobian!(out, tp, x)
    test_approx(out, DiffResults.jacobian(test))

    result = DiffResults.JacobianResult(y, x)
    ReverseDiff.jacobian!(result, tp, x)
    test_approx(DiffResults.value(result), DiffResults.value(test))
    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))

    # with compiled JacobianTape

    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
        ctp = ReverseDiff.compile(tp)

        out = ReverseDiff.jacobian!(ctp, x)

        test_approx(out, DiffResults.jacobian(test))

        out = similar(DiffResults.jacobian(test))
        ReverseDiff.jacobian!(out, ctp, x)
        test_approx(out, DiffResults.jacobian(test))

        result = DiffResults.JacobianResult(y, x)
        ReverseDiff.jacobian!(result, ctp, x)
        test_approx(DiffResults.value(result), DiffResults.value(test))
        test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
    end
end

function test_binary_jacobian(f, a, b)
    test_val = f(a, b)
    test_a = ForwardDiff.jacobian(x -> f(x, b), a)
    test_b = ForwardDiff.jacobian(x -> f(a, x), b)

    # without JacobianConfig

    Ja, Jb = ReverseDiff.jacobian(f, (a, b))
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    Ja = similar(a, length(test_val), length(a))
    Jb = similar(b, length(test_val), length(b))
    ReverseDiff.jacobian!((Ja, Jb), f, (a, b))
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    Ja = DiffResults.JacobianResult(test_val, a)
    Jb = DiffResults.JacobianResult(test_val, b)
    ReverseDiff.jacobian!((Ja, Jb), f, (a, b))
    test_approx(DiffResults.value(Ja), test_val)
    test_approx(DiffResults.value(Jb), test_val)
    test_approx(DiffResults.gradient(Ja), test_a)
    test_approx(DiffResults.gradient(Jb), test_b)

    # with JacobianConfig

    cfg = ReverseDiff.JacobianConfig((a, b))

    Ja, Jb = ReverseDiff.jacobian(f, (a, b), cfg)
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    Ja = similar(a, length(test_val), length(a))
    Jb = similar(b, length(test_val), length(b))
    ReverseDiff.jacobian!((Ja, Jb), f, (a, b), cfg)
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    Ja = DiffResults.JacobianResult(test_val, a)
    Jb = DiffResults.JacobianResult(test_val, b)
    ReverseDiff.jacobian!((Ja, Jb), f, (a, b), cfg)
    test_approx(DiffResults.value(Ja), test_val)
    test_approx(DiffResults.value(Jb), test_val)
    test_approx(DiffResults.jacobian(Ja), test_a)
    test_approx(DiffResults.jacobian(Jb), test_b)

    # with JacobianTape

    tp = ReverseDiff.JacobianTape(f, (rand(eltype(a), size(a)), rand(eltype(b), size(b))))

    Ja, Jb = ReverseDiff.jacobian!(tp, (a, b))
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    Ja = similar(a, length(test_val), length(a))
    Jb = similar(b, length(test_val), length(b))
    ReverseDiff.jacobian!((Ja, Jb), tp, (a, b))
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    Ja = DiffResults.JacobianResult(test_val, a)
    Jb = DiffResults.JacobianResult(test_val, b)
    ReverseDiff.jacobian!((Ja, Jb), tp, (a, b))
    test_approx(DiffResults.value(Ja), test_val)
    test_approx(DiffResults.value(Jb), test_val)
    test_approx(DiffResults.gradient(Ja), test_a)
    test_approx(DiffResults.gradient(Jb), test_b)

    # with compiled JacobianTape

    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
        ctp = ReverseDiff.compile(tp)

        Ja, Jb = ReverseDiff.jacobian!(ctp, (a, b))
        test_approx(Ja, test_a)
        test_approx(Jb, test_b)

        Ja = similar(a, length(test_val), length(a))
        Jb = similar(b, length(test_val), length(b))
        ReverseDiff.jacobian!((Ja, Jb), ctp, (a, b))
        test_approx(Ja, test_a)
        test_approx(Jb, test_b)

        Ja = DiffResults.JacobianResult(test_val, a)
        Jb = DiffResults.JacobianResult(test_val, b)
        ReverseDiff.jacobian!((Ja, Jb), ctp, (a, b))
        test_approx(DiffResults.value(Ja), test_val)
        test_approx(DiffResults.value(Jb), test_val)
        test_approx(DiffResults.gradient(Ja), test_a)
        test_approx(DiffResults.gradient(Jb), test_b)
    end
end

for f in (DiffTests.ARRAY_TO_ARRAY_FUNCS..., DiffTests.MATRIX_TO_MATRIX_FUNCS...)
    test_println("ARRAY_TO_ARRAY_FUNCS + MATRIX_TO_MATRIX_FUNCS", f)
    test_unary_jacobian(f, rand(5, 5))
end

for f! in DiffTests.INPLACE_ARRAY_TO_ARRAY_FUNCS
    test_println("INPLACE_ARRAY_TO_ARRAY_FUNCS", f!)
    test_unary_jacobian(f!, rand(25), rand(25))
end

for f in DiffTests.BINARY_MATRIX_TO_MATRIX_FUNCS
    test_println("BINARY_MATRIX_TO_MATRIX_FUNCS", f)
    test_binary_jacobian(f, rand(5, 5), rand(5, 5))
end

for f in (DiffTests.ARRAY_TO_ARRAY_FUNCS..., DiffTests.MATRIX_TO_MATRIX_FUNCS...)
    test_println("ARRAY_TO_ARRAY_FUNCS + MATRIX_TO_MATRIX_FUNCS", f)

    x = rand(5, 5)
    test = ForwardDiff.jacobian(y -> ForwardDiff.jacobian(f, y), x)

    # without JacobianTape

    J = ReverseDiff.jacobian(y -> ReverseDiff.jacobian(f, y), x)
    test_approx(J, test)

    # with JacobianTape

    tp = ReverseDiff.JacobianTape(y -> ReverseDiff.jacobian(f, y), rand(eltype(x), size(x)))
    J = ReverseDiff.jacobian!(tp, x)
    test_approx(J, test)
end

for f in DiffTests.BINARY_MATRIX_TO_MATRIX_FUNCS
    test_println("BINARY_MATRIX_TO_MATRIX_FUNCS", f)

    a, b = rand(5, 5), rand(5, 5)

    test_val = f(a, b)
    test_a = ForwardDiff.jacobian(y -> ForwardDiff.jacobian(x -> f(x, b), y), a)
    test_b = ForwardDiff.jacobian(y -> ForwardDiff.jacobian(x -> f(a, x), y), b)

    # without JacobianTape

    Ja = ReverseDiff.jacobian(y -> ReverseDiff.jacobian(x -> f(x, b), y), a)
    Jb = ReverseDiff.jacobian(y -> ReverseDiff.jacobian(x -> f(a, x), y), b)
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    # with JacobianTape

    ra = ReverseDiff.JacobianTape(y -> ReverseDiff.jacobian(x -> f(x, b), y), rand(eltype(a), size(a)))
    rb = ReverseDiff.JacobianTape(y -> ReverseDiff.jacobian(x -> f(a, x), y), rand(eltype(b), size(b)))
    Ja = ReverseDiff.jacobian!(ra, a)
    Jb = ReverseDiff.jacobian!(rb, b)
    test_approx(Ja, test_a)
    test_approx(Jb, test_b)

    # The below will fail until support for the Jacobian of
    # functions with multiple output arrays is implemented

    # Ja, Jb = ReverseDiff.jacobian((x, y) -> ReverseDiff.jacobian(f, (x, y)), (a, b))
    # test_approx(Ja test_a)
    # test_approx(Jb test_b)
end

end # module


text:
export rbspsurf

"""
UNTESTED	rbspsurf(b[], ordx, ordy, npts, mpts, p1, p2)

Calculate a Cartesian product rational B-spline surface (NURBS) using an open uniform knot vector.

Call: basis, knot, sumrbas

---

# Arguments
- `B::Array{Float64}`: .
- `ordx::Int64`: order in the x direction
- `ordy::Int64`: order in the y direction
- `npts::Int64`: the number of control net vertices in x direction
- `mpts::Int64`: the number of control net vertices in y direction
- `p1::Int64`: number of parametric lines in the x direction
- `p2::Int64`: number of parametric lines in the y direction

---

# Examples
```jldoctest
julia>
```

```jldoctest
julia>
```

```jldoctest
julia>
```

---

_By Paolo Macciacchera, Elia Onofri_
"""

function rbspsurf(B::Array{Float64,2}, ordx::Int64, ordy::Int64, npts::Int64, mpts::Int64, p1::Int64, p2::Int64)::Tuple{Array{Float64, 2}, Array{Array{Int64, 1}, 1}, Array{Array{Int64, 1}, 1}}

    nplusc = npts + ordx
    mplusc = mpts + ordy
    x = zeros(nplusc)
    y = zeros(mplusc)
    nbasis = zeros(1, npts)
    mbasis = zeros(1, mpts)
    q = zeros(3, p1 * p2)

    #generate the open uniform knot vectors
    x = knot(npts, ordx)
    y = knot(mpts, ordy)

    icount = 0
    #calculate the points on the B-spline surface
    stepu = x[nplusc] / (p1 - 1)
    stepw = y[mplusc] / (p2 - 1)
    for u = 0 : stepu : x[nplusc]
        nbasis = rbasis(npts, ordx, u, x, B[4,:])
        for w = 0 : stepw : y[mplusc]
            mbasis = rbasis(mpts, ordy, w, y, B[4,:])
            sum = sumrbas(B, nbasis, mbasis, npts, mpts)
            icount = icount + 1
            for i = 1 : npts
                for j = 1 : mpts
		    if mbasis[j]!=0
                        for s = 1 : 3
                            j1 = mpts * (i - 1) + j
                            qtemp = (B[4, j1] * B[s, j1] * nbasis[i] * mbasis[j]) / sum
                            q[s, icount] = q[s, icount] + qtemp
                        end
		    end
                end
            end
        end
    end

    EV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for i = 1 : p1 * p2 - 1
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+1
        push!(EV,C)
    end

    for i = 1 : p2 * (p1 - 1)
        C = Array{Int64}(2)
        C[1] = i
        C[2] = i+p2
        push!(EV,C)
    end

    FV = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

    for j = 0 : p1 - 2
        for i = 1 : p2 - 1
            C = Array{Int64}(4)
            C[1] = j * p2 + i
            C[2] = j * p2 + i + 1
            C[3] = (j + 1) * p2 + i
            C[4] = (j + 1) * p2 + i + 1
            push!(FV,C)
        end
    end
    return(q, EV, FV)
end

#-----------------------------------------------------------------------

"""
UNTESTED	sumrbas(B[], nbasis, mbasis, npts, mpts)

Calculate the sum of the nonrational basis functions.

---

# Arguments
- `B::Array{Float64}`: array containing the control net vertices
- `nbasis::Array{Float64}`: array containing the nonrational basis functions for x
- `mbasis::Array{Float64}`: array containing the nonrational basis functions for y
- `npts::Int64`: the number of control net vertices in x direction
- `mpts::Int64`: the number of control net vertices in y direction

---

# Examples
```jldoctest
julia>
```

```jldoctest
julia>
```

```jldoctest
julia>
```

---

_By Paolo Macciacchera, Elia Onofri_
"""

function sumrbas(B::Array{Float64,2}, nbasis::Array{Float64}, mbasis::Array{Float64}, npts::Int64, mpts::Int64)::Float64
    sum = 0
    for i = 1 : npts
        for j = 1 : mpts
            j1 = mpts * (i - 1) + j
            sum = sum + B[4, j1] * nbasis[i] * mbasis[j]
        end
    end
    return(sum)
end


text:
function CLIMA_slicer(solver_config)
    ## Slicing routine for structured CLIMA meshes
    # Define number nodal points along each axis of an element and number of elements along each
    # axis of topology
    Np=[4,4,4]
    Nn=[Np[1]+1,Np[2]+1,Np[3]+1]
    Ne=[20,20,20]

    # Get X, Y, Z coords and make sure they are in a CPU array to make life easier.
    Gnd=reshape(solver_config.dg.grid.vgeo,(Nn[1],Nn[2],Nn[3],16,Ne[1],Ne[2],Ne[3]))
    xvals_raw=Gnd[:,:,:,13,:,:,:]
    yvals_raw=Gnd[:,:,:,14,:,:,:]
    zvals_raw=Gnd[:,:,:,15,:,:,:]
    xvals_cpu=ones(size(xvals_raw))
    yvals_cpu=ones(size(yvals_raw))
    zvals_cpu=ones(size(zvals_raw))
    copyto!(xvals_cpu,xvals_raw)
    copyto!(yvals_cpu,yvals_raw)
    copyto!(zvals_cpu,zvals_raw)

    # Get field of interest and put on CPU too
    Qnd=reshape(solver_config.Q.realdata,(Nn[1],Nn[2],Nn[3],4,Ne[1],Ne[2],Ne[3]));
    iFld=1 # u      ( m/s)
    iFld=2 # v      ( m/s)
    iFld=3 # Eta    (   m)
    iFld=4 # Theta  (  oC)
    phi1_raw=Qnd[:,:,:,iFld,:,:,:]
    phi1_cpu=ones(size(phi1_raw))
    copyto!(phi1_cpu,phi1_raw)

    # Have coords and values in CPU arrays
    # Now stick things in a DataFrame to make manipulation easier
    df = DataFrame(X=xvals_cpu[:],
                   Y=yvals_cpu[:],
                   Z=zvals_cpu[:],
                   F=phi1_cpu[:])
    println(" ")
    println( first(df,6) )
    println(" ")

    # Get underlying coords by finding unique vals to some tolerance
    # Check lengths ( should be Np*Ne+1 )
    deps=1.e-2
    dedup(x,d)=(
     # x must be sorted before this
     for i=2:length(x); dx=x[i]-x[i-1];
      if abs(dx) < d;
       x[i]=x[i-1];
      end;
     end;
     return unique(x);
    )
    lx_exp=Np[1]*Ne[1]+1
    xlocs=sort!( unique(df,:X), [:X]              )[!,:X]
    xlocs=dedup(xlocs,deps)
    println(" ")
    println("xlocs = ",xlocs)
    println(" ")
    @assert length(xlocs) == lx_exp
    ly_exp=Np[2]*Ne[2]+1
    ylocs=sort!( unique(df,:Y), [:Y]              )[!,:Y]
    ylocs=dedup(ylocs,deps)
    println(" ")
    println("ylocs = ",ylocs)
    println(" ")
    @assert length(ylocs) == ly_exp
    lz_exp=Np[3]*Ne[3]+1
    zlocs=sort!( unique(df,:Z), [:Z], rev=(false) )[!,:Z]
    zlocs=dedup(zlocs,deps)
    println(" ")
    println("zlocs = ",zlocs)
    println(" ")
    @assert length(zlocs) == lz_exp

    # Now get slices by finding points that fit mesh within tolerance geps
    geps=1.e-5

    # Example 1.
    # Line in Z with Z descending retaining jump points.
    line_z_val=df[ ( isapprox.( df.X , xlocs[1];atol=geps) ) .&
		   ( isapprox.( df.Y , ylocs[1];atol=geps) ), : ]
    println(" ")
    println(" line_z_val")
    println( sort!(line_z_val, [:Z], rev=(true) ) )
    println(" ")

    s_choice="YZ"

    # XZ slice with X fastest and plot 2d heatmap
    # Set slicing params
    if s_choice == "XZ"
     const_ax  = df.Y
     const_val = ylocs[5]
     xc=xlocs
     yc=zlocs
     slice_2d_val=df[ ( isapprox.( const_ax , const_val;atol=geps) ), : ]
     ax1       = slice_2d_val.X
     xsym      = :X
     ax2       = slice_2d_val.Z
     ysym      = :Z
    end

    # YZ slice with Y fastest and plot 2d heatmap
    # Set slicing params
    if s_choice == "YZ"
     const_ax  = df.X
     const_val = xlocs[5]
     xc=ylocs
     yc=zlocs
     slice_2d_val=df[ ( isapprox.( const_ax , const_val;atol=geps) ), : ]
     ax1       = slice_2d_val.Y
     xsym      = :Y
     ax2       = slice_2d_val.Z
     ysym      = :Z
    end

    # XY slice with X fastest and plot 2d heatmap
    # Set slicing params
    if s_choice == "XY"
     const_ax  = df.Z
     const_val = zlocs[end]
     xc=xlocs
     yc=ylocs
     slice_2d_val=df[ ( isapprox.( const_ax , const_val;atol=geps) ), : ]
     ax1       = slice_2d_val.X
     xsym      = :X
     ax2       = slice_2d_val.Y
     ysym      = :Y
    end

    # Now have a slice specified and with points in native order
    # whatever that may be

    # Normalize column and row coords to fit standard locs in xc and yc
    # - this helps in dealing with jump points, and is also a check
    # - that things are working.
    # x-dimension of slice
    i=0
    for loc in ax1
	    l_diff=loc.-xc
	    i=i+1
	    v,j=findmin(abs.(l_diff))
	    # if we understood grid coorectly all points should be within
	    # geps of a standard loc coordinate value.
	    @assert v <= geps
	    if v != 0
              slice_2d_val[i,xsym]=xc[j]
            end
    end
    # y-dimension of slice
    i=0
    for loc in ax2
	    l_diff=loc.-yc
      	    i=i+1
	    v,j=findmin(abs.(l_diff))
	    # if we understood grid coorectly all points should be within
	    # geps of a standard loc coordinate value.
	    @assert v <= geps
	    if v != 0
              slice_2d_val[i,ysym]=yc[j]
            end
    end

    println(" ")
    println(" slice_2d_val")
    println( first( sort(slice_2d_val, [xsym,ysym], rev=(false,false) ), 400) )
    println(" ")

      # Make a scatter plot - include "jump" points
      # TBD

      # Make a contourf plot - filtering "jump" points.
      # Here we filter by selecting one of them.
      # xc=xlocs
      lx=length(xc)
      # yc=zlocs
      ly=length(yc)
      # Take first jump point at element edges
      # ( Needed for a structured plot since xc and yc are Np*Ne+1, whereas )
      # ( slice_xz_val will be (Np+1)*Ne                                    )
      slice_2d_ct=unique(slice_2d_val,[xsym,ysym])
      # Now need to sort so valuees are ordered with coord
      sort!(slice_2d_ct,[xsym,ysym])
      println( size(slice_2d_ct) )
      zvals=reshape(slice_2d_ct[!,:F],(lx,ly))
      println( lx )
      println( ly )
      # savefig( contour(yc,xc,zvals), "fooo.png" )
      savefig( heatmap(xc,yc,zvals,c=:gist_ncar), "fooo.png" )

      # Save Dataframe to csv - include "jump" points.

   exit()

end


text:
using Test
using SMTPClient

@testset "SMTPClient" begin
  for t ∈ (:send, :error)
    @info "testset: $t..."
    include("./$t.jl")
  end
end


text:
function solve()
    n = parse(Int, readline())
    as = zeros(Int, n)
    bs = zeros(Int, n)
    for i=1:n
        a, b = [parse(Int, x) for x in split(readline())]
        as[i] = a
        bs[i] = b
    end

    ret = 10^6
    for i=1:n, j=1:n
        t = if i == j
            as[i] + bs[j]
        else
            max(as[i], bs[j])
        end
        ret = min(ret, t)
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
type ArterialTerminal # lumped peripheral vasculature
    C::Vector{Float64}
    R::Vector{Float64}
    L::Vector{Float64}
    V0::Vector{Float64}
    P::Any
    V::Any
    Q::Any

    function ArterialTerminal()
        this = new()
        this.C = Vector{Float64}[];
        this.R = Vector{Float64}[];
        this.L = Vector{Float64}[];
        this.V0 = Vector{Float64}[];
        this.P = Array{Float64,2}[];
        this.V = Array{Float64,2}[];
        this.Q = Array{Float64,2}[];
        return this
    end
end

type ArterialBranches # 1D arterial domain
    name::Vector{String}
    parentname::Vector{String}
    ID::Vector{Int64}
    parentID::Vector{Int64}
    lengthincm::Vector{Float64}
    radiusincm::Vector{Float64}
    thicknessincm::Vector{Float64}
    YoungsModinMPa::Vector{Float64}
    A0::Any
    beta::Any
    c0::Any
    k::Vector{Float64}
    group::Vector{String}
    A::Any
    Q::Any
    P::Any
    Fp::Any
    Fbarforward::Any
    Fbarbackward::Any
    Abackward::Any
    Aforward::Any
    Qbackward::Any
    Qforward::Any
    W1end::Vector{Float64}
    W1::Vector{Float64}
    W2::Vector{Float64}
    W1root::Float64
    W2root::Float64
    children::Any
    term::Any
    # Rao::Float64

    function ArterialBranches(filename="test.csv",old=Dict("a"=>0),restart="no")
        this = new()
        this.name = Vector{String}[];
        this.parentname = Vector{String}[];
        this.ID = Vector{Int64}[];
        this.parentID = Vector{Int64}[];
        this.children = Vector{Int64}[];
        this.lengthincm = Vector{Float64}[];
        this.radiusincm = Vector{Float64}[];
        this.thicknessincm = Vector{Float64}[];
        this.YoungsModinMPa = Vector{Float64}[];
        this.A0 = Array{Float64,1}[];
        this.beta = Array{Float64,1}[];
        this.c0 = Array{Float64,1}[];
        this.k = Vector{Float64}[];
        this.group = Vector{String}[];
        this.A = Array{Float64,2}[];
        this.Q = Array{Float64,2}[];
        this.P = Array{Float64,2}[];
        this.Fp = Array{Float64,1}[];
        this.Fbarforward = Array{Float64,1}[];
        this.Fbarbackward = Array{Float64,1}[];
        this.Abackward = Array{Float64,1}[];
        this.Aforward = Array{Float64,1}[];
        this.Qbackward = Array{Float64,1}[];
        this.Qforward = Array{Float64,1}[];
        this.W1end = Vector{Float64}[];
        this.W1 = Vector{Float64}[];
        this.W2 = Vector{Float64}[];
        if restart == "no"
            this.W1root = Float64(0);
            this.W2root = Float64(0);
            # this.Rao = 0.01*mmHgToPa/cm3Tom3;

            # pull in artery data from text file
            temp = CVModule.loadtexttree(filename);
            this.name = [string(temp[1,:Name])]
            this.parentname = [string(temp[1,:ParentName])]
            this.ID = [temp[1,:ID]]
            this.parentID = [temp[1,:parentID]]
            this.lengthincm = [temp[1,:Length_cm]]
            this.radiusincm = [temp[1,:Radius_cm]]
            this.thicknessincm = [temp[1,:Thickness_cm]]
            this.YoungsModinMPa = [temp[1,:YoungsModulus_MPa]]
            for i = 2:length(temp[:Name])
                push!(this.name,string(temp[i,:Name]))
                push!(this.parentname,string(temp[i,:ParentName]))
                push!(this.ID,temp[i,:ID])
                push!(this.parentID,temp[i,:parentID])
                push!(this.lengthincm,temp[i,:Length_cm])
                push!(this.radiusincm,temp[i,:Radius_cm])
                push!(this.thicknessincm,temp[i,:Thickness_cm])
                push!(this.YoungsModinMPa,temp[i,:YoungsModulus_MPa])
            end
            for i in 1:size(temp[:children_1],1)
                if !isa(temp[i,:group],Missings.Missing)
                    push!(this.group,string(temp[i,:group]))
                else
                    push!(this.group,"NA")
                end
                if !isa(temp[i,:children_1],Missings.Missing)
                    push!(this.children,[temp[i,:children_1]])
                else
                    push!(this.children,[])
                end
                if !isa(temp[i,:children_2],Missings.Missing)
                    push!(this.children[i],temp[i,:children_2])
                end
                if !isa(temp[i,:children_3],Missings.Missing)
                    push!(this.children[i],temp[i,:children_3])
                end
                if !isa(temp[i,:children_4],Missings.Missing)
                    push!(this.children[i],temp[i,:children_4])
                end
            end
        elseif restart == "yes"
            this.W1root = old["W1root"];
            this.W2root = old["W2root"];
            this.name = [old["name"][1]];
            this.parentname = [old["parentname"][1]];
            this.ID = [old["ID"][1]];
            this.parentID = [old["parentID"][1]];
            this.group = [old["group"][1]];
            this.lengthincm = [old["lengthincm"][1]];
            this.radiusincm = [old["radiusincm"][1]];
            this.thicknessincm = [old["thicknessincm"][1]];
            this.YoungsModinMPa = [old["YoungsModinMPa"][1]];
            for i = 2:length(old["ID"])
                push!(this.name,old["name"][i]);
                push!(this.parentname,old["parentname"][i]);
                push!(this.ID,old["ID"][i]);
                push!(this.parentID,old["parentID"][i]);
                push!(this.group,old["group"][i]);
                push!(this.lengthincm,old["lengthincm"][i]);
                push!(this.radiusincm,old["radiusincm"][i]);
                push!(this.thicknessincm,old["thicknessincm"][i]);
                push!(this.YoungsModinMPa,old["YoungsModinMPa"][i]);
            end
            for i = 1:length(old["children"])
                temp = Array{Int64}(1);
                if length(old["children"][i]) == 1
                    temp = [old["children"][i]];
                elseif length(old["children"][i]) > 1
                    temp = old["children"][i];
                else
                    temp = [0];
                end
                push!(this.children,temp)
                deleteat!(this.children[i],findin(this.children[i],0))
            end
        end

        this.term = Vector{CVModule.ArterialTerminal}(length(this.ID));

        return this
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using DataFrames
using LinearAlgebra
using PhysicalConstants.CODATA2018
using Unitful

import Unitful:
    nm, μm, mm, cm, ns, μs, ms, ps, s

"""
    set_units(df::DataFrame)

Return a dataframe with the same structure as the input but with length and
time units set.
"""
function set_units(df::DataFrame)
    ## mm columns, assumes structure as 2021-10-05, should be generalised
    len_cols = [:r1, :r2, :r1x, :r2x, :x1, :x2, :xb1, :xb2, :xr1, :xr2,
        :xs, :xt1, :xt2, :y1, :y2, :yb1, :yb2, :yr1, :yr2, :ys, :yt1, :yt2, :z1,
        :z2, :zb1, :zb2, :zr1, :zr2, :zs, :zt1, :zt2]

    # ns cols:
    time_cols = [:t1, :t2, :ta1, :ta2, :tr1, :tr2]
    if :ti1 in propertynames(df)
        ## For backwards compatibility of data runs.
        push!(time_cols, :ti1, :ti2)
    end

    units_df = transform(df, len_cols .=> (x -> x * mm), renamecols=false)
    transform!(units_df, time_cols .=> (x -> x * ns), renamecols=false)
    return units_df
end

"""
	radial_correction(x::Vector{<:Real}, y::Vector{<:Real},
                      z::Vector{<:Real}, r::Vector{<:Real})
Computes the x, y position of the interaction from a radius
of interaction and the x, y position from the SiPM barycentre
Kept for backwards compatibility, z argument pointless
"""
function radial_correction(x::Vector{<:Real}, y::Vector{<:Real},
                  z::Vector{<:Real}, r::Vector{<:Real})
    xy = calculate_xy.(x, y, r)
    xy = reinterpret(reshape, eltype(xy[1]), xy)
    xy[1,:], xy[2,:], z
end

"""
    calculate_xy(x1::Real, y1::Real, r::Real)
Calculate the xy position of a ray at r given x1 and y1 at a different r
"""
function calculate_xy(x1::Real, y1::Real, r::Real)
    ϕ = atan(y1, x1)
    x = r * cos(ϕ)
    y = r * sin(ϕ)
    return x, y
end

"""
    predict_interaction_radius
Given a calibration funcion and the physical limits of the Detector
and function return the radius of interaction (maxr - DOI).
"""
function predict_interaction_radius(cal_func::Function,
        minr::Real, maxr::Real, bias::Real=zero(Real))
  function prediction(sigmas::Vector{<:Real})
    rpred = cal_func.(sigmas) .+ bias
    rpred[rpred .< minr] .= minr
    rpred[rpred .> maxr] .= maxr
    return rpred
  end
end

"""
    time_of_flight(df::DataFrame, ri::Symbol, rsipm::Real, n::Real)

Calculate the time of flight between an interaction radius (ri) and the
detection radius (rsipm).
"""
function time_of_flight(df::DataFrame, ri::Symbol, rsipm::Real, n::Real=1.6)
    cc       = SpeedOfLightInVacuum / n
    path_len = rsipm * mm .- df[!, ri]
    return uconvert.(ps, path_len ./ cc)
end

"""
    time_of_flight(df::DataFrame, source::Vector{Symbol},
                    interaction::Vector{Symbol}, n::Real)

Time of flight between two arbitrary points in a medium of refraction index n.
"""
function time_of_flight(df::DataFrame, source_xyz::Vector{Symbol},
        int_xyz::Vector{Symbol}, n::Real=1.0)
    function displacement(x1, y1, z1, x2, y2, z2)
        norm.(eachrow(hcat(x1, y1, z1) - hcat(x2, y2, z2)))
    end
    cc       = SpeedOfLightInVacuum / n
    path_len = combine(df, vcat(int_xyz, source_xyz) => displacement => :path)
    return uconvert.(ps, path_len[!, :path] ./ cc)
end

"""
    interaction_time(df::DataFrame, ri::Symbol, td::Symbol, rsipm::Real, n::Real)

Calculate the interaciton time given radius of interaction (ri) and detection time.
"""
function interaction_time(df::DataFrame, ri::Symbol, td::Symbol,
        rsipm::Real, n::Real=1.69)
    flight_time = time_of_flight(df, ri, rsipm, n)
    return uconvert.(ps, df[!, td]) .- flight_time
end

"""
    CRT(df::DataFrame, times::Vector{Symbol}, radii::Vector{Symbol}, det_rad::Real)

Calculate the difference between the calculated and expected interaction
time differences.
"""
#TODO: Protections for which symbols valid?
function CRT(df::DataFrame, times::Vector{Symbol},
        radii::Vector{Symbol}, det_rad::Real, n::Real=1.69)::Vector{<:Real}
    int1_flight = time_of_flight(df, [:xs, :ys, :zs], [:xt1, :yt1, :zt1])
    int2_flight = time_of_flight(df, [:xs, :ys, :zs], [:xt2, :yt2, :zt2])
    true_dt     = int2_flight ./ ps .- int1_flight ./ ps

    int1_det    = interaction_time(df, radii[1], times[1], det_rad, n)
    int2_det    = interaction_time(df, radii[2], times[2], det_rad, n)
    calc_dt     = int2_det ./ ps .- int1_det ./ ps

    return calc_dt - true_dt
end

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
################################################################################
# Components - Processors, Memory, IO Handlers
################################################################################

# Task classes are used to control which tasks are mapped to which cores.
# each core has fulfill one or more task classes. Each task will be assigned a
# class and will be mapped to a core that supports that task.
"All task classes currently recognized by the Mapper."
const _mapper_task_classes = Set([
      "processor",
      "memory_processor",
      "input_handler",
      "output_handler",
      "memory_1port",
      "memory_2port",
    ])

# Function to unify mapper types cleanly between tasks and cores.
memory_meta(ports::Int) = "memory_$(ports)port"
ismemory(s::String) = startswith(s, "memory") && endswith(s, "port")
ismemory(x) = false

# Global NamedTuple for dealing with attributes.
const MTypes = (
    proc        = "processor",
    memoryproc  = "memory_processor",
    input       = "input_handler",
    output      = "output_handler",
    memory      = memory_meta,
)

typekey() = "mapper_type"

ismappable(c::AbstractComponent)    = haskey(c.metadata, typekey())

isproc(c::AbstractComponent)        = ismappable(c) && in(MTypes.proc, c.metadata[typekey()])
ismemoryproc(c::AbstractComponent)  = ismappable(c) && in(MTypes.memoryproc, c.metadata[typekey()])
isinput(c::AbstractComponent)       = ismappable(c) && in(MTypes.input, c.metadata[typekey()])
isoutput(c::AbstractComponent)      = ismappable(c) && in(MTypes.output, c.metadata[typekey()])
islowpower(c::AbstractComponent)     = ismappable(c) && in(MTypes.lowpower, c.metadata[typekey()])
ishighperformance(c::AbstractComponent) = ismappable(c) && in(MTypes.highperformance, c.metadata[typekey()])

function ismemory(c::AbstractComponent)
    ismappable(c) || return false
    for i in c.metadata[typekey()]
        if ismemory(i)
            return true
        end
    end
    return false
end

################################################################################
# Assignment of metadata to Components.
################################################################################
proc_metadata() = Dict{String,Any}(typekey() => [MTypes.proc])
mem_proc_metadata() = Dict{String,Any}(typekey() => [MTypes.proc, MTypes.memoryproc])

function mem_nport_metadata(nports::Int)
    attrs = [MTypes.memory(i) for i in 1:nports]
    return Dict{String,Any}(
        typekey() => attrs,
        # Options for plotting
        "shadow_offset" => [Address(0,i) for i in 1:nports-1],
        "fill"          => "memory"
    )
end

input_handler_metadata() = Dict{String,Any}(typekey() => [MTypes.input])
output_handler_metadata() = Dict{String,Any}(typekey() => [MTypes.output])

add_lowpower(c::Component) = push!(c.metadata[typekey()], MTypes.lowpower)
add_highperformance(c::Component) = push!(c.metadata[typekey()], MTypes.highperformance)

################################################################################
# Metadata for Taskgraph Nodes
################################################################################

const TN = TaskgraphNode
# Setting task -> core requirements.
make_input!(t::TN)       = (t.metadata[typekey()] = MTypes.input)
make_output!(t::TN)      = (t.metadata[typekey()] = MTypes.output)
make_proc!(t::TN)        = (t.metadata[typekey()] = MTypes.proc)
make_memoryproc!(t::TN)  = (t.metadata[typekey()] = MTypes.memoryproc)
make_memory!(t::TN, ports::Integer) = t.metadata[typekey()] = MTypes.memory(ports)

# Getting task -> core requirements.
ismappable(t::TN)   = true
isinput(t::TN)      = t.metadata[typekey()] == MTypes.input
isoutput(t::TN)     = t.metadata[typekey()] == MTypes.output
isproc(t::TN)       = t.metadata[typekey()] == MTypes.proc
ismemoryproc(t::TN) = t.metadata[typekey()] == MTypes.memoryproc
ismemory(t::TN)     = ismemory(t.metadata[typekey()])

################################################################################
# Routing Resources
################################################################################

"Routng link types currently implemented in the Mapper."
const _mapper_link_classes = Set([
        "circuit_link",
        "memory_request_link",
        "memory_response_link",
    ])

# Generic function for making a metadata dictionary for ports/links where only
# the link-class is needed.
function routing_metadata(link_class)
    if !in(link_class, _mapper_link_classes)
        throw(KeyError(link_class))
    end

    return Dict{String,Any}(
        "link_class" => link_class,
    )
end

# Processor component metadata.
function proc_output_metadata(ndirs, nlayers)
    metadata_vec = map(1:(ndirs * nlayers)) do i
        return Dict{String,Any}(
            "link_class"    => "circuit_link",
            "index"         => i-1,
            "network_id"    => div(i-1, ndirs),
        )
    end
    return metadata_vec
end

function proc_fifo_metadata(nfifos)
    metadata_vec = map(1:nfifos) do i
        return Dict{String,Any}(
            "link_class" => "circuit_link",
            "index" => i-1
        )
    end
    return metadata_vec
end

function proc_memory_request_metadata()
    return Dict{String,Any}(
        "link_class" => "memory_request_link",
        "index" => 0
    )
end

function proc_memory_return_metadata()
    return Dict{String,Any}(
        "link_class" => "memory_response_link",
        "index" => 0
    )
end

# Memories
function mem_memory_request_metadata(nports)
    metadata_vec = map(1:nports) do i
        return Dict{String,Any}(
            "link_class" => "memory_request_link",
            "index" => i-1,
        )
    end
    return metadata_vec
end

function mem_memory_return_metadata(nports)
    metadata_vec = map(1:nports) do i
        return Dict{String,Any}(
            "link_class" => "memory_response_link",
            "index" => i-1,
        )
    end
    return metadata_vec
end

# io handlers
function input_handler_port_metadata(nlinks)
    metadata_vec = map(1:nlinks) do i
        return Dict{String,Any}(
            "link_class" => "circuit_link",
            "index" => i-1,
        )
    end
end

function output_handler_port_metadata(nlinks)
    metadata_vec = map(1:nlinks) do i
        return Dict{String,Any}(
            "link_class" => "circuit_link",
            "index" => i-1,
        )
    end
    return metadata_vec
end

interpolate(x0, y0, x1, y1, t) = (x0, y0) .+ (t .* (x1 - x0, y1 - y0))
interpolate(x::Tuple, y::Tuple, t) = interpolate(x..., y..., t)

################################################################################
# Metadata for top level ports - used for pretty post-route plotting.
function top_level_port_metadata(style, orientation, direction, class, num_links)
    # Get a base dictionary for this port class.
    base = routing_metadata(class)
    # Parameters controlling offset generation.
    spacing = 0.4/num_links

    a, b = port_boundaries(style, orientation)
    start = initial_offset(style, orientation, direction)

    offset_dicts = map(0:num_links - 1) do i
        x, y = interpolate(a, b, start + spacing * i)
        d = Dict(
            "x" => x,
            "y" => y,
        )
        return merge(d, base)
    end

    return offset_dicts
end


text:
# EFTfitter.jl - Empty Template
# ============= Parameters =============================================#
parameters = BAT.NamedTupleDist(
    p1 = -2..2,
)

# ============= Observables ============================================#
function observable1(params)
    return params.p1
end

# ============= Measurements ===========================================#
measurements = (
    Meas1 = Measurement(observable1, 0.0, uncertainties = (unc1 = 0.1,), active=true),

    #MeasDist = MeasurementDistribution(obs_array, values_array, uncertainties = (unc1 = unc1_array,), active=false),
)

# ============= Correlations ===========================================#
correlations = (
    unc1 = NoCorrelation(active=true),
)

#corr_matrix = to_correlation_matrix(measurements,
#  (:Meas1, :Meas2, 0.1),
#)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using Distributions, Random
using Test, LinearAlgebra, PDMats

#  dimensions

n = 3
p = 6

#  parameters

M    = randn(n, p)
U    = rand(InverseWishart(n + 2, Matrix(1.0I, n, n)))
V    = rand(InverseWishart(p + 2, Matrix(1.0I, p, p)))
UF32 = Matrix{Float32}(U)
VBF  = Matrix{BigFloat}(V)
PDU  = PDMat(U)
PDV  = PDMat(V)
m    = randn(n)
w    = randn(p)
μ    = randn()
σ    = 2.0rand()
σ²   = σ ^ 2

#  distribution instances

D = MatrixNormal(M, U, V)  #  n x p
G = MatrixNormal(p, n)     #  p x n
L = MatrixNormal(reshape(m, n, 1), U, reshape([σ²], 1, 1))  #  n x 1
H = MatrixNormal(reshape(w, 1, p), reshape([σ²], 1, 1), V)  #  1 x p
K = MatrixNormal(reshape([μ], 1, 1), reshape([σ], 1, 1), reshape([σ], 1, 1))  #  1 x 1

d = vec(D) # MvNormal(vec(M), V ⊗ U)
g = MvNormal( Matrix(1.0I, p*n, p*n) )
l = MvNormal(m, σ² * U)
h = MvNormal(w, σ² * V)
k = Normal(μ, σ)

#  draws

A = rand(D)
B = rand(G)
X = rand(L)
Y = rand(H)
Z = rand(K)

a = vec(A)
b = vec(B)
x = vec(X)
y = vec(Y)
z = Z[1, 1]

@testset "Check all MatrixNormal constructors" begin
    @test MatrixNormal(M, U,        PDV.chol) isa MatrixNormal
    @test MatrixNormal(M, U,        PDV)      isa MatrixNormal
    @test MatrixNormal(M, PDU.chol, V)        isa MatrixNormal
    @test MatrixNormal(M, PDU.chol, PDV.chol) isa MatrixNormal
    @test MatrixNormal(M, PDU.chol, PDV)      isa MatrixNormal
    @test MatrixNormal(M, PDU,      V)        isa MatrixNormal
    @test MatrixNormal(M, PDU,      PDV.chol) isa MatrixNormal
    @test MatrixNormal(M, PDU,      PDV)      isa MatrixNormal
end

@testset "MatrixNormal promotion during construction" begin
    R = MatrixNormal(M, UF32, VBF)
    @test partype(R) == BigFloat
end

@testset "MatrixNormal construction errors" begin
    @test_throws ArgumentError MatrixNormal(M, V, U)
    @test_throws ArgumentError MatrixNormal(M, U, U)
end

@testset "MatrixNormal params" begin
    MM, PDUU, PDVV = params(D)
    @test MM == M
    @test U == PDUU.mat
    @test V == PDVV.mat
end

@testset "MatrixNormal size" begin
    @test size(D) == (n, p)
    @test size(G) == (p, n)
    @test size(L) == (n, 1)
    @test size(H) == (1, p)
    @test size(K) == (1, 1)
end

@testset "MatrixNormal rank" begin
    @test rank(D) == min(n, p)
    @test rank(G) == min(n, p)
    @test rank(L) == 1
    @test rank(H) == 1
    @test rank(K) == 1

    @test rank(D) == rank(rand(D))
    @test rank(G) == rank(rand(G))
    @test rank(L) == rank(rand(L))
    @test rank(H) == rank(rand(H))
    @test rank(K) == rank(rand(K))
end

@testset "MatrixNormal insupport" begin
    @test insupport(D, rand(D))
    @test insupport(G, rand(G))
    @test insupport(L, rand(L))
    @test insupport(H, rand(H))
    @test insupport(K, rand(K))

    @test !insupport(D, rand(G))
    @test !insupport(D, M + M * im)
    @test !insupport(L, rand(H))
    @test !insupport(H, rand(L))
end

@testset "MatrixNormal mean" begin
    @test mean(D) == M
    @test mean(G) == zeros(p, n)
    @test mean(L) == reshape(m, n, 1)
    @test mean(H) == reshape(w, 1, p)
    @test mean(K) == reshape([μ], 1, 1)
end

@testset "MatrixNormal mode" begin
    @test mode(D) == M
    @test mode(G) == zeros(p, n)
    @test mode(L) == reshape(m, n, 1)
    @test mode(H) == reshape(w, 1, p)
    @test mode(K) == reshape([μ], 1, 1)
end

@testset "MatrixNormal logpdf" begin
    @test logpdf(D, A) ≈ logpdf(d, a)
    @test logpdf(G, B) ≈ logpdf(g, b)
    @test logpdf(L, X) ≈ logpdf(l, x)
    @test logpdf(H, Y) ≈ logpdf(h, y)
    @test logpdf(K, Z) ≈ logpdf(k, z)
end

@testset "MatrixNormal sample moments" begin
    @test isapprox(mean(rand(D, 100000)), mean(D) , atol = 0.1)
    @test isapprox(cov(hcat(vec.(rand(D, 100000))...)'), kron(V, U) , atol = 0.1)
end

@testset "MatrixNormal conversion" for elty in (Float32, Float64, BigFloat)
    Del1 = convert(MatrixNormal{elty}, D)
    Del2 = convert(MatrixNormal{elty}, M, PDU, PDV, D.logc0)

    @test partype(Del1) == elty
    @test partype(Del2) == elty
end


text:
module Diffusion3D_Threads
    using ParallelStencil
    using ParallelStencil.FiniteDifferences3D
    @init_parallel_stencil(Threads, Float64, 3)
    include(joinpath(@__DIR__, "..", "..", "shared", "diffusion3D.jl"))
end


text:
# Julia wrapper for header: /usr/include/scip/cons_and.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrAnd(scip)
    ccall((:SCIPincludeConshdlrAnd, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsAnd(scip, cons, name, resvar, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsAnd, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Cint, Ptr{Ptr{SCIP_VAR}}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, resvar, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicAnd(scip, cons, name, resvar, nvars, vars)
    ccall((:SCIPcreateConsBasicAnd, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Cint, Ptr{Ptr{SCIP_VAR}}), scip, cons, name, resvar, nvars, vars)
end

function SCIPgetNVarsAnd(scip, cons)
    ccall((:SCIPgetNVarsAnd, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsAnd(scip, cons)
    ccall((:SCIPgetVarsAnd, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetResultantAnd(scip, cons)
    ccall((:SCIPgetResultantAnd, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPisAndConsSorted(scip, cons)
    ccall((:SCIPisAndConsSorted, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPsortAndCons(scip, cons)
    ccall((:SCIPsortAndCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPchgAndConsCheckFlagWhenUpgr(scip, cons, flag)
    ccall((:SCIPchgAndConsCheckFlagWhenUpgr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, flag)
end

function SCIPchgAndConsRemovableFlagWhenUpgr(scip, cons, flag)
    ccall((:SCIPchgAndConsRemovableFlagWhenUpgr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, flag)
end


text:
#------------------------------------------------------------------------------
"""
    excise(x...)

Remove all lines where the is a NaN/missing in any of the x arrays

# Examples
- `x1 = excise(x)`
- `(y1,x1) = excise(y,x)`

"""
function excise(x...)
  n  = length(x)
  vv = FindNNPs(x...)         #find rows with NaN/missing
  z = ntuple(i->copy(selectdim(x[i],1,vv)),n)    #create a tuple of arrays
  (n==1) && (z = z[1])                           #if a single array in the tuple
  return z
end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    FindNNPs(x...;Keepdim=1)

Find rows (if Keepdim==1) which have no NaNs missing in other dimensions (eg. in no columns).

# Input
- `z::Array`: one or several numerical arrays
- `Keepdim::Int`: (keyword) 1 if check rows, 2 if check columns, etc

# Output
- `vvb::BitVector`: vector, element t is true if row (if Keepdim==1) t has no NaN or missing

# Notice
- Set Keepdim=2 if we should instead look for NaNs/missings along rows (and other dimensions).
- For heterogenous arrays like `x=[x1,x1]`, use `FindNNPs(x...)`

Paul.Soderlind@unisg.ch

"""
function FindNNPs(x...;Keepdim=1)

  N  = length(x)
  T     = size(x[1],Keepdim)                 #length of output

  xDims = maximum(ndims,x)                   #max of ndims(x[i]) for i=1:N
  dims  = setdiff(1:xDims,Keepdim)           #dimensions to check

  vvM = falses(T,N)
  for i = 1:N                             #loop over inputs
    vvM[:,i] = any(z -> ismissing(z) || isnan(z),x[i],dims=dims)   #this order is important
  end

  vvb = vec(.!any(vvM,dims=2))      #rows witout NaN/missing in any of the x matrices

  return vvb

end
#------------------------------------------------------------------------------

max_stars_repo_name: duzaichuan/FinancialEconometrics

max_issues_repo_name: duzaichuan/FinancialEconometrics

max_forks_repo_name: duzaichuan/FinancialEconometrics


text:
using Pkg
using LinearAlgebra
using DataFrames
using Printf
using Random

println("test pf")

# load functions
include(pwd()*"/src/SDEMEM OU process/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU process/mcmc.jl")

N_time = 200
M_subjects = 40

nbr_particles_pf = 3000
ρ = 0.0
run_sort = false

y,x,t_vec,dt,η,σ_ϵ,ϕ,prior_parameters_η,prior_parameters_σ_ϵ = set_up(M=M_subjects,N=N_time,seed=100)

################################################################################
# Run pf at one # of particles
################################################################################

u_prop_old_block_1 = randn(nbr_particles_pf,N_time + 1,M_subjects)
u_resample_old_block_1 = randn(N_time,2,M_subjects)

u_prop_new = randn(nbr_particles_pf,N_time + 1,M_subjects)
u_resample_new = randn(N_time,2,M_subjects)

# correlat with old standard normal numbers
u_prop_tilde = ρ*u_prop_old_block_1 + sqrt(1-ρ^2)*u_prop_new
u_resample_tilde = ρ*u_resample_old_block_1 + sqrt(1-ρ^2)*u_resample_new

@time  cpf(y, σ_ϵ, ϕ, dt, u_prop_old_block_1, u_resample_old_block_1, nbr_particles_pf, run_sort)

function time_pf()

    for i = 1:100

        u_prop_old_block_1 = randn(nbr_particles_pf,N_time + 1,M_subjects)
        u_resample_old_block_1 = randn(N_time,2,M_subjects)

        u_prop_new = randn(nbr_particles_pf,N_time + 1,M_subjects)
        u_resample_new = randn(N_time,2,M_subjects)

        # correlat with old standard normal numbers
        u_prop_tilde = ρ*u_prop_old_block_1 + sqrt(1-ρ^2)*u_prop_new
        u_resample_tilde = ρ*u_resample_old_block_1 + sqrt(1-ρ^2)*u_resample_new

        ll =  cpf(y, σ_ϵ, ϕ, dt, u_prop_old_block_1, u_resample_old_block_1, nbr_particles_pf, run_sort)

    end

end

runtime_pf = @elapsed time_pf()

runtime_pf = runtime_pf/100

################################################################################
# Find correct nbr of particels
################################################################################

function loglik_est(nbr_particles_pf, nbr_pf_eval=100)

    ll = zeros(nbr_pf_eval)

    for i = 1:nbr_pf_eval
        u_prop_old_block_1 = randn(nbr_particles_pf,N_time + 1,M_subjects)
        u_resample_old_block_1 = randn(N_time,2,M_subjects)

        u_prop_new = randn(nbr_particles_pf,N_time + 1,M_subjects)
        u_resample_new = randn(N_time,2,M_subjects)

        # correlat with old standard normal numbers
        u_prop_tilde = ρ*u_prop_old_block_1 + sqrt(1-ρ^2)*u_prop_new
        u_resample_tilde = ρ*u_resample_old_block_1 + sqrt(1-ρ^2)*u_resample_new

        ll[i] = sum(cpf(y, σ_ϵ, ϕ, dt, u_prop_old_block_1, u_resample_old_block_1, nbr_particles_pf, run_sort))

    end

    return ll

end

nbr_pf_eval = 100
nbr_particles_pf = [100, 500, 1000, 1500, 2000, 2500, 3000, 4000, 5000]

loglik_matrix = zeros(nbr_pf_eval, length(nbr_particles_pf))
run_times = zeros(length(nbr_particles_pf))

run_times[1]  = @elapsed loglik_matrix[:,1] = loglik_est(nbr_particles_pf[1], nbr_pf_eval)

Random.seed!(1)

for i in 1:length(nbr_particles_pf)
    @printf "Current nbr of particels: %.2f\n" nbr_particles_pf[i]
    run_times[i]  = @elapsed loglik_matrix[:,i]  = loglik_est(nbr_particles_pf[i], nbr_pf_eval)
end

using PyPlot

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, mean(loglik_matrix, dims = 1)[:], "--*")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Avg. loglik est. val.")
PyPlot.savefig("figures/avg_loglik_val_vs_nbr_particels_pf.png")

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, var(loglik_matrix, dims = 1)[:], "--*")
PyPlot.plot(nbr_particles_pf, 2*ones(length(nbr_particles_pf),1), "k")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Var of loglik est.")
PyPlot.savefig("figures/loglik_var_vs_nbr_particles_pf.png")

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, run_times/100, "--*")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Run time (sec)")
PyPlot.savefig("figures/run_time_vs_nbr_particls_pf.png")


text:
using Dates: DateFormat, DatePart, min_width, max_width

function tryparsenext_fixedtz(str, i, len, min_width::Int=1, max_width::Int=0)
    i == len && str[i] == 'Z' && return ("Z", i+1)

    tz_start, tz_end = i, 0
    min_pos = min_width <= 0 ? i : i + min_width - 1
    max_pos = max_width <= 0 ? len : min(chr2ind(str, ind2chr(str,i) + max_width - 1), len)
    state = 1
    @inbounds while i <= max_pos
        c, ii = iterate(str, i)::Tuple{Char, Int}
        if state == 1 && (c == '-' || c == '+')
            state = 2
            tz_end = i
        elseif (state == 1 || state == 2) && '0' <= c <= '9'
            state = 3
            tz_end = i
        elseif state == 3 && c == ':'
            state = 4
            tz_end = i
        elseif (state == 3 || state == 4) && '0' <= c <= '9'
            tz_end = i
        else
            break
        end
        i = ii
    end

    if tz_end <= min_pos
        return nothing
    else
        tz = SubString(str, tz_start, tz_end)
        return tz, i
    end
end

function tryparsenext_tz(str, i, len, min_width::Int=1, max_width::Int=0)
    tz_start, tz_end = i, 0
    min_pos = min_width <= 0 ? i : i + min_width - 1
    max_pos = max_width <= 0 ? len : min(chr2ind(str, ind2chr(str,i) + max_width - 1), len)
    @inbounds while i <= max_pos
        c, ii = iterate(str, i)::Tuple{Char, Int}
        if c == '/' || c == '_' || isletter(c)
            tz_end = i
        else
            break
        end
        i = ii
    end

    if tz_end == 0
        return nothing
    else
        name = SubString(str, tz_start, tz_end)

        # If the time zone is recognized make sure that it is well-defined. For our
        # purposes we'll treat all abbreviations except for UTC and GMT as ambiguous.
        # e.g. "MST": "Mountain Standard Time" (UTC-7) or "Moscow Summer Time" (UTC+3:31).
        if occursin("/", name) || name in ("UTC", "GMT")
            return name, i
        else
            return nothing
        end
    end
end

function Dates.tryparsenext(d::DatePart{'z'}, str, i, len)
    tryparsenext_fixedtz(str, i, len, min_width(d), max_width(d))
end

function Dates.tryparsenext(d::DatePart{'Z'}, str, i, len)
    tryparsenext_tz(str, i, len, min_width(d), max_width(d))
end

function Dates.format(io::IO, d::DatePart{'z'}, zdt, locale)
    write(io, string(zdt.zone.offset))
end

function Dates.format(io::IO, d::DatePart{'Z'}, zdt, locale)
    write(io, string(zdt.zone))  # In most cases will be an abbreviation.
end

# Note: ISOZonedDateTimeFormat is defined in the module __init__ which means that this
# function can not be called from within this module. TODO: Ignore linting for this line
function ZonedDateTime(str::AbstractString, df::DateFormat=ISOZonedDateTimeFormat)
    try
        parse(ZonedDateTime, str, df)
    catch e
        if e isa ArgumentError
            rethrow(ArgumentError(
                "Unable to parse string \"$str\" using format $df. $(e.msg)"
            ))
        else
            rethrow()
        end
    end
end

function ZonedDateTime(str::AbstractString, format::AbstractString; locale::AbstractString="english")
    ZonedDateTime(str, DateFormat(format, locale))
end

Dates.default_format(::Type{ZonedDateTime}) = ISOZonedDateTimeFormat