text:
using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.BalanceLaws
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods.FVReconstructions: FVConstant, FVLinear
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Geometry
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.Orientations
using ClimateMachine.SystemSolvers
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using ClimateMachine.VTK

using CLIMAParameters
using CLIMAParameters.Planet: kappa_d
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates, Test

include("isentropicvortex_setup.jl")

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const output_vtk = false

function main()
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    numlevels = integration_testing ? 4 : 1

    expected_error = Dict()

    # Just to make it shorter and aligning
    Roe = RoeNumericalFlux

    # Float64, Dim 2, degree 4 in the horizontal, FV order 1, refinement level
    expected_error[Float64, FVConstant(), Roe, 1] = 3.5317756615538940e+01
    expected_error[Float64, FVConstant(), Roe, 2] = 2.5104217086071472e+01
    expected_error[Float64, FVConstant(), Roe, 3] = 1.6169569358521223e+01
    expected_error[Float64, FVConstant(), Roe, 4] = 9.4731749125284708e+00

    expected_error[Float64, FVLinear(), Roe, 1] = 2.6132907774912638e+01
    expected_error[Float64, FVLinear(), Roe, 2] = 8.8198392537283006e+00
    expected_error[Float64, FVLinear(), Roe, 3] = 2.4517109427575416e+00
    expected_error[Float64, FVLinear(), Roe, 4] = 7.4154384427579900e-01

    dims = 2
    @testset "$(@__FILE__)" begin
        for FT in (Float64,)
            for NumericalFlux in (Roe,)
                setup = IsentropicVortexSetup{FT}()
                for fvmethod in (FVConstant(), FVLinear())
                    @info @sprintf """Configuration
                                      FT                = %s
                                      ArrayType         = %s
                                      FV Reconstruction = %s
                                      NumericalFlux     = %s
                                      dims              = %d
                                      """ FT ArrayType fvmethod NumericalFlux dims
                    errors = Vector{FT}(undef, numlevels)

                    for level in 1:numlevels

                        # Match element numbers
                        numelems = (
                            2^(level - 1) * 5,
                            2^(level - 1) * 5 * polynomialorder,
                        )

                        errors[level] = test_run(
                            mpicomm,
                            ArrayType,
                            fvmethod,
                            polynomialorder,
                            numelems,
                            NumericalFlux,
                            setup,
                            FT,
                            dims,
                            level,
                        )

                        @test errors[level] ≈
                              expected_error[FT, fvmethod, NumericalFlux, level]
                    end
                    @info begin
                        msg = ""
                        for l in 1:(numlevels - 1)
                            rate = log2(errors[l]) - log2(errors[l + 1])
                            msg *= @sprintf(
                                "\n  rate for level %d = %e\n",
                                l,
                                rate
                            )
                        end
                        msg
                    end
                end

            end
        end
    end
end

function test_run(
    mpicomm,
    ArrayType,
    fvmethod,
    polynomialorder,
    numelems,
    NumericalFlux,
    setup,
    FT,
    dims,
    level,
)
    brickrange = ntuple(dims) do dim
        range(
            -setup.domain_halflength;
            length = numelems[dim] + 1,
            stop = setup.domain_halflength,
        )
    end
    connectivity = dims == 3 ? :full : :face

    topology = StackedBrickTopology(
        mpicomm,
        brickrange;
        periodicity = ntuple(_ -> true, dims),
        connectivity = connectivity,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = (polynomialorder, 0),
    )

    problem =
        AtmosProblem(boundaryconditions = (), init_state_prognostic = setup)

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = NoReferenceState(),
        turbulence = ConstantDynamicViscosity(FT(0)),
        moisture = DryModel(),
    )
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        problem = problem,
        orientation = NoOrientation(),
        source = (),
    )

    dgfvm = DGFVModel(
        model,
        grid,
        fvmethod,
        NumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    timeend = FT(2 * setup.domain_halflength / 10 / setup.translation_speed)

    # Determine the time step
    elementsize = minimum(step.(brickrange))
    dt = elementsize / soundspeed_air(param_set, setup.T∞) / polynomialorder^2
    nsteps = ceil(Int, timeend / dt)
    dt = timeend / nsteps

    Q = init_ode_state(dgfvm, FT(0))
    lsrk = LSRK54CarpenterKennedy(dgfvm, Q; dt = dt, t0 = 0)

    eng0 = norm(Q)
    dims == 2 && (numelems = (numelems..., 0))
    @info @sprintf """Starting refinement level %d
                      numelems  = (%d, %d, %d)
                      dt        = %.16e
                      norm(Q₀)  = %.16e
                      """ level numelems... dt eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            runtime = Dates.format(
                convert(DateTime, now() - starttime[]),
                dateformat"HH:MM:SS",
            )
            @info @sprintf """Update
                              simtime = %.16e
                              runtime = %s
                              norm(Q) = %.16e
                              """ gettime(lsrk) runtime energy
        end
    end
    callbacks = (cbinfo,)

    if output_vtk
        # Create vtk dir
        vtkdir =
            "vtk_isentropicvortex" *
            "_poly$(polynomialorder)_dims$(dims)_$(ArrayType)_$(FT)_level$(level)"
        mkpath(vtkdir)

        vtkstep = 0
        # Output initial step
        do_output(mpicomm, vtkdir, vtkstep, dgfvm, Q, Q, model)

        # Setup the output callback
        outputtime = timeend
        cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
            vtkstep += 1
            Qe = init_ode_state(dgfvm, gettime(lsrk), setup)
            do_output(mpicomm, vtkdir, vtkstep, dgfvm, Q, Qe, model)
        end
        callbacks = (callbacks..., cbvtk)
    end

    solve!(Q, lsrk; timeend = timeend, callbacks = callbacks)

    # Final statistics
    Qe = init_ode_state(dgfvm, timeend, setup)
    engf = norm(Q)
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished refinement level %d
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ level engf engf / eng0 engf - eng0 errf errf / engfe
    errf
end

function do_output(
    mpicomm,
    vtkdir,
    vtkstep,
    dgfvm,
    Q,
    Qe,
    model,
    testname = "isentropicvortex",
)
    ## Name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
    exactnames = statenames .* "_exact"

    writevtk(filename, Q, dgfvm, statenames, Qe, exactnames)

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## Name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(
            pvtuprefix,
            prefixes,
            (statenames..., exactnames...),
            eltype(Q),
        )

        @info "Done writing VTK: $pvtuprefix"
    end
end

main()


text:
using MarketData

facts("logic integrity") do

  context("things happen in order") do
      @fact 1 => roughly(1.0000000001)
  end

  context("transitions are defined") do
      @fact 1 => roughly(1.0000000001)
  end
end

ext: jl


text:
using Colors, GLVisualize
using DataFrames, Gadfly, GLVisualize.ComposeBackend

if !isdefined(:runtests)
    window = glscreen()
    composebackend = ComposeBackend.GLVisualizeBackend(window)
end
const static_example = true

xs = 0:0.1:20

df_cos = DataFrame(
    x=xs,
    y=cos(xs),
    ymin=cos(xs) .- 0.5,
    ymax=cos(xs) .+ 0.5,
    f="cos"
)
df_sin = DataFrame(
    x=xs,
    y=sin(xs),
    ymin=sin(xs) .- 0.5,
    ymax=sin(xs) .+ 0.5,
    f="sin"
)
df = vcat(df_cos, df_sin)
p = plot(df, x=:x, y=:y, ymin=:ymin, ymax=:ymax, color=:f, Geom.line, Geom.ribbon)

draw(composebackend, p)

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


text:
# ---
# title: 591. Tag Validator
# id: problem591
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Stack
# link: <https://leetcode.com/problems/tag-validator/description/>
# hidden: true
# ---
#
# Given a string representing a code snippet, you need to implement a tag
# validator to parse the code and return whether it is valid. A code snippet is
# valid if all the following rules hold:
#
#   1. The code must be wrapped in a **valid closed tag**. Otherwise, the code is invalid.
#   2. A **closed tag** (not necessarily valid) has exactly the following format : `<TAG_NAME>TAG_CONTENT</TAG_NAME>`. Among them, `<TAG_NAME>` is the start tag, and `</TAG_NAME>` is the end tag. The TAG_NAME in start and end tags should be the same. A closed tag is **valid** if and only if the TAG_NAME and TAG_CONTENT are valid.
#   3. A **valid** `TAG_NAME` only contain **upper-case letters** , and has length in range [1,9]. Otherwise, the `TAG_NAME` is **invalid**.
#   4. A **valid** `TAG_CONTENT` may contain other **valid closed tags** , **cdata** and any characters (see note1) **EXCEPT** unmatched `<`, unmatched start and end tag, and unmatched or closed tags with invalid TAG_NAME. Otherwise, the `TAG_CONTENT` is **invalid**.
#   5. A start tag is unmatched if no end tag exists with the same TAG_NAME, and vice versa. However, you also need to consider the issue of unbalanced when tags are nested.
#   6. A `<` is unmatched if you cannot find a subsequent `>`. And when you find a `<` or `</`, all the subsequent characters until the next `>` should be parsed as TAG_NAME (not necessarily valid).
#   7. The cdata has the following format : `<![CDATA[CDATA_CONTENT]]>`. The range of `CDATA_CONTENT` is defined as the characters between `<![CDATA[` and the **first subsequent** `]]>`.
#   8. `CDATA_CONTENT` may contain **any characters**. The function of cdata is to forbid the validator to parse `CDATA_CONTENT`, so even it has some characters that can be parsed as tag (no matter valid or invalid), you should treat it as **regular characters**.
#
# **Valid Code Examples:**
#
#
#
#     Input: "
This is the first line <![CDATA[
]]>
"
#
#     Output: True
#
#     Explanation:
#
#     The code is wrapped in a closed tag :
 and
.
#
#     The TAG_NAME is valid, the TAG_CONTENT consists of some characters and cdata.
#
#     Although CDATA_CONTENT has unmatched start tag with invalid TAG_NAME, it should be considered as plain text, not parsed as tag.
#
#     So TAG_CONTENT is valid, and then the code is valid. Thus return true.
#
#
#     Input: "
>>  ![cdata[]] <![CDATA[
]>]]>]]>>]
"
#
#     Output: True
#
#     Explanation:
#
#     We first separate the code into : start_tag|tag_content|end_tag.
#
#     start_tag -> **"
"**
#
#     end_tag -> **"
"**
#
#     tag_content could also be separated into : text1|cdata|text2.
#
#     text1 -> **">>  ![cdata[]] "**
#
#     cdata -> **" <![CDATA[
]>]]>"**, where the CDATA_CONTENT is **"
]>"**
#
#     text2 -> **"]]>>]"**
#
#
#     The reason why start_tag is NOT **"
>>"** is because of the rule 6.
#     The reason why cdata is NOT **" <![CDATA[
]>]]>]]>"** is because of the rule 7.
#
#
# **Invalid Code Examples:**
#
#
#
#     Input: "      "
#     Output: False
#     Explanation: Unbalanced. If "" is closed, then "" must be unmatched, and vice versa.
#
#     Input: "
  div tag is not closed
"
#     Output: False
#
#     Input: "
  unmatched <
"
#     Output: False
#
#     Input: "
 closed tags with invalid tag name  123
"
#     Output: False
#
#     Input: "
 unmatched tags with invalid tag name  </1234567890> and <CDATA[[]]>
"
#     Output: False
#
#     Input: "
  unmatched start tag   and unmatched end tag </C>
"
#     Output: False
#
#
# **Note:**
#
#   1. For simplicity, you could assume the input code (including the **any characters** mentioned above) only contain `letters`, `digits`, `'<'`,`'>'`,`'/'`,`'!'`,`'['`,`']'` and `' '`.
#
#
## @lc code=start
using LeetCode

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


text:
using Jecco
using Jecco.AdS5_3_1
using FFTW
using Plots
gr()

dirname   = "/home/mikel/Documentos/Jecco.jl/data/end_data/"
outdir    = "/home/mikel/Documentos/Jecco.jl/data/new_data/"

grid = SpecCartGrid3D(
    x_min            = -15.0,
    x_max            =  15.0,
    x_nodes          =  32,
    y_min            = -15.0,
    y_max            =  15.0,
    y_nodes          =  32,
    u_outer_min      =  0.1,
    u_outer_max      =  1.005,
    u_outer_domains  =  1,
    u_outer_nodes    =  64,
    u_inner_nodes    =  12,
    fd_order         =  4,
    sigma_diss       =  0.2,
)
io = InOut(recover_dir = dirname, out_dir = outdir, checkpoint_dir = outdir, out_boundary_every=1,out_gauge_every=1,out_bulk_every=1,remove_existing = true,)

parameters = AdS5_3_1.new_parameters(
#    e_new   = 1.5,
#    a4_ampy = 1.0,
#    a4_ky   = 2,
#    boostx = true,
    #fx20   = 1.0,
    #u_AH   = 0.9,
)

parameters_collision =AdS5_3_1.new_parameters_coll(
    dirname1  = dirname,
    dirname2  = dirname,
    x1_center = -15.0,
    y1_center = 0.0,
    x2_center = 15.0,
    y2_center = 0.0,
    fx21      = -0.1,
    fy21      = 0.0,
    fx22      = 0.1,
    fy22      = 0.0,
    u_AH      = 1.0,
)

AdS5_3_1.create_new_data(grid, io, parameters)
#AdS5_3_1.design_collision(grid, io, parameters_collision)

#=
phi11 = BulkTimeSeries(dirname,:phi,1)
phi12 = BulkTimeSeries(dirname,:phi,2)
phi21 = BulkTimeSeries(outdir,:phi,1)
phi22 = BulkTimeSeries(outdir,:phi,2)
=#

e  = VEVTimeSeries(outdir,:energy)
Jx = VEVTimeSeries(outdir,:Jx)
Jy = VEVTimeSeries(outdir,:Jy)

t,x,y = get_coords(e,:,:,:)
plan  = plan_rfft(e[1,:,:])
Nx    = length(x)
Ny    = length(y)
e0    = real(1/(Nx*Ny) * (plan * e[1,:,:]))[1]
Jx0   = real(1/(Nx*Ny) * (plan * Jx[1,:,:]))[1]
Jy0   = real(1/(Nx*Ny) * (plan * Jy[1,:,:]))[1]

#=
plan = plan_fft(e_old[end,:,:]);
Nx = length(x)
Ny = length(y)
dx = x[2]-x[1]
dy = y[2]-y[1]
ek = 1/(Nx*Ny)*(plan*e_old[end,:,:])
kx = 2*pi*fftfreq(Nx,1/dx)
ky = 2*pi*fftfreq(Ny,1/dy)

x0 = x[1]
y0 = y[1]

sum = 0.0
for j in 1:length(ky)
    for i in 1:length(kx)
        global sum += ek[i,j]*exp(im*(kx[i]*x0+ky[j]*y0))
    end
end
=#
println("New set up is:")
println("(xmin, xmax) = ($(x[1]),$(x[end]+x[2]-x[1]))")
println("(ymin, ymax) = ($(y[1]),$(y[end]+y[2]-y[1]))")
println("Nx, Ny = $(length(x)), $(length(y))")
println("Average Energy Density = $e0")
println("Average x momenta = $Jx0")
println("Average y momenta = $Jy0")
println("Maximum x momenta = $(maximum(abs.(Jx[end,:,:])))")
println("Maximum y momenta = $(maximum(abs.(Jy[end,:,:])))")
#=
plot(x,y,e[1,:,:],st=:surface, camera=(50,65))
xlabel!("x")
ylabel!("y")
=#


text:
module MATLABPlots

import MATLAB

#Basic One Liners
#TODO: SAVE
figure(n::Int) = MATLAB.eval_string("figure($n)")
figure() = MATLAB.eval_string("figure()")
function hold_on()
    MATLAB.mat"hold on"
end
function hold_off()
    MATLAB.mat"hold off"
end
clf(n::Int) = MATLAB.eval_string("clf($n)")
clf() = MATLAB.eval_string("clf()")
subplot(a,b,c) = MATLAB.eval_string("subplot($a,$b,$c)")
export figure,clf,hold_on,hold_off,subplot

#Utilities for Plotting
include("utils.jl")

#(1D and 2D) Plotting
include("2DMatlabPlots.jl")
export plot,semilogx,loglog,semilogy,stairs

#Timeseries
using Dates
include("timeseries.jl")

#(3D) Plotting
include("3DMatlabPlots.jl")
export plot3,surf,contourf,contour,surface

#(Statistical) Plotting
include("StatisticalMatlabPlotting.jl")
export histogram

#Geo Plotting
include("GeoPlots.jl")
export geoplot

#Legends, Axis Labels, GCA
include("Style.jl")
export xlabel,ylabel,zlabel,title,legend,setgca

# Saving Plots
function saveas(filename)
	MATLAB.eval_string("saveas(gcf,\"$filename\")")
end
export saveas
end


text:
function solve!(solver::ProjectedNewtonSolver)
    reset!(solver)

    update_constraints!(solver)
    copy_constraints!(solver)
    constraint_jacobian!(solver)
    copy_jacobians!(solver)
    cost_expansion!(solver)
    update_active_set!(solver)
    copy_active_set!(solver)

    if solver.opts.verbose
        println("\nProjection:")
    end
    viol = projection_solve!(solver)
    copyto!(solver.Z, solver.P)
end

function projection_solve!(solver::ProjectedNewtonSolver)
    ϵ_feas = solver.opts.constraint_tolerance
    viol = norm(solver.d[solver.active_set], Inf)
    max_projection_iters = solver.opts.n_steps

    count = 0
    while count < max_projection_iters && viol > ϵ_feas
        viol = _projection_solve!(solver)
        count += 1
        record_iteration!(solver, viol)
    end
    return viol
end

function record_iteration!(solver::ProjectedNewtonSolver, viol)
    solver.stats.iterations += 1
    i = solver.stats.iterations
    solver.stats.cost[i] = cost(solver)
    solver.stats.c_max[i] = viol
end

function reset!(solver::ProjectedNewtonSolver)
    solver.stats.iterations = 0
    solver.stats.cost .*= 0
    solver.stats.c_max .*= 0
    return nothing
end

function _projection_solve!(solver::ProjectedNewtonSolver)
    Z = primals(solver)
    a = solver.active_set
    max_refinements = 10
    convergence_rate_threshold = solver.opts.r_threshold
    ρ = solver.opts.ρ

    # Assume constant, diagonal cost Hessian (for now)
    H = Diagonal(solver.H)

    # Update everything
    update_constraints!(solver)
    constraint_jacobian!(solver)
    update_active_set!(solver)
    cost_expansion!(solver)

    # Copy results from constraint sets to sparse arrays
    copyto!(solver.P, solver.Z)
    copy_constraints!(solver)
    copy_jacobians!(solver)
    copy_active_set!(solver)

    # Get active constraints
    D,d = active_constraints(solver)

    viol0 = norm(d,Inf)
    if solver.opts.verbose
        println("feas0: $viol0")
    end

    HinvD = H\D'
    S = Symmetric(D*HinvD)
    Sreg = cholesky(S + ρ*I)
    viol_prev = viol0
    count = 0
    while count < max_refinements
        viol = _projection_linesearch!(solver, (S,Sreg), HinvD)
        convergence_rate = log10(viol) / log10(viol_prev)
        viol_prev = viol
        count += 1

        if solver.opts.verbose
            println("conv rate: $convergence_rate")
        end

        if convergence_rate < convergence_rate_threshold ||
                       viol < solver.opts.constraint_tolerance
            break
        end
    end
    return viol_prev
end

function _projection_linesearch!(solver::ProjectedNewtonSolver,
        S, HinvD)
    conSet = get_constraints(solver)
    a = solver.active_set
    d = solver.d[a]
    viol0 = norm(d,Inf)
    viol = Inf

    P = solver.P
    Z = solver.Z
    P̄ = solver.P̄
    Z̄ = solver.Z̄

    solve_tol = 1e-8
    refinement_iters = 25
    α = 1.0
    ϕ = 0.5
    count = 1
    while true
        δλ = reg_solve(S[1], d, S[2], solve_tol, refinement_iters)
        δZ = -HinvD*δλ
        P̄.Z .= P.Z + α*δZ

        copyto!(Z̄, P̄)
        update_constraints!(solver, Z̄)
        max_violation!(conSet)
        viol_ = maximum(conSet.c_max)
        copy_constraints!(solver)
        d = solver.d[a]
        viol = norm(d,Inf)

        if solver.opts.verbose
            println("feas: ", viol, " (", viol_, ")")
        end
        if viol < viol0 || count > 10
            break
        else
            count += 1
            α *= ϕ
        end
    end
    copyto!(P.Z, P̄.Z)
    # copyto!(solver.Z, P.Z)
    return viol
end

reg_solve(A, b, reg::Real, tol=1e-10, max_iters=10) = reg_solve(A, b, A + reg*I, tol, max_iters)
function reg_solve(A, b, B, tol=1e-10, max_iters=10)
    x = B\b
    count = 0
    while count < max_iters
        r = b - A*x
        # println("r_norm = $(norm(r))")

        if norm(r) < tol
            break
        else
            x += B\r
            count += 1
        end
    end
    # println("iters = $count")

    return x
end

function active_constraints(solver::ProjectedNewtonSolver)
    return solver.D[solver.active_set, :], solver.d[solver.active_set]  # this allocates
end

function cost_expansion!(solver::ProjectedNewtonSolver)
    Z = get_trajectory(solver)
    E = solver.E
    obj = get_objective(solver)
    cost_expansion!(E, obj, Z)

    xinds, uinds = solver.P.xinds, solver.P.uinds
    H = solver.H
    g = solver.g
    copy_expansion!(H, g, E, xinds, uinds)
    return nothing
end

function copy_expansion!(H, g, E, xinds, uinds)
    N = length(E)

    for k = 1:N-1
        H[xinds[k],xinds[k]] .= E[k].xx
        H[uinds[k],uinds[k]] .= E[k].uu
        H[uinds[k],xinds[k]] .= E[k].ux
        g[xinds[k]] .= E[k].x
        g[uinds[k]] .= E[k].u
    end
    H[xinds[N],xinds[N]] .= E[N].xx
    g[xinds[N]] .= E[N].x
    return nothing
end

@inline copy_constraints!(solver::ProjectedNewtonSolver) = copy_constraints!(solver.d, solver)
@inline copy_jacobians!(solver::ProjectedNewtonSolver) = copy_jacobians!(solver.D, solver)
@inline copy_active_set!(solver::ProjectedNewtonSolver) = copy_active_set!(solver.active_set, solver)


text:
"""A `Bank1D` is a one-dimensional wavelet filter bank, parametrized by
* `T`: numeric type of input, e.g. Float32, Float64.
* `D`: transform domain. Either `FourierDomain{1}` or `SpatialDomain{1}`.
* `G`: point group. Either `TrivialGroup` or `ReflectionGroup`.
* `W`: wavelet class, e.g. `Morlet` or `Gammatone`.
Its fields are
* `ϕ`: low-pass filter, also called scaling function.
* `ψs`: 3d array of wavelets, indexed by spin, chroma, and octave.
* `behavior`: mutable behavior, e.g. oversampling.
* `spec`: immutable specifications, e.g. number of filters per octave.
To create a `Bank1D`
1. define a `Spec1D`,
2. define a `PathKey` on which to apply the wavelet filterbank,
3. if needed, provide behavior characteristics as keyword arguments.
Example:
spec = Spec1D(n_filters_per_octave = 12, n_octaves = 10)
bank = Bank1D(spec, j_range = 3:10)"""
immutable Bank1D{
        T<:Number,
        D<:LineDomains,
        G<:LineGroups,
        W<:RedundantWaveletClass} <: AbstractBank{T,D,G,W}
    ϕ::AbstractFilter{T,D}
    ψs::Array{AbstractFilter{T,D},3}
    behavior::Behavior{T}
    spec::Spec1D{T,D,G,W}
    function (::Type{Bank1D}){T,D,G,W}(
            spec::Spec1D{T,D,G,W},
            pathkey::PathKey = PathKey(:time) ;
            is_ϕ_applied::Bool = false,
            j_range::UnitRange{Int} = 0:(spec.n_octaves-1),
            log2_oversampling::Int = 0,
            max_log2_stride::Int = spec.n_octaves-1,
            weighting::AbstractWeighting = EqualWeighting())
        (nΘs, nΧs, nJs) = size(spec.ψmetas)
        ψs = Array(AbstractFilter{T,D}, (nΘs, nΧs, nJs))
        ψs[1, :, :] =
            pmap(AbstractFilter, spec.ψmetas[1, :, :], fill(spec, nΧs, nJs))
        (nΘs > 1) && spin!(ψs)
        ϕ = AbstractFilter(spec.ϕmeta, spec)
        renormalize!(ϕ, ψs, spec)
        behavior = Behavior(ϕ, ψs, spec, is_ϕ_applied, j_range,
            log2_oversampling, max_log2_stride, pathkey, weighting)
        new{T,D,G,W}(ϕ, ψs, behavior, spec)
    end
end

"""Collects band-pass filters `ψ`s and low-pass filter `ϕ` in the Fourier domain. The band-pass filters are organized as a 4-D tensor, where the
dimensions correspond to frequency `ω`, orientation `θ`, chroma `χ`,
and octave index `j`."""
function Base.collect{T}(bank::Bank1D{T,FourierDomain{1}})
    N = 1 << bank.spec.log2_size
    halfN = N >> 1
    (nΘs, nΧs, nJs) = size(bank.spec.ψmetas)
    ψs = zeros(T, N, nΘs, nΧs, nJs)
    ωs = (-halfN):(halfN-1)
    for j in 0:(nJs-1), χ in 0:(nΧs-1), θ in 0:(nΘs-1)
        ψs[:, 1+θ, 1+χ, 1+j] =
            fftshift(bank.ψs[1+θ, 1+χ, 1+j][ωs])
    end
    ϕ = fftshift(bank.ϕ[ωs])
    return (ψs, ϕ)
end

Base.ndims(::Bank1D) = 1


text:
using DPC, CairoMakie
include("utils.jl")
# set_theme!(presentation_theme)
set_theme!(mytheme)

config0 = """
neurons:
  - id: C
    θ_syn: 0
    branches:
      - id: A
        θ_syn: 5
"""

config1 = """
neurons:
  - id: C
    θ_syn: 0
    θ_seg: 1
    branches:
      - id: A
        θ_syn: 5
        branches:
          - id: B
            θ_syn: 5
"""

config2 = """
neurons:
  - id: C
    θ_syn: 0
    θ_seg: 1
    branches:
      - id: A
        θ_syn: 5
      - id: B
        θ_syn: 5
"""

num_synapses = 10
prob_A = fill(0.0, num_synapses)
prob_B = fill(0.0, num_synapses)
volleys = Dict(:A => (prob_A, ones(num_synapses)), :B => (prob_B, ones(num_synapses)))

probabilities = 0:0.01:1

firing_probabilities0 = Vector{Vector{Float64}}(undef, num_synapses)
(net0,objects0) = load_network(YAML_source=config0)
for θ in 1:num_synapses
    println("Computing for single segment, θ=$(θ)")
    objects0[:A].θ_syn = θ
    firing_probabilities0[θ] = [(prob_A .= a; P_fire(objects0[:C], volleys)) for a ∈ probabilities]
end

println("Computing for AND neuron")
(net1,objects1) = load_network(YAML_source=config1)
firing_probabilities1 = [(prob_A .= a; prob_B .= b; P_fire(objects1[:C], volleys)) for a ∈ probabilities, b ∈ probabilities]

println("Computing for OR neuron")
(net2,objects2) = load_network(YAML_source=config2)
firing_probabilities2 = [(prob_A .= a; prob_B .= b; P_fire(objects2[:C], volleys)) for a ∈ probabilities, b ∈ probabilities]

## Plotting

fig = Figure(resolution=(0.75textwidth,0.7textwidth))
ax1a = fig[1,1] = Axis(fig, backgroundcolor=:transparent,
    title="a.    single\n       segment", titlealign=:left, tellwidth=false
)

ax2a = fig[1,2] = Axis(fig, backgroundcolor=:transparent,
    title="b.    sequential\n       segments", titlealign=:left, tellwidth=false
)

ax3a = fig[1,3] = Axis(fig, backgroundcolor=:transparent,
    title="c.    parallel\n       segments", titlealign=:left, tellwidth=false
)

ax1b = fig[2,1] = Axis(fig, xlabel="p₁", ylabel="Plateau probability")#, aspect=DataAspect()
ax2b = fig[2,2] = Axis(fig, xlabel="p₁", ylabel="p₂", )
ax3b = fig[2,3] = Axis(fig, xlabel="p₁", ylabel="p₂", )
cell_col_1 = fig[3,1]
cell_col_2 = fig[3,2:3]

colors = cgrad(:viridis,length(firing_probabilities0), categorical=true)
for (prob,col) in zip(firing_probabilities0, colors)
    lines!(ax1b, probabilities, prob, color=col, linewidth=3)
end
ax_col_1 = Colorbar(cell_col_1, height=10, colormap=colors, limits=(0.5,10.5), ticks=1:10, label="TS₁", vertical=false, flipaxis = false)

contourf!(ax2b, probabilities, probabilities, firing_probabilities1, colormap=:viridis, colorrange=(0,1))

contourf!(ax3b, probabilities, probabilities, firing_probabilities2, colormap=:viridis, colorrange=(0,1))

p_n=plot!(ax3a, objects2[:C], ports=Dict(:C=>[],:A=>[:A, :dummy],:B=>[:dummy, :B]), angle_between=20/180*π, branch_width=0.2, branch_length=1.0, color=Dict(:C=>RGBAf0(0.5,0.5,0.5,0.5), :A=>color_1, :B=>color_2))
ports = Dict(p_n.attributes[:ports][])
arrows!(ax3a, [ports[:A]-Point2f0(0.6,0),ports[:B]-Point2f0(0.45,0)], [Point2f0(0.55,0),Point2f0(0.4,0)], linewidth=2, arrowsize = [10,10])
text!.(ax3a, "A~B(10,P₁)", position=ports[:A]-Point2f0(0.6,0), align=(:left, :top), textsize=14, color=:black)
text!.(ax3a, "B~B(10,P₂)", position=ports[:B]-Point2f0(0.45,0), align=(:left, :bottom), textsize=14, color=:black)

p_n=plot!(ax2a, objects1[:C], ports=Dict(:C=>[],:A=>[:A],:B=>[:B]), angle_between=20/180*π, branch_width=0.2, branch_length=1.0, color=Dict(:C=>RGBAf0(0.5,0.5,0.5,0.5), :A=>color_1, :B=>color_2))
ports = Dict(p_n.attributes[:ports][])
arrows!(ax2a, [ports[:A]-Point2f0(0.35,0),ports[:B]-Point2f0(0.35,0)] , Node([Point2f0(0.3,0),Point2f0(0.3,0)]), linewidth=2, arrowsize = [10,10])
text!.(ax2a, "A~B(10,P₁)", position=ports[:A]-Point2f0(0.35,0), align=(:left, :bottom), textsize=14, color=:black)
text!.(ax2a, "B~B(10,P₂)", position=ports[:B]-Point2f0(0.35,0), align=(:left, :bottom), textsize=14, color=:black)

p_n=plot!(ax1a, objects0[:C], ports=Dict(:C=>[:C],:A=>[:A]), angle_between=20/180*π, branch_width=0.2, branch_length=1.0, color=Dict(:C=>RGBAf0(0.5,0.5,0.5,0.5), :A=>color_1, :B=>color_2))
ports = Dict(p_n.attributes[:ports][])
arrows!(ax1a, [ports[:A]-Point2f0(0.35,0)], [Point2f0(0.3,0)], linewidth=2, arrowsize = 10)
text!.(ax1a, "A~B(10,P₁)", position=ports[:A]-Point2f0(0.35,0), align=(:left, :bottom), textsize=14, color=:black)

ax_col_2 = Colorbar(cell_col_2, height=10, colormap=:viridis, limits=(0,1), label="Plateau probability for TS₁=TS₂=5", vertical=false, flipaxis = false)

hidedecorations!(ax3a)
hidespines!(ax3a)
hidedecorations!(ax2a)
hidespines!(ax2a)
hidedecorations!(ax1a)
hidespines!(ax1a)

hideydecorations!(ax3b)
colgap!(fig.layout, 1, Fixed(45))
colgap!(fig.layout, 2, Fixed(45))
rowsize!(fig.layout, 1, Aspect(1, 1.5))
rowsize!(fig.layout, 2, Aspect(1,1))

ylims!(ax1a, -0.25,2.25)
ylims!(ax2a, -0.25,2.25)
ylims!(ax3a, -0.25,2.25)
xlims!(ax1a, -0.8333,0.8333)
xlims!(ax2a, -0.8333,0.8333)
xlims!(ax3a, -0.8333,0.8333)

xlims!(ax1b, 0,1)
xlims!(ax2b, 0,1)
xlims!(ax3b, 0,1)
ylims!(ax1b, 0,1)
ylims!(ax2b, 0,1)
ylims!(ax3b, 0,1)
ax2b.yticks[] = [0,0.5,1.0]

# save(joinpath("figures","probabilistic.pdf"), fig)
save(joinpath("figures","probabilistic.svg"), fig)
# save(joinpath("figures","probabilistic.png"), fig)
fig


text:
#-------------------------------
#Allows reading in N
#-------------------------------
include("QuadraturePoints.jl")

#Some Constants
DFloat = Float16
N=parse(Int,ARGS[1])

#Some Output
@show (N)

# {{{ main
function main()

#Chebyshev Points
(ξ,ω) = QuadraturePoints.chebyshev_gauss(N+1)
println("Chebyshev Points ") #output
@show (ξ,ω)

#Legendre Points
(ξ,ω) = QuadraturePoints.legendre_gauss(N+1)
println("Legendre Points ") #output
@show (ξ,ω)

#Lobatto Points
(ξ,ω) = QuadraturePoints.lobatto_gauss(N+1)
println("Lobatto Points ") #output
@show (ξ,ω)

end

main()


text:
module RK3TVD
export Rk3tvd!,RK3TVDData
# RK explicit method, TVD.
struct RK3TVDData
    size::Int64
    c21::Float64
    c22::Float64
    c31::Float64
    c32::Float64
    u1::Array{Float64,1}
    u2::Array{Float64,1}
RK3TVDData(Size)=new(Size,3./4.,1./4,1./3.,2./3.,zeros(Size),zeros(Size))
end
function Rk3tvd!(R::RK3TVDData,W,dt,In::Array{Float64},Out::Array{Float64})
    size=R.size

    W(In,R.u1)
    #R.u1=In + dt*R.u1
    @simd for i=1:size
        R.u1[i]=In[i] + dt*R.u1[i]
    end

    W(R.u1,R.u2)
    #R.u2= R.c21*In+R.c22*(R.u1+dt*R.u2)
    @simd for i=1:size
        R.u2[i]= R.c21*In[i]+R.c22*(R.u1[i]+dt*R.u2[i])
    end

    W(R.u2,R.u1)
    @simd for i=1:size
        Out[i]=R.c31*In[i]+R.c32*(R.u2[i]+dt*R.u1[i])
    end
    #Out= @. R.c31*In+R.c32*(R.u2+dt*R.u1)
    nothing
end
end

max_stars_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_issues_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_forks_repo_name: wolfv/BenchmarksPythonJuliaAndCo


text:
uciisready() = println("readyok")

ucinewgame() = engineGame{SimpleGame}(SimpleGame())

function position!(eg::engineGame, commands)
    g = SimpleGame()
    endindex = 1
    if commands[1] == "fen"
        deleteat!(commands, 1)
        endindex = findlast(s -> isdigit(s[1]), commands)
        s = join(commands[1:endindex], ' ')
        println(s)
        g = SimpleGame(fromfen(s))
        #deleteat!(commands, 1)
    end
    for i = 1:endindex
        deleteat!(commands, 1)
    end
    eg.g = g
    g
end

function moves!(eg::engineGame, commands)
    lastmove = Move(5)
    g = eg.g
    while !isempty(commands) && 'a' <= commands[1][1] <= 'h'
        lastmove = movefromstring(String(commands[1]))
        if isnothing(lastmove)
            break
        end
        domove!(g, lastmove)
        deleteat!(commands, 1)
    end
end

function go!(eg::engineGame, commands; stop = Ref(false))
    availableTime = Millisecond(0)
    depthSet = false
    @debug(commands)
    while !isempty(commands)
        if (commands[1] == "wtime")
            deleteat!(commands, 1)
            eg.wtime = Millisecond(parse(Int, commands[1]))
            deleteat!(commands, 1)
        elseif (commands[1] == "btime")
            deleteat!(commands, 1)
            eg.btime = Millisecond(parse(Int, commands[1]))
            deleteat!(commands, 1)
        elseif (commands[1] == "winc")
            deleteat!(commands, 1)
            eg.winc = Millisecond(parse(Int, commands[1]))
            deleteat!(commands, 1)
        elseif (commands[1] == "binc")
            deleteat!(commands, 1)
            eg.binc = Millisecond(parse(Int, commands[1]))
            deleteat!(commands, 1)
        elseif (commands[1] == "movestogo")
            deleteat!(commands, 1)
            #eg.movestogo = parse(Int, commands[1])
            deleteat!(commands, 1)
        elseif (commands[1] == "depth")
            deleteat!(commands, 1)
            eg.depth = parse(Int, commands[1])
            deleteat!(commands, 1)
            depthSet = true
        elseif (commands[1] == "movetime")
            deleteat!(commands, 1)
            availableTime = Millisecond(parse(Int, commands[1]))
            deleteat!(commands, 1)
        end
    end
    buffer = false
    if availableTime == Millisecond(0)
        buffer = true
        if sidetomove(board(eg.g)) == WHITE
            availableTime += min((eg.wtime * 7) ÷ 10, eg.winc + eg.wtime ÷ 60)
        else
            availableTime += min((eg.btime * 7) ÷ 10, eg.binc + eg.btime ÷ 60)
        end
        eg.depth = max(1, eg.depth - 1)
    end
    iDepth = Ref(eg.depth)
    if depthSet
        move = @async negamax(
            eg.g,
            eg.depth,
            1,
            true;
            stop = stop,
            uci = true,
            evalFunc = pstFull,
        )
    else
        move = @async iterativeDeepening(
            eg.g,
            availableTime,
            iDepth;
            buffer = buffer,
            stop = stop,
            uci = true,
            evalFunc = pstAndMate,
        )
    end
    return move
    println("bestmove $(tostring(move))")
end

function stop(eg, st = Channel{Bool}(5))
    #put!(st, true)
    st[] = true
    yield()
    #move = rand(moves(board(eg.g)))
end

function uci_listen(; id = true)
    if id
        println("id name $myId")
        println("id author David Weingut")
        println("uciok")
    end
    while true
        ln = readline()
        ln == "isready" && break
        ln == "quit" && return
    end
    @time iterativeDeepening(
        Game(),
        Millisecond(100),
        Ref(1);
        evalFunc = pstFull,
    )
    eg = engineGame{SimpleGame}(SimpleGame())
    st = Ref{Bool}(false) #Channel{Bool}(5)
    fenArr = ["fen", split("rnb1kbnr/pp1p1ppp/4p3/1Bp3q1/8/2P1P3/PP1P1PPP/RNBQK1NR w KQkq - 0 1")...]
    @time position!(eg, ["fen", split("rnb1kbnr/pp1p1ppp/4p3/1Bp3q1/8/2P1P3/PP1P1PPP/RNBQK1NR w KQkq - 0 1")...])
    println(board(eg.g))
    println(moves(board(eg.g)))
    moves!(eg, ["a2a3"])

    uciisready()
    while true
        @debug(st[])#=isready(st)=#
        #uciisready()
        ln = readline()
        @debug("read line")
        commands = split(ln)
        #println("info string $commands")
        startLength = length(commands) + 1
        while !isempty(commands) && length(commands) < startLength
            startLength = length(commands)
            com = Ref(commands, 1)
            if com[] == "quit"
                return
            end
            if com[] == "isready"
                uciisready()
                break
            end
            if com[] == "ucinewgame"
                eg = ucinewgame()
                #uciisready()
                break
            end
            if com[] == "position"
                deleteat!(commands, 1)
                eg.g = position!(eg, commands)
                #println("position returned")
                if isempty(commands)
                    @debug "commands emptied"
                    #uciisready()
                    break
                end
            end
            if com[] == "moves"
                deleteat!(commands, 1)
                moves!(eg, commands)
                if isempty(commands)
                    #uciisready()
                    break
                end
            end
            if com[] == "go"
                @debug("eg.depth = $(eg.depth)")
                deleteat!(commands, 1)
                st[] = false
                go!(eg, commands; stop = st)
                @debug("go rturned")
                isempty(commands) && break
            end
            if com[] == "stop"
                stop(eg, st)
                @debug("returned from stop()")
            end
        end
    end
end


text:
module ExpiringCaches

using Dates

export Cache, @cacheable

struct TimestampedValue{T}
    value::T
    timestamp::DateTime
end

TimestampedValue(x::T) where {T} = TimestampedValue{T}(x, Dates.now(Dates.UTC))
TimestampedValue{T}(x) where {T} = TimestampedValue{T}(x, Dates.now(Dates.UTC))

"""
    ExpiringCaches.Cache{K, V}(timeout::Dates.Period)

Create a thread-safe, expiring cache where values older than `timeout`
are "invalid" and will be deleted.

An `ExpiringCaches.Cache` is an `AbstractDict` and tries to emulate a regular
`Dict` in all respects. It is most useful when the cost of retrieving or
calculating a value is expensive and is able to be "cached" for a certain
amount of time. To avoid using the cache (i.e. to invalidate the cache),
a `Cache` supports the `delete!` and `empty!` methods to remove values
manually.
"""
struct Cache{K, V, P <: Dates.Period} <: AbstractDict{K, V}
    lock::ReentrantLock
    cache::Dict{K, TimestampedValue{V}}
    timeout::P
end
Cache{K, V}(timeout::Dates.Period=Dates.Minute(1)) where {K, V} = Cache(ReentrantLock(), Dict{K, TimestampedValue{V}}(), timeout)

expired(x::TimestampedValue, timeout) = (Dates.now(Dates.UTC) - x.timestamp) > timeout

function Base.iterate(x::Cache)
    lock(x.lock)
    state = iterate(x.cache)
    if state === nothing
        unlock(x.lock)
        return nothing
    end
    while expired(state[1][2], x.timeout)
        state = iterate(x.cache, state[2])
        if state === nothing
            unlock(x.lock)
            return nothing
        end
    end
    return (state[1][1], state[1][2].value), state[2]
end

function Base.iterate(x::Cache, st)
    state = iterate(x.cache, st)
    if state === nothing
        unlock(x.lock)
        return nothing
    end
    while expired(state[1][2], x.timeout)
        state = iterate(x.cache, state[2])
        if state === nothing
            unlock(x.lock)
            return nothing
        end
    end
    return (state[1][1], state[1][2].value), state[2]
end

function Base.haskey(cache::Cache{K, V}, k::K) where {K, V}
    lock(cache.lock) do
        if haskey(cache.cache, k)
            x = cache.cache[k]
            if !expired(x, cache.timeout)
                return true
            else
                delete!(cache.cache, k)
                return false
            end
        end
        return false
    end
end

function Base.setindex!(cache::Cache{K, V}, val::V, key::K) where {K, V}
    lock(cache.lock) do
        cache.cache[key] = TimestampedValue{V}(val)
        return val
    end
end

function Base.get(cache::Cache{K, V}, key::K, default::V) where {K, V}
    lock(cache.lock) do
        if haskey(cache.cache, key)
            x = cache.cache[key]
            if expired(x, cache.timeout)
                delete!(cache.cache, key)
                return default
            else
                return x.value
            end
        else
            return default
        end
    end
end

function Base.get!(cache::Cache{K, V}, key::K, default::V) where {K, V}
    lock(cache.lock) do
        if haskey(cache.cache, key)
            x = cache.cache[key]
            if expired(x, cache.timeout)
                return setindex!(cache, default, key)
            else
                return x.value
            end
        else
            return setindex!(cache, default, key)
        end
    end
end

function Base.get!(f::Function, cache::Cache{K, V}, key::K) where {K, V}
    lock(cache.lock) do
        if haskey(cache.cache, key)
            x = cache.cache[key]
            if expired(x, cache.timeout)
                return setindex!(cache, f()::V, key)
            else
                return x.value
            end
        else
            return setindex!(cache, f()::V, key)
        end
    end
end

Base.delete!(cache::Cache{K}, key::K) where {K} = lock(() -> delete!(cache.cache, key), cache.lock)
Base.empty!(cache::Cache) = lock(() -> empty!(cache.cache), cache.lock)
Base.length(cache::Cache) = length(cache.cache)

"""
    @cacheable timeout function_definition::ReturnType

For a function definition (`function_definition`, either short-form
or full), create an `ExpiringCaches.Cache` and store results for `timeout`
(hashed by the exact input arguments obviously).

Note that the function definition _MUST_ include the `ReturnType` declartion
as this is used as the value (`V`) type in the `Cache`.
"""
macro cacheable(timeout, func)
    @assert func.head == :function
    func.args[1].head == :(::) || throw(ArgumentError("@cacheable function must specify return type: $func"))
    returnType = func.args[1].args[2]
    sig = func.args[1].args[1]
    functionBody = func.args[2]
    funcName = sig.args[1]
    internalFuncName = Symbol("__$funcName")
    sig.args[1] = internalFuncName
    funcArgs = sig.args[2:end]
    argTypes = map(x->x.args[2], funcArgs)
    internalFunction = Expr(:function, sig, functionBody)
    cacheName = gensym()
    return esc(quote
        const $cacheName = ExpiringCaches.Cache{Tuple{$(argTypes...)}, $returnType}($timeout)
        $internalFunction
        function $funcName(args...)::$returnType
            return get!($cacheName, args) do
                $internalFuncName(args...)
            end
        end
        ExpiringCaches.getcache(f::typeof($funcName)) = $cacheName
        $funcName
    end)
end

function getcache end

# @cacheable Dates.Minute(2) function foo(arg1::Int, arg2::String)::ReturnType
#     x = x + 1
#     y = y * 2
#     return foobar
# end
# @cacheable Dates.Minute(2) foo(arg1::Int, arg2::Int)::String = # ...

# const CACHE_foo_Int_String = Cache{Tuple{Int, String}, ReturnType}(timeout)
# function foo(args...)::ReturnType
#     return get!(CACHE_foo_Int_String, args) do
#         _foo(args...)
#     end
# end
# function _foo(arg1::Int, arg2::String)::ReturnType
#     # ...
# end

end # module


text:
function rotation_from_inertial(obj::LOFType, arg0::Field, arg1::FieldPVCoordinates)
    return jcall(obj, "rotationFromInertial", FieldRotation, (Field, FieldPVCoordinates), arg0, arg1)
end

function rotation_from_inertial(obj::LOFType, arg0::PVCoordinates)
    return jcall(obj, "rotationFromInertial", Rotation, (PVCoordinates,), arg0)
end

function transform_from_inertial(obj::LOFType, arg0::AbsoluteDate, arg1::PVCoordinates)
    return jcall(obj, "transformFromInertial", Transform, (AbsoluteDate, PVCoordinates), arg0, arg1)
end

function transform_from_inertial(obj::LOFType, arg0::FieldAbsoluteDate, arg1::FieldPVCoordinates)
    return jcall(obj, "transformFromInertial", FieldTransform, (FieldAbsoluteDate, FieldPVCoordinates), arg0, arg1)
end

function value_of(::Type{LOFType}, arg0::JString)
    return jcall(LOFType, "valueOf", LOFType, (JString,), arg0)
end

function values(::Type{LOFType})
    return jcall(LOFType, "values", Vector{LOFType}, ())
end


text:
using Pkg.Artifacts

# This is the path to the Artifacts.toml we will manipulate
artifact_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
isfile(artifact_toml) && rm(artifact_toml)

# Query the `Artifacts.toml` file for the hash bound to the name "blpapi"
# (returns `nothing` if no such binding exists)
let
    blpapi_hash = artifact_hash("blpapi", artifact_toml)
    @assert blpapi_hash == nothing
end

function create_artifact_windows()
    blpapi_hash = create_artifact() do artifact_dir
        dll_filename = "blpapi3_64.dll"
        local_dll = joinpath(@__DIR__, dll_filename)
        @assert isfile(local_dll) "$local_dll not found"
        cp(local_dll, joinpath(artifact_dir, dll_filename))
    end

    bind_artifact!(artifact_toml, "blpapi", blpapi_hash)
    blpapi_artifact_path = artifact_path(blpapi_hash)
    @info("Artifact folder: $blpapi_artifact_path")
    tar_sha = archive_artifact(blpapi_hash, joinpath(@__DIR__, "blpapi_cpp_3.12.3.1-windows-x64.tar.gz"))
    println("tar SHA: `$tar_sha`")
end

function create_artifact_linux()
    blpapi_hash = create_artifact() do artifact_dir
        dll_filename = "libblpapi3_64.so"
        local_dll = joinpath(@__DIR__, dll_filename)
        @assert isfile(local_dll) "$local_dll not found"
        cp(local_dll, joinpath(artifact_dir, dll_filename))
    end

    bind_artifact!(artifact_toml, "blpapi", blpapi_hash)
    blpapi_artifact_path = artifact_path(blpapi_hash)
    @info("Artifact folder: $blpapi_artifact_path")
    tar_sha = archive_artifact(blpapi_hash, joinpath(@__DIR__, "blpapi_cpp_3.12.3.1-linux-x64.tar.gz"))
    println("tar SHA: `$tar_sha`")
end

function create_artifact_macos()
    blpapi_hash = create_artifact() do artifact_dir
        dll_filename = "libblpapi3_64.so"
        local_dll = joinpath(@__DIR__, dll_filename)
        @assert isfile(local_dll) "$local_dll not found"
        cp(local_dll, joinpath(artifact_dir, dll_filename))
    end

    bind_artifact!(artifact_toml, "blpapi", blpapi_hash)
    blpapi_artifact_path = artifact_path(blpapi_hash)
    @info("Artifact folder: $blpapi_artifact_path")
    tar_sha = archive_artifact(blpapi_hash, joinpath(@__DIR__, "blpapi_cpp_3.14.3.1-macos-x64.tar.gz"))
    println("tar SHA: `$tar_sha`")
end

#create_artifact_windows()
#create_artifact_linux()
#create_artifact_macos()

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
#collection of bias correction routines

ma(K,N) = ma(K,N,2.0)
function ma(K, N, base)
    return (K-1)/(2*N*log(base))
end

blank(args...) = 0


text:
@doc raw"""
    sumBassLocal(p,q,m,c0,alpha,F,T,Pi,xstart)

This function locally maximizes a sum of Bass functions with a hot start at `xstart`
and returns the maximizer and maximum. It uses JuMP's interface with the IpOpt solver.
"""
function sumBassLocal(p,q,m,c0,alpha,F,T,Pi::Number,xstart)
    n = length(p) # Number of sites
    model = Model(optimizer_with_attributes(Ipopt.Optimizer,"print_level" => 0))
    @variable(model,x[1:n])
    set_start_value.(x,xstart)
    @constraint(model, F .<= x .<= T)
    # for i in 1:n
    #     set_start_value(x[i],node.x[i]) # Starting point at the maximizer found by SigmidalProgramming.jl
    #     @constraint(model, F[i] <= x[i] <= T[i])
    # end
    C = ones(1,n)
    @constraint(model,sum(C[i]*x[i] for i=1:n) == Pi)
    SumBass(x...) = sum(alpha[i]*Bass(x[i],p[i],q[i],m[i],c0[i]) for i=1:length(x))
    function SumBass_derivative(g,x...)
        for i in 1:length(x)
            g[i] = alpha[i]*Bass_derivative(x[i],p[i],q[i],m[i],c0[i])
            nothing
        end
    end
    register(model,:SumBass,n,SumBass,SumBass_derivative)
    @NLobjective(model,Max,SumBass(x...))
    # register(model,:Bass,5,Bass,autodiff=true)
    # @NLobjective(model,Max,sum(alpha[i]*Bass(x[i],p[i],q[i],m[i],c0[i]) for i=1:n))
    optimize!(model)
    maximizer = value.(x)
    maximum = objective_value(model)
    # if termination_status(model) == MOI.LOCALLY_SOLVED
    #     return value.(x), objective_value(model)
    # else
    #     error("The model was not solved correctly.")
    # end
    return maximizer,maximum
end

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
function Coeffgamm2_dVdq0(k, T, mfun)
    (
        k * (coth(Epi(k, mfun(k)) / (2 * T)) / Epi(k, mfun(k))^2)
    ) / (16 * pi^2)
end

function VImintqs_delta1_dq0(p0, ps, k, T, Npi,IRScale,UVScale, mfun, lamfun)
    -(2 + Npi) *
    π *
    3 *
    (
        deltasumkAll_dp0(p0 + Epi(k, mfun(k)), ps, k, T, Npi,IRScale, UVScale, mfun, lamfun) -
        deltasumkAll_dp0(p0 - Epi(k, mfun(k)), ps, k, T, Npi,IRScale, UVScale, mfun, lamfun)
    )
end

function VImintqs_delta1_dVdq0(
    p0,
    ps,
    k,
    T,
    Npi,
    IRScale,
    UVScale,
    mfun,
    lamfun,
)
    return -(2 + Npi) *
           π *
           3 *
           central_fdm(5, 1)(
               z -> (
                   deltasumkAllSmooth(
                       p0 + z,
                       ps,
                       k,
                       T,
                       Npi,
                       IRScale,
                       UVScale,
                       mfun,
                       lamfun,
                   ) + deltasumkAllSmooth(
                       p0 - z,
                       ps,
                       k,
                       T,
                       Npi,
                       IRScale,
                       UVScale,
                       mfun,
                       lamfun,
                   )
               ),
               Epi(k, mfun(k)),
           )
end

function fastpropImintqs_dVdq0(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            central_fdm(5, 1)(
                z -> dkVImintqs(
                    p0,
                    ps,
                    z,
                    x[1],
                    x[2] * (x[1] - UVScale) + UVScale,
                    mfun(x[2] * (x[1] - UVScale) + UVScale),
                    T,
                    Npi,
                    lamfun(x[2] * (x[1] - UVScale) + UVScale),
                ),
                Epi(x[1], mfun(x[1])),
            ) *
            Coeffgamm2_dVdq0(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1]
end

function propImintqs_delta1_dVdq0(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
    return -hquadrature(
        k ->
            2 *
            VImintqs_delta1_dVdq0(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun) *
            Coeffgamm2_dVdq0(k, T, mfun),
        IRScale,
        UVScale,
        rtol = 1e-6,
        atol = 1e-6,
        initdiv=20,
        maxevals = 8000,
    )[1]
end

function fastpropImintqs_All_dVdq0(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return fastpropImintqs_dVdq0(
        p0,
        ps,
        T,
        IRScale,
        UVScale,
        Npi,
        mfun,
        lamfun;
        kwarg...,
    ) + propImintqs_delta1_dVdq0(
        p0,
        ps,
        T,
        IRScale,
        UVScale,
        Npi,
        mfun,
        lamfun,
    )
end

@doc raw"""
    propImintqs_delta1_dq0_delta(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)

Compute the delta function part appears
"""
function propImintqs_delta1_dq0_delta1(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun,
)
    deltaf1(x) = Epi(x, mfun(x)) - p0
    if deltaf1(IRScale) * deltaf1(UVScale) >= 0
        return 0.0
    else
        k01 = find_zero(deltaf1, (IRScale, UVScale))
        return (2 + Npi) *sqrt(k01^2+mfun(k01))/(2*k01+derivative(mfun,k01))*
               π *
               3 *
               2 *deltasumkfix_dp0_deltafun(p0+Epi(k01,mfun(k01)),
                   ps,
                   k01,
                   k01,
                   T,
                   Npi,
                   IRScale,
                   UVScale,
                   mfun,
                   lamfun,
               ) * Coeffgamm2_dVdq0(k01, T, mfun)
    end
end

function propImintqs_delta1_dq0_delta2(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun,
)
    deltaf2(x) = 3 * Epi(x, mfun(x)) - p0
    if deltaf2(IRScale) * deltaf2(UVScale) >= 0
        return 0.0
    else
        k02 = find_zero(deltaf2, (IRScale, UVScale))
        return -(2 + Npi) *sqrt(k02^2+mfun(k02))/(2*k02+derivative(mfun,k02))*
               π *
               3 *
               2 *deltasumkfix_dp0_deltafun(p0-Epi(k02,mfun(k02)),
                   ps,
                   k02,
                   k02,
                   T,
                   Npi,
                   IRScale,
                   UVScale,
                   mfun,
                   lamfun,
               ) * Coeffgamm2_dVdq0(k02, T, mfun)
    end
end

function propImintqs_delta1_dq0(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
    return -hquadrature(
        k ->
            2 *
            VImintqs_delta1_dq0(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun) *
            Coeffgamm2_dVdq0(k, T, mfun),
        IRScale,
        UVScale,
        rtol = 1e-6,
        atol = 1e-6,
        maxevals = 2000,
    )[1]
end

function fastpropImintqs_All_dq0(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return fastpropImintqs_dVdq0(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun;
               kwarg...,
           ) +
           propImintqs_delta1_dq0(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun,
           ) +
           propImintqs_delta1_dq0_delta1(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun,
           )+propImintqs_delta1_dq0_delta2(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun,
           )
end


text:
import Base.read, Base.eof, Base.close

type OpenSecretsBuffer <: IO
    io::IO

    OpenSecretsBuffer(io::IO) = new(io)
end

function Base.read(os::OpenSecretsBuffer, ::Type{Uint8})
    c = Base.read(os.io, Uint8)
    c == '\\' ? ' ' : c
end

Base.eof(os::OpenSecretsBuffer) = eof(os.io)

Base.close(os::OpenSecretsBuffer) = close(os.io)


text:
function ueberpruefeFormel(Formel,Aufloesung)
#save(File(format"JLD","SuS/AB3.ipynb"));
Phantom = echteCTDaten(Aufloesung);
nrows,ncols = size(Phantom);
PhantomVec = reshape(Phantom,(Int(nrows)*Int(ncols),1));

if Aufloesung != 2 && Aufloesung != 4 && Aufloesung != 8 && Aufloesung != 16 && Aufloesung != 32 && Aufloesung != 64 && Aufloesung != 128 && Aufloesung != 256 && Aufloesung != 512
    error("Die Aufloesung darf nur zwischen 3 und 100 liegen.")
end

x0,y0,N,JJ,nx,ny,x1,y1,verschx,verschy,a,b,A,theta1,counter = Dateneinlesen(Aufloesung);
A2 = MartinStrahlen(Aufloesung,x0,x1,y0,y1,N,JJ,nx,ny,verschx,verschy);
Messwerte =  A2*PhantomVec;

Arows,Acols = size(A2);
f_nah = zeros(Int(Acols),1).+1;

if almostequal(Formel(A2,f_nah,Messwerte),norm(A2*f_nah-Messwerte)^2)
    println("Eure Formel für den Abstand ist korrekt!");
elseif almostequal(Formel(A2,f_nah,Messwerte),norm(A2*f_nah-Messwerte))
    println("Eure Formel ist noch nicht ganz korrekt. Da ihr hier einen Abstand berechnet und eigentliche dadurch eine Wurzel beim auflösen der Formel auftritt, quadriert man die rechte Seite.");
elseif almostequal(Formel(A2,f_nah,Messwerte),A2*f_nah-Messwerte)
    println("Eure Formel ist noch nicht ganz korrekt. Ihr müsst den Abstand zwischen den beiden Vektoren bestimmen.");
else
    println("Eure Formel stimmt noch nicht.");
end
Schritt = 1.1;
return  Messwerte,A2,PhantomVec,Arows,Acols,Schritt;

end


text:
using Plots

N = 1000
dt = 0.001
x0 = 1.0
p = scatter()
for i in 1:N - 1
    t = i * dt
    x = (1 - dt) * x0
    scatter!((t, x), label=:none)
    x0 = x
end

p
savefig("Julia-file/output/result_non-array.png")


text:
@use GADT

import Base: convert

struct Fun{T, R}
    fn :: Function
end

function (typed_fn :: Fun{T, R})(arg :: T) :: R where {T, R}
    typed_fn.fn(arg)
end

function convert(::Type{Fun{T, R}}, fn :: Function) where {T, R}
    Fun{T, R}(fn)
end

function convert(::Type{Fun{T, R}}, fn :: Fun{C, D}) where{T, R, C <: T, D <: R}
    Fun{T, R}(fn.fn)
end

⇒(::Type{A}, ::Type{B}) where {A, B} = Fun{A, B}

@data public Exp{T} begin
    Sym :: Symbol => Exp{A} where {A}
    Val{A} :: A => Exp{A}
    App{A, B, A_ <: A} :: (Exp{Fun{A, B}}, Exp{A_}) => Exp{B}
    Lam{A, B} :: (Symbol, Exp{B}) => Exp{Fun{A, B}}
    If{A} :: (Exp{Bool}, Exp{A}, Exp{A}) => Exp{A}
end

function substitute(template :: Exp{T}, pair :: Tuple{Symbol, Exp{G}}) where {T, G}
    (sym, exp) = pair
    @match template begin
        Sym(&sym) => exp
        Val(_) => template
        App(f, a) => App(substitute(f, pair), substitute(a, pair)) :: Exp{T}
        Lam(&sym, exp) => template
        If(cond, exp1, exp2) =>
            let (cond, exp1, exp2) = map(substitute, (cond, exp1, exp2))
                If(cond, exp1, exp2) :: Exp{T}
            end
    end
end

function eval_exp(exp :: Exp{T}, ctx :: Dict{Symbol, Any}) where T
    @match exp begin
        Sym(a) => (ctx[a] :: T, ctx)
        Val(a :: T) => (a, ctx)
        App{A, T, A_}(f :: Exp{Fun{A, T}}, arg :: Exp{A_}) where {A, A_ <: A} =>
            let (f, ctx) = eval_exp(f, ctx),
                (arg, ctx) = eval_exp(arg, ctx)
                (f(arg), ctx)
            end
        Lam{A, B}(sym, exp::Exp{B}) where {A, B} =>
            let f(x :: A) = begin
                    A
                    eval_exp(substitute(exp, sym => Val(x)), ctx)[1]
                end

                (f, ctx)
            end
        If(cond, exp1, exp2) =>
            let (cond, ctx) = eval_exp(cond, ctx)
                eval_exp(cond ? exp1 : exp2, ctx)
            end
    end
end

add = Val{Number ⇒ Number ⇒ Number}(x -> y -> x + y)
sub = Val{Number ⇒ Number ⇒ Number}(x -> y -> x - y)
gt = Val{Number ⇒ Number ⇒ Bool}(x -> y -> x > y)
ctx = Dict{Symbol, Any}()

@assert 3 == eval_exp(App(App(add, Val(1)), Val(2)), ctx)[1]
@assert -1 == eval_exp(App(App(sub, Val(1)), Val(2)), ctx)[1]
@assert 1 == eval_exp(
    If(
        App(App(gt, Sym{Int}(:x)), Sym{Int}(:y)),
        App(App(sub, Sym{Int}(:x)), Sym{Int}(:y)),
        App(App(sub, Sym{Int}(:y)), Sym{Int}(:x))
    ), Dict{Symbol, Any}(:x => 1, :y => 2))[1]


text:
# AUTO GENERATED FILE - DO NOT EDIT

export cardsection

"""
    cardsection(;kwargs...)
    cardsection(children::Any;kwargs...)
    cardsection(children_maker::Function;kwargs...)

A CardSection component.
Card.Section is a special component that is used to remove Card padding from its children while other elements still have horizontal spacing. For more information, see: https://mantine.dev/core/card/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Section content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
"""
function cardsection(; kwargs...)
        available_props = Symbol[:children, :id, :className]
        wild_props = Symbol[]
        return Component("cardsection", "CardSection", "dash_mantine_components", available_props, wild_props; kwargs...)
end

cardsection(children::Any; kwargs...) = cardsection(;kwargs..., children = children)
cardsection(children_maker::Function; kwargs...) = cardsection(children_maker(); kwargs...)


text:
export greedysolver

function greedysolver(algo::String, ψ::AbstractMatrix{T}, u::AbstractVector{T}; invert::Bool=true, verbose::Bool = true, ϵrel::Float64 = 1e-1, maxterms::Int64=typemax(Int64)) where {T}
    if algo ∈ ["pivotedqr", "pivot"]
        return pivotedqr(ψ, u; invert = invert, verbose = verbose, ϵrel = ϵrel, maxterms = maxterms)
    elseif algo ∈ ["omp", "lsomp"]
        return lsomp(ψ, u; invert = invert, verbose = verbose, ϵrel = ϵrel, maxterms = maxterms)
    elseif algo ∈ ["qromp"]
        return qromp(ψ, u; invert = invert, verbose = verbose, ϵrel = ϵrel, maxterms = maxterms)
    end
end

max_forks_repo_name: mleprovost/QROMP


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libwebp_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libwebp")
JLLWrappers.@generate_main_file("libwebp", UUID("c5f90fcd-3b7e-5836-afba-fc50a0988cb2"))
end  # module libwebp_jll


text:
module CollabUtils2MiniHeartDoorUnlockCutsceneTrigger

using ..Ahorn, Maple

@mapdef Trigger "CollabUtils2/MiniHeartDoorUnlockCutsceneTrigger" MiniHeartDoorUnlockCutsceneTrigger(x::Integer, y::Integer,
    width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, doorID::String="")

const placements = Ahorn.PlacementDict(
    "Mini Heart Door Unlock Cutscene (Collab Utils 2 / READ DOCS)" => Ahorn.EntityPlacement(
        MiniHeartDoorUnlockCutsceneTrigger,
        "rectangle",
    ),
)

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Psi4 $\leftrightarrow$ Julia Data Sharing
#
# The heart of the Psi4Julia project its the ability to easily share and manipulate quantities in Julia. While Psi4 offers the ability to manipulate most objects and perform tensor operations at the Python layer, it is often much easier to use Julia. Fortunately, Psi4 offers seemless integration with Julia. More details on the underlying functions can be found in the Psi4 [documentation](http://psicode.org/psi4manual/master/numpy.html).
#
# As before, let us start off with importing Psi4 and NumPy while also creating a random `5 x 5` Julia array:

# +
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy")

# Random number array
array = rand(5, 5)
# -

# Converting this to a Psi4 Matrix, which is an instance of the [`psi4.core.Matrix`](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Matrix
# "Go to API") class, and back again is as simple as:

# +
psi4_matrix = psi4.core.Matrix.from_array(array)
new_array = np.array(psi4_matrix)

println("Allclose new_array, array: ", np.allclose(new_array, array))
# -

# ## Views
# Because both of these objects have the same in-memory data layout, the conversion is accomplished using the NumPy
# [array_interface](https://docs.scipy.org/doc/numpy/reference/arrays.interface.html). This also opens the opportunity
# to manipulate the Psi4 Matrix and Vector classes directly in memory.  To do this, we wrap the PyObject with a Julia Matrix object:

# +
matrix = psi4.core.Matrix(3, 3)
print("Zero Psi4 Matrix:")
display(np.array(matrix))

#matrix.np[:] .= 1 # does not work as in Python, it is read-only
#matrix.set(1)     # works but it uses psi.core, not julia
function psi4view(psi4matrix)
   # Limitations: you have to check the type manually, it may not be a matrix
   # Use numpy __array_interface__
   # ["data"][1] is the pointer to the first element of the array
   # ["data"][2] is true when it is read-only
   # ["strides"] is nothing if has C ordering (default)
   # ["shape"]   has the shape in (C ordering by default) [needs reversing]
   # ["typestr"] has the element type
   array_interface = psi4matrix.__array_interface__
   array_interface["data"][2] == false   || @warn "Not writable"
   array_interface["strides"] == nothing || @warn "Different ordering than C"
   ptr   = array_interface["data"][1]
   shape = reverse(array_interface["shape"])
   unsafe_wrap(Matrix{Float64}, Ptr{Float64}(ptr), shape)
end
jlmatrix = psi4view(matrix) # This is a Julia matrix using psi4 memory
jlmatrix .= 1
print("Matrix updated to ones:")
display(np.asarray(matrix))
# -

# Theoretically, `PyArray(matrix)` (from PyCall) should return a Julia array. But here it does not work.
# The `psi4view` function created here effectively returns a Julia matrix that uses the memory of the Psi4 object. This view can then be manipulated as a conventional NumPy array and the underlying Psi4 Matrix data will be modified.
#
# **Warning!** The following operation operation is incorrect and can potenitally lead to confusion:

display(psi4.core.Matrix(3, 3).np)

# While the above operation works about ~90% of the time, occasionally you will notice extremely large and small values. This is due to the fact that when you create the Psi4 Matrix and grab its view, the Psi4 Matrix is no longer bound to anything, and Python will attempt to "garbage collect" or remove the object. This sometimes happens *before* Python prints out the object so the NumPy view is pointing to a random piece of data in memory. However, this has not been observed in Julia. A safe way to do this would be:

# +
mat = psi4.core.Matrix(3, 3)
display(mat.np)

# or
display(np.asarray(psi4.core.Matrix(3, 3)))
# -

# Similar to the `.np` attribute, one can use `np.asarray` to create a Julia copy of a Psi4 object. `np.array` which will copy the data too.

# +
mat = psi4.core.Matrix(3, 3)
mat_view = np.asarray(mat)

mat_view .= rand(mat.shape...)
display(mat.np)

# +
mat = psi4.core.Matrix(3, 3)
mat_view = psi4view(mat)

mat_view .= rand(mat.shape...)
display(mat.np)
# -

# Keep in mind that you must *update* this view using the `.=` syntax and not replace it (`=`). The following example should demonstrate the difference:

# +
mat_view = zeros(3, 3)

# Mat is not updated as we replaced the mat_view with a new Julia matrix.
display(mat.np)
# -

# ## Vector class
# Like the Psi4 Matrix class, the [`psi4.core.Vector`](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.Vector "Go to API")
# class has similar accessors:

arr = rand(5)
vec = psi4.core.Vector.from_array(arr)
display(vec.np)


text:
#=
    2019.3.20
    wechat:  Cyvadra
    mail:    dj@lifesdk.com
    License: Apache 2.0
    版权所有: 上海禀土科技发展有限公司

=#

using Dates
using SoftGlobalScope

function generate_date(t_patch::Dict{} = time_patch)
    Dates.DateTime( t_patch["year"],t_patch["month"],t_patch["day"],t_patch["hour"],t_patch["minute"] )
end

function time_standardization(t_patch::Dict{})
    while t_patch["minute"] >= 60
        t_patch["hour"]   += 1
        t_patch["minute"] -= 60
    end
    while t_patch["hour"] >= 24
        t_patch["day"]  += 1
        t_patch["hour"] -= 24
    end
    @softscope while t_patch["day"] > Dates.day( Dates.lastdayofmonth( Dates.Date(t_patch["year"],t_patch["month"],1) ) )
        t_patch["day"]   -= Dates.day( Dates.lastdayofmonth( Dates.Date(t_patch["year"],t_patch["month"],1) ) )
        t_patch["month"] += 1
    end
    t_patch
end

function _into_num(pure_number_pattern)
    try
        if string( parse(Int,pure_number_pattern) ) == string(pure_number_pattern)
            return parse(Int64,pure_number_pattern)
        end
    catch
        tmp_first, tmp_sum = 1, 0
        d_nums  = Dict( '半' => 0.5, '零' => 0,  '一' => 1,  '二' => 2, '两' => 2, '三' => 3,  '四' => 4, '五' => 5,  '六' => 6, '七' => 7,  '八' => 8, '九' => 9,'0' => 0,  '1' => 1,'2' => 2,'3' => 3, '4' => 4,'5' => 5,'6' => 6, '7' => 7,'8' => 8, '9' => 9 )
        d_units = Dict( '个' => 1, '十' => 10, '百' => 100 )
        for char in pure_number_pattern
            if haskey(d_units,char)
                tmp_sum += d_units[char]*tmp_first
                tmp_first= 0
            end
            if haskey(d_nums,char)
                tmp_first = d_nums[char]
            end
        end
        if tmp_first !== 0
            tmp_sum += tmp_first
        end
        return tmp_sum
    end
end

function calc_target_time(input_time_string="明天八点半叫我起床")

    # data patch
        current_time       = now()
        time_patch         = Dict(  "year" => Dates.year(current_time),
                                    "month"=> Dates.month(current_time),
                                    "day"  => Dates.day(current_time),
                                    "hour" => Dates.hour(current_time),
                                    "minute"=>Dates.minute(current_time)
                                )
        current_time       = generate_date(time_patch)
        default_hour       = 12
        default_morning    = 9
        default_afternoon  = 14
        default_evening    = 19

    # correction
        input_time_string  = replace(input_time_string,r"([早晚天])+儿" => s"\g<0>")
        input_time_string  = replace(input_time_string,"明早" => "明天早上")
        input_time_string  = replace(input_time_string,"明晚" => "明天晚上")
        input_time_string  = replace(input_time_string,"今晚" => "今天晚上")

    # month & date
        pattern_next_month = r"(下[个]?月)+([0-9一二三]?十?[0-9一二三四五六七八九]?[号日])+"
            if match(pattern_next_month,input_time_string) !== nothing
                time_patch["month"] += 1
                time_patch["day"]    = 1
                time_patch["hour"] = default_hour; time_patch["minute"] = 0;
            end
        pattern_month      = r"(([0-9一二三四五六七八九十]+[一二]?)月)?([0-9一二三]?十?[0-9一二三四五六七八九十]+)[号日]+"
            if match(pattern_month,input_time_string) !== nothing
                tmp_month = match(pattern_month,input_time_string)[2]
                tmp_day   = match(pattern_month,input_time_string)[4]
                tmp_month !== nothing ? time_patch["month"] = _into_num(tmp_month) : false
                tmp_day   !== nothing ? time_patch["day"]  = _into_num(tmp_day)    : false
                time_patch["hour"] = default_hour; time_patch["minute"] = 0;
            end

    # week
        pattern_next_week  = r"下个?[周礼拜星期]+"
            if match(pattern_next_week,input_time_string) !== nothing
                @softscope while Dates.week( Dates.Date(time_patch["year"],time_patch["month"],time_patch["day"]) ) == Dates.week(Dates.Date(current_time)) # should == current+1
                    time_patch["day"] += 1
                    time_patch = time_standardization(time_patch)
                end
            end
        pattern_week       = r"[周礼拜星期]+([一二三四五六天日])+"
            if match(pattern_week,input_time_string) !== nothing
                tmp_weekday = match(pattern_week,input_time_string)[1]
                weekday_list= Dict( "一"=>Dates.ismonday,"二"=>Dates.istuesday,"三"=>Dates.iswednesday,"四"=>Dates.isthursday,"五"=>Dates.isfriday,"六"=>Dates.issaturday,"天"=>Dates.issunday,"日"=>Dates.issunday )
                @softscope while weekday_list[tmp_weekday]( Dates.Date(time_patch["year"],time_patch["month"],time_patch["day"]) ) == false
                    time_patch["day"] += 1
                    time_patch = time_standardization(time_patch)
                end
                time_patch["hour"] = default_hour; time_patch["minute"] = 0;
            end

    # day afterwards
        if match(r"[明后]天",input_time_string) !== nothing
            occursin("明天",input_time_string)   ? time_patch["day"] += 1 : false
            occursin("后天",input_time_string)   ? time_patch["day"] += 2 : false
            occursin("大后天",input_time_string) ? time_patch["day"] += 1 : false # extra
            time_patch["hour"] = default_hour; time_patch["minute"] = 0;
        end
        pattern_day_later  = r"([0-9零一二三四五六七八九十]?[0-9零一二三四五六七八九十]?)+[天日]+后"
            if match(pattern_day_later,input_time_string) !== nothing
                time_patch["day"] += _into_num( match(pattern_day_later,input_time_string)[1] )
                time_patch = time_standardization(time_patch)
            end

    # hour
        pattern_next_hour  = r"([0-9一二三四五六七八九十]个?半?)小时后"
            if match(pattern_next_hour,input_time_string) !== nothing
                time_patch["minute"] += _into_num( match(pattern_next_hour,input_time_string)[1] )*60
                #time_standardization()
            end
        pattern_hour       = r"([0-9零一二三四五六七八九十]+)点钟?"
            if match(r"中午",input_time_string) !== nothing
                time_patch["hour"] = 12
            end
        pattern_morning    = r"(早上|上午|白天|凌晨)+"
        pattern_afternoon  = r"(下午|晚上|夜里)+"
            if match(pattern_morning,input_time_string) !== nothing
                time_patch["hour"]   = default_morning
                time_patch["minute"] = 0;
            elseif match(r"(下午)+",input_time_string) !== nothing
                time_patch["hour"]   = default_afternoon
                time_patch["minute"] = 0;
            elseif match(r"(晚上|夜里)+",input_time_string) !== nothing
                time_patch["hour"]   = default_evening
                time_patch["minute"] = 0;
            end
        if match(pattern_hour,input_time_string) !== nothing
            if match(pattern_morning,input_time_string) !== nothing
                tmp_hour   = _into_num( match(pattern_hour,input_time_string)[1] )
            elseif match(pattern_afternoon,input_time_string) !== nothing
                tmp_hour   = _into_num( match(pattern_hour,input_time_string)[1] ) + 12
                tmp_hour >= 24 ? tmp_hour -= 12 : tmp_hour
            else
                tmp_hour   = _into_num( match(pattern_hour,input_time_string)[1] )
            end
            time_patch["hour"]   = tmp_hour
            time_patch["minute"] = 0;
            time_patch           = time_standardization(time_patch)
        end

    # minute
        pattern_hour_half = r"点半"
            if match(pattern_hour_half,input_time_string) !== nothing
                time_patch["minute"] = 30
            end
        pattern_quarter   = r"([1-3一二三])刻钟?^后"
            if match(pattern_quarter,input_time_string) !== nothing
                time_patch["minute"] = _into_num( match(pattern_quarter,input_time_string) )*15
            end
        pattern_minute    = r"[0-9零一二三四五六七八九十]+点钟?([1-9一二三四五六七八九十]?[0-9一二三四五六七八九十]+)"
            if match(pattern_minute,input_time_string)  !== nothing
                time_patch["minute"] = _into_num( match(pattern_minute,input_time_string)[1] )
            end
        pattern_min_later = r"([0-9一二三四五六七八九十]+)分钟后"
            if match(pattern_min_later,input_time_string)  !== nothing
                time_patch["minute"] += _into_num( match(pattern_min_later,input_time_string)[1] )
            end

    # standardization
        target_time_patch = time_standardization( time_patch )
        target_time = generate_date(target_time_patch)

    if target_time == current_time
        return Dict( "hasTime"=>false,"DateTime"=>current_time )
    else
        return Dict( "hasTime"=>true, "DateTime"=>target_time  )
    end

end

max_stars_repo_name: Cyvadra/julia-cnnlp-toolkit

max_issues_repo_name: Cyvadra/julia-cnnlp-toolkit

max_forks_repo_name: Cyvadra/julia-cnnlp-toolkit


text:
## Show the calculation of Gaussian and Airy PSF models

#include("../src/BayesMultiFit.jl")
#import BAMF: RJMCMC
using BayesMultiFit
BAMF=BayesMultiFit
using ImageView
using Plots
ImageView.closeall()

## Setup
sz=Int32(32)

## setup point source
n=Int32(1)
x=[Float32(sz/2f0)]
y=[Float32(sz/2f0)]
photons=[1f0]
bg=0f0
pointsource=BAMF.StateFlatBg(n,x,y,photons,bg)

## 2D Gaussian PSF
σ=Float32(1.3)
psf=BAMF.PSF_gauss2D(σ)
gausspsf=BAMF.ArrayDD(sz)
BAMF.genmodel!(pointsource,psf,gausspsf)
imshow(gausspsf.data)
println(("integrated value: " ,sum(gausspsf.data)))

## 2D Airy PSF
pixelsize=.05
nₐ=1.4
λ=.6
ν=Float32(2π*nₐ/λ)*pixelsize
psf=BAMF.PSF_airy2D(ν)
airypsf=BAMF.ArrayDD(sz)
BAMF.genmodel!(pointsource,psf,airypsf)
imshow(airypsf.data)
println(("integrated value: " ,sum(airypsf.data)))


text:
function LatitudeExtremumDetector(arg0::OneAxisEllipsoid)
    return LatitudeExtremumDetector((OneAxisEllipsoid,), arg0)
end

function LatitudeExtremumDetector(arg0::jdouble, arg1::jdouble, arg2::OneAxisEllipsoid)
    return LatitudeExtremumDetector((jdouble, jdouble, OneAxisEllipsoid), arg0, arg1, arg2)
end

function g(obj::LatitudeExtremumDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_body(obj::LatitudeExtremumDetector)
    return jcall(obj, "getBody", BodyShape, ())
end


text:
type BinaryNode{T}
    left::Nullable{T}
    right::Nullable{T}
    BinaryNode()                  = new(Nullable{T}(), Nullable{T}())
    BinaryNode(left::T, right::T) = new(Nullable{T}(a), Nullable{T}(b))
end
type RectanglePacker{T}
    children::BinaryNode{RectanglePacker{T}}
    area::SimpleRectangle{T}
end

left(a::RectanglePacker)                                = get(a.children.left)
left{T}(a::RectanglePacker{T}, r::RectanglePacker{T})   = (a.children.left = Nullable(r))
right(a::RectanglePacker)                               = get(a.children.right)
right{T}(a::RectanglePacker{T}, r::RectanglePacker{T})  = (a.children.right = Nullable(r))
RectanglePacker{T}(area::SimpleRectangle{T})            = RectanglePacker{T}(BinaryNode{RectanglePacker{T}}(), area)
isleaf(a::RectanglePacker)                              = isnull(a.children.left) && isnull(a.children.right)

# This is rather append, but it seems odd to use another function here.
# Maybe its a bad idea, to call it push regardless!?
function Base.push!{T}(node::RectanglePacker{T}, areas::Vector{SimpleRectangle{T}})
    sort!(areas)
    RectanglePacker{T}[push!(node, area) for area in areas]
end
function Base.push!{T}(node::RectanglePacker{T}, area::SimpleRectangle{T})
    if !isleaf(node)
        l = push!(left(node), area)
        l == nothing && return push!(right(node), area) # if left does not have space, try right
        return l
    end
    newarea = RectanglePacker(area).area
    if newarea.w <= node.area.w && newarea.h <= node.area.h
        oax,oay,oaxw,oayh   = node.area.x+newarea.w, node.area.y, xwidth(node.area), node.area.y + newarea.h
        nax,nay,naxw,nayh   = node.area.x, node.area.y+newarea.h, xwidth(node.area), yheight(node.area)
        rax,ray,raxw,rayh   = node.area.x, node.area.y, node.area.x+newarea.w, node.area.y+newarea.h
        left(node,  RectanglePacker(SimpleRectangle(oax, oay, oaxw-oax, oayh-oay)))
        right(node, RectanglePacker(SimpleRectangle(nax, nay, naxw-nax, nayh-nay)))
        return RectanglePacker(SimpleRectangle(rax,ray,raxw-rax,rayh-ray))
    end
end


text:
"""
	grdtrack(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

Interpolates the grid(s) at the positions in the table and returns the table with the
interpolated values added as (one or more) new columns.

Full option list at [`grdtrack`]($(GMTdoc)/grdtrack.html)

Parameters
----------

- **A** | **interp_path** :: [Type => Str]		`Arg = f|p|m|r|R[+l]`

    For track resampling (if `crossprofile` or `profile` are set) we can select how this is to be performed.
    ($(GMTdoc)grdtrack.html#a)
- **C** | **crossprofile** :: [Type => Str]		`Arg = length/ds[/spacing][+a|+v][l|r]`

    Use input line segments to create an equidistant and (optionally) equally-spaced set of crossing
    profiles along which we sample the grid(s)
    ($(GMTdoc)grdtrack.html#c)
- **D** | **dfile** :: [Type => Str]

    In concert with `crossprofile` we can save the (possibly resampled) original lines to the file dfile
    ($(GMTdoc)grdtrack.html#d)
- **E** | **profile** :: [Type => Str]

    Instead of reading input track coordinates, specify profiles via coordinates and modifiers.
    ($(GMTdoc)grdtrack.html#e)
- **F** | **critical** :: [Type => Str]

    Find critical points along each cross-profile as a function of along-track distance. Requires
    `crossprofile` and a single input grid.
    ($(GMTdoc)grdtrack.html#e)
- **G** | **grid** :: [Type => Str | GMTgrid | Tuple(GMTgrid's)]

    ($(GMTdoc)grdtrack.html#g)
- **N** | **no_skip** :: [Type => Bool]

    ($(GMTdoc)grdtrack.html#n)
- $(GMT.opt_R)
- **S** | **stack** :: [Type => Str]

    ($(GMTdoc)grdtrack.html#s)
- **T** | **radius** :: [Type => Number, Str | []]

    ($(GMTdoc)grdtrack.html#t)
- **Z** | **z_only** :: [Type => Bool]

    ($(GMTdoc)grdtrack.html#z)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_bo)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_n)
- $(GMT.opt_o)
- $(GMT.opt_s)
- $(GMT.opt_w)
- $(GMT.opt_swap_xy)

When using two numeric inputs and no G option, the order of the x,y and grid is not important.
That is, both of this will work: D = grdtrack([0 0], G);  or  D = grdtrack(G, [0 0]);
"""
function grdtrack(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	length(kwargs) == 0 && arg1 === nothing && return monolitic("grdtrack", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:R :V_params :bi :bo :di :e :f :g :h :i :n :o :s :w :yx])
	cmd  = parse_these_opts(cmd, d, [[:A :interp_path], [:C :crossprofile :equidistant], [:D :dfile], [:E :profile],
	                                 [:F :critical], [:N :no_skip], [:S :stack], [:T :radius], [:Z :z_only]])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)
	cmd, grid_tuple, arg1, arg2 = parse_G_grdtrk(d, [:G, :grid], cmd, arg1, arg2)

	# Because we allow arg1 and arg2 to either exist or not and also contain data & grid in any order
	if (arg1 !== nothing && arg2 !== nothing)
		arg2_is_table = false;	arg1_is_grid = false
		(isa(arg1, GMTgrid)) && (arg1_is_grid = true)
		(isa(arg2, Array) || isa(arg2, GMTdataset)) && (arg2_is_table = true)
		if (arg2_is_table && arg1_is_grid)			# Swap the arg1, arg2
			tmp = arg1;		arg1 = arg2;	arg2 = tmp
		end
	end

	(isa(arg1, GMTgrid) || isa(arg2, GMTgrid) && !occursin("-G", cmd)) && (cmd *= " -G")
	if (cmd0 != "" && !isa(arg1, GMTgrid) && !isa(arg1, Vector{GMTgrid}) && !occursin("-G", cmd))
		cmd = " -G" * cmd
		(isa(arg1, String)) && (cmd *= " " * arg1)
	end

	if (isa(grid_tuple, Tuple))
		return common_grd(d, "grdtrack " * cmd, (got_fname != 0) ? grid_tuple : tuple(arg1, grid_tuple...))
	else
		return common_grd(d, "grdtrack " * cmd, arg1, arg2)
	end

end

# ---------------------------------------------------------------------------------------------------
function parse_G_grdtrk(d::Dict, symbs::Vector{<:Symbol}, cmd::String, arg1, arg2)

	(show_kwargs[1]) && return (print_kwarg_opts(symbs, "GMTgrid | Tuple | String"), nothing,arg1,arg2)

	if ((grid_tuple = find_in_dict(d, symbs)[1]) !== nothing)
		if (isa(grid_tuple, Tuple))
			for k = 1:length(grid_tuple)  cmd *= " -G"  end		# Need as many -G as numel(grid_tuple)
		elseif (isa(grid_tuple, GMTgrid))
			cmd *= " -G"
			if     (arg1 === nothing)  arg1 = grid_tuple;
			elseif (arg2 === nothing)  arg2 = grid_tuple;
			else   error("Can't send the Grid data via G and input array")
			end
		else
			cmd = string(cmd, " -G", arg2str(grid_tuple))
		end
	end
	return cmd, grid_tuple, arg1, arg2
end

# ---------------------------------------------------------------------------------------------------
grdtrack(arg1, arg2=nothing, cmd0::String=""; kw...) = grdtrack(cmd0, arg1, arg2; kw...)


text:
module LazIO

using Libdl
using FileIO
using LasIO
using FixedPointNumbers

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LazIO not installed properly, run Pkg.build(\"LazIO\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
    # temporary _ until LasIO defers this key
    add_format(format"LAZ_", (), ".laz", [:LazIO])
end

include("utils.jl")
include("laszip_h.jl")
include("laszip.jl")
include("fileio.jl")
include("convert.jl")
include("dataset.jl")
include("write.jl")
include("table.jl")

end # module


text:
"""
    flattenbatch(x::AbstractArray)
flatten a multi dimensional array to keep only the last dimension.
It returns a 2 dimensional array of size (flatten_dim, batch_size)
"""
function flattenbatch(x::AbstractArray)
    reshape(x, (:, size(x)[end]))
end

"""
    huber_loss(x)
Compute the Huber Loss (from ReinforcementLearning.jl)
"""
function huber_loss(x)
    abserror = abs.(x)
    quadratic = min.(abserror, one(x))
    linear = abserror .- quadratic
    return 0.5f0 .*quadratic .* quadratic .+ linear
end

"""
    isrecurrent(m)
returns true if m contains a recurrent layer
"""
function isrecurrent(m)
    for layer in m
        if layer isa Flux.Recur
            return true
        end
    end
    return false
end

"""
    globalnorm(p::Params, gs::Flux.Zygote.Grads)
returns the maximum absolute values in the gradients of W
"""
function globalnorm(ps::Flux.Params, gs::Flux.Zygote.Grads)
    gnorm = 0f0
    for p in ps
        gs[p] === nothing && continue
        curr_norm = maximum(abs.(gs[p]))
        gnorm =  curr_norm > gnorm  ? curr_norm : gnorm
    end
    return gnorm
end

"""
    batch_trajectories(s::AbstractArray, traj_length::Int64, batch_size::Int64)
converts multidimensional arrays into batches of trajectories to be process by a Flux recurrent model.
It takes as input an array of dimension state_dim... x traj_length x batch_size
"""
function batch_trajectories(s::AbstractArray, traj_length::Int64, batch_size::Int64)
    return Flux.batchseq([[s[axes(s)[1:end-2]..., j, i] for j=1:traj_length] for i=1:batch_size], 0)
end

"""
    hiddenstates(m)
returns the hidden states of all the recurrent layers of a model
"""
function hiddenstates(m)
    return [l.state for l in m if l isa Flux.Recur]
end

"""
    sethiddenstates!(m, hs)
Given a list of hiddenstate, set the hidden state of each recurrent layer of the model m
to what is in the list.
The order of the list should match the order of the recurrent layers in the model.
"""
function sethiddenstates!(m, hs)
    i = 1
    for l in m
        if isa(l, Flux.Recur)
            l.state = hs[i]
            i += 1
        end
    end
end


text:
using PolyChaosODE,
    BenchmarkTools, DifferentialEquations, Distributions, LaTeXStrings, Plots

### Test based on the example provided at https://julia.quantecon.org/continuous_time/seir_model.html ###

# Values randomly sampled from Uniform(0.0005, 0.001)
const _e_1 = 0.0005187461238971071
const _e_2 = 0.0007294624156978437
const _e_3 = 0.000980192256832965

const VAR_INDEX = 3

μ = 0.0
σ = sqrt(2.0) / 2

function seir!(du, u, p, x, t)
    γ = p[1] + _e_1 * x[1]
    σ = p[2] + _e_2 * x[2]
    R_0 = p[3] + _e_3 * x[3]
    du[1] = -γ * R_0 * u[1] * u[3]
    du[2] = γ * R_0 * u[1] * u[3] - σ * u[2]
    du[3] = σ * u[2] - γ * u[3]
    du[4] = γ * u[3]
end

vars = Tuple(Normal(μ, σ) for i = 1:3)
dim_ = 4
stoch_galerkin_ode = StochGalerkinODE(seir!, dim_, vars)

i_0 = 1E-7
e_0 = 4.0 * i_0
s_0 = 1.0 - i_0 - e_0
r_0 = 0.0

u0 = [s_0, e_0, i_0, r_0]
tspan = (0.0, 350.0)
p = [1 / 18, 1 / 5.2, 3.0]

interval_t_start = 300.0
tstep = 0.1
interval_t = interval_t_start:tstep:tspan[2]

sol = stoch_galerkin_ode(u0, tspan, p; alg = VCABM())

sobol_indices_pce = mapreduce(
    transpose,
    vcat,
    compute_total_order_sobol_indices(stoch_galerkin_ode, sol, interval_t, VAR_INDEX),
)[
    2:end,
    :,
]

sobol_indices_mc_1k = mc_gsa_estimate_total_indices(
    stoch_galerkin_ode,
    VAR_INDEX,
    1000,
    interval_t,
    u0,
    tspan,
    p,
)

sobol_indices_mc_20k = mc_gsa_estimate_total_indices(
    stoch_galerkin_ode,
    VAR_INDEX,
    20000,
    interval_t,
    u0,
    tspan,
    p,
)

plot(
    interval_t[2:end],
    ((x, y) -> log10.(abs.(x .- y) ./ y)).(sobol_indices_pce, sobol_indices_mc_20k);
    label = [L"\log_{10} err\left(\hat{S}_{1, 3}^{(PCE)}\right)" L"\log_{10} err\left(\hat{S}_{2, 3}^{(PCE)}\right)" L"\log_{10} err\left(\hat{S}_{3, 3}^{(PCE)}\right)" L"\log_{10} err\left(\hat{S}_{4, 3}^{(PCE)}\right)"],
    legend = :outertopright,
    size = (500, 400),
    ylims = (-5, 0),
)

xlabel!(L"t")
ylabel!(L"\log_{10}\left\vert \frac{\hat{S} - S}{S} \right\vert")

# savefig("seir_pce.png")

plot(
    interval_t[2:end],
    ((x, y) -> log10.(abs.(x .- y) ./ y)).(sobol_indices_mc_1k, sobol_indices_mc_20k),
    label = [L"\log_{10} err\left(\hat{S}_{1, 3}^{(MC)}\right)" L"\log_{10} err\left(\hat{S}_{2, 3}^{(MC)}\right)" L"\log_{10} err\left(\hat{S}_{3, 3}^{(MC)}\right)" L"\log_{10} err\left(\hat{S}_{4, 3}^{(MC)}\right)"],
    legend = :outertopright,
    size = (500, 400),
    ylims = (-5, 0),
)

xlabel!(L"t")
ylabel!(L"\log_{10}\left\vert \frac{\hat{S} - S}{S} \right\vert")

# savefig("seir_mc.png")

@btime stoch_galerkin_ode(u0, tspan, p; alg = VCABM())
@btime compute_total_order_sobol_indices(stoch_galerkin_ode, sol, interval_t, VAR_INDEX)
@btime mc_gsa_estimate_total_indices(
    stoch_galerkin_ode,
    VAR_INDEX,
    1000,
    interval_t,
    u0,
    tspan,
    p,
)

max_stars_repo_name: d-netto/PolyChaosODE

max_issues_repo_name: d-netto/PolyChaosODE

max_forks_repo_name: d-netto/PolyChaosODE


text:
using Logging
logger = SimpleLogger(stdout, Logging.Info)
global_logger(logger)

@info "Learn Julia 2D Array"
row_count = 5
col_count = 10
A = Array{Int}(undef, row_count, col_count)
B = fill(Int(0), row_count, col_count)

@info "Learn Julia Array of Array"
frac_count = 10
C = Array{Array{Int}}(undef, frac_count)
C[1] = [1, 2, 3]
C[2] = [-1,0,1]
for i in 3:length(C)
    C[i] = [-2,-1,0,1,2]
end

@info "Learn comprehension"
D = [[-13,-7,0,7,13] for i in 1:length(C)]

nothing


text:
using Dates

# Test data is 3 days after a starting reference on Jan 1:
# The actual phi phase is [0, 1, 2, 3]
#
#
geolist = [Date(2019, 1, 1), Date(2019, 1, 3), Date(2019, 1, 6), Date(2019, 1, 10)]
intlist = [(geolist[i], geolist[j]) for i = 1:length(geolist)-1 for j = i+1:length(geolist)]

# julia> intlist
#  (2019-01-01, 2019-01-03)
#  (2019-01-01, 2019-01-06)
#  (2019-01-01, 2019-01-10)
#  (2019-01-03, 2019-01-06)
#  (2019-01-03, 2019-01-10)
#  (2019-01-06, 2019-01-10)

@test span(geolist) == 9
@test span(intlist) == 9