text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler multicomponent equations

# 1) Dry Air  2) Helium + 28% Air
equations = CompressibleEulerMulticomponentEquations2D(gammas        = (1.4, 1.648),
                                                       gas_constants = (0.287, 1.578))

"""
    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})

A shock-bubble testcase for multicomponent Euler equations
- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
"""
function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})
  # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
  # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares
  @unpack gas_constants = equations

  # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
  delta   = 0.03

  # Region I
  rho1_1  = delta
  rho2_1  = 1.225 * gas_constants[1]/gas_constants[2] - delta
  v1_1    = zero(delta)
  v2_1    = zero(delta)
  p_1     = 101325

  # Region II
  rho1_2  = 1.225-delta
  rho2_2  = delta
  v1_2    = zero(delta)
  v2_2    = zero(delta)
  p_2     = 101325

  # Region III
  rho1_3  = 1.6861 - delta
  rho2_3  = delta
  v1_3    = -113.5243
  v2_3    = zero(delta)
  p_3     = 159060

  # Set up Region I & II:
  inicenter = SVector(zero(delta), zero(delta))
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)

  if (x[1] > 0.50)
    # Set up Region III
    rho1    = rho1_3
    rho2    = rho2_3
    v1      = v1_3
    v2      = v2_3
    p       = p_3
  elseif (r < 0.25)
    # Set up Region I
    rho1    = rho1_1
    rho2    = rho2_1
    v1      = v1_1
    v2      = v2_1
    p       = p_1
  else
    # Set up Region II
    rho1    = rho1_2
    rho2    = rho2_2
    v1      = v1_2
    v2      = v2_2
    p       = p_2
  end

  return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
end
initial_condition = initial_condition_shock_bubble

surface_flux        = flux_lax_friedrichs
volume_flux         = flux_ranocha
basis               = LobattoLegendreBasis(3)
indicator_sc        = IndicatorHennemannGassner(equations, basis,
                                                alpha_max=0.5,
                                                alpha_min=0.001,
                                                alpha_smooth=true,
                                                variable=density_pressure)
volume_integral     = VolumeIntegralShockCapturingHG(indicator_sc;
                                                     volume_flux_dg=volume_flux,
                                                     volume_flux_fv=surface_flux)
solver              = DGSEM(basis, surface_flux, volume_integral)

coordinates_min     = (-2.25, -2.225)
coordinates_max     = ( 2.20,  2.225)
mesh                = TreeMesh(coordinates_min, coordinates_max,
                               initial_refinement_level=3,
                               n_cells_max=1_000_000)

semi                = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan               = (0.0, 0.01)
ode                 = semidiscretize(semi, tspan)

summary_callback    = SummaryCallback()

analysis_interval   = 300
analysis_callback   = AnalysisCallback(semi, interval=analysis_interval,
                                       extra_analysis_integrals=(Trixi.density,))

alive_callback      = AliveCallback(analysis_interval=analysis_interval)

save_solution       = SaveSolutionCallback(interval=300,
                                           save_initial_solution=true,
                                           save_final_solution=true,
                                           solution_variables=cons2prim)

stepsize_callback   = StepsizeCallback(cfl=0.3)

callbacks           = CallbackSet(summary_callback,
                                  analysis_callback,
                                  alive_callback,
                                  save_solution,
                                  stepsize_callback)

###############################################################################
# run the simulation
sol                 = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
                                                      dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                                                      save_everystep=false,
                                                      callback=callbacks,
                                                      maxiters=1e5);
summary_callback() # print the timer summary


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

#
# The Computer Language Benchmarks Game
# nbody benchmark
# http://shootout.alioth.debian.org/u32/performance.php?test=nbody
#
# A straight port from the Java version
#

module NBody

# Constants
const solar_mass = 4 * pi * pi
const days_per_year = 365.24

# A heavenly body in the system
mutable struct Body
    x::Float64
    y::Float64
    z::Float64
    vx::Float64
    vy::Float64
    vz::Float64
    mass::Float64
end

function offset_momentum(b::Body, px, py, pz)
    b.vx = -px / solar_mass
    b.vy = -py / solar_mass
    b.vz = -pz / solar_mass
end

function init_sun(bodies)
    local px::Float64 = 0.0
    local py::Float64 = 0.0
    local pz::Float64 = 0.0
    for b in bodies
        px += b.vx * b.mass
        py += b.vy * b.mass
        pz += b.vz * b.mass
    end
    offset_momentum(bodies[1], px, py, pz)
end

function advance(bodies, dt)
    for i = 1:length(bodies)
        for j = i+1:length(bodies)
            dx = bodies[i].x - bodies[j].x
            dy = bodies[i].y - bodies[j].y
            dz = bodies[i].z - bodies[j].z
            dsq = dx^2 + dy^2 + dz^2
            distance = sqrt(dsq)
            mag = dt / (dsq * distance)

            bodies[i].vx -= dx * bodies[j].mass * mag
            bodies[i].vy -= dy * bodies[j].mass * mag
            bodies[i].vz -= dz * bodies[j].mass * mag

            bodies[j].vx += dx * bodies[i].mass * mag
            bodies[j].vy += dy * bodies[i].mass * mag
            bodies[j].vz += dz * bodies[i].mass * mag
        end
    end

    for b in bodies
        b.x += dt * b.vx
        b.y += dt * b.vy
        b.z += dt * b.vz
    end
end

function energy(bodies)
    local e::Float64 = 0.0
    for i = 1:length(bodies)
        e += 0.5 * bodies[i].mass *
             (bodies[i].vx^2 + bodies[i].vy^2 + bodies[i].vz^2)
        for j = i+1:length(bodies)
            dx = bodies[i].x - bodies[j].x
            dy = bodies[i].y - bodies[j].y
            dz = bodies[i].z - bodies[j].z
            distance = sqrt(dx^2 + dy^2 + dz^2)
            e -= (bodies[i].mass * bodies[j].mass) / distance
        end
    end
    e
end

function nbody(N::Int=1000)
    jupiter = Body( 4.84143144246472090e+00,                   # x
                   -1.16032004402742839e+00,                   # y
                   -1.03622044471123109e-01,                   # z
                   1.66007664274403694e-03 * days_per_year,   # vx
                   7.69901118419740425e-03 * days_per_year,   # vy
                   -6.90460016972063023e-05 * days_per_year,   # vz
                   9.54791938424326609e-04 * solar_mass)      # mass

    saturn = Body( 8.34336671824457987e+00,
                  4.12479856412430479e+00,
                  -4.03523417114321381e-01,
                  -2.76742510726862411e-03 * days_per_year,
                  4.99852801234917238e-03 * days_per_year,
                  2.30417297573763929e-05 * days_per_year,
                  2.85885980666130812e-04 * solar_mass)

    uranus = Body( 1.28943695621391310e+01,
                  -1.51111514016986312e+01,
                  -2.23307578892655734e-01,
                  2.96460137564761618e-03 * days_per_year,
                  2.37847173959480950e-03 * days_per_year,
                  -2.96589568540237556e-05 * days_per_year,
                  4.36624404335156298e-05 * solar_mass)

    neptune = Body( 1.53796971148509165e+01,
                   -2.59193146099879641e+01,
                   1.79258772950371181e-01,
                   2.68067772490389322e-03 * days_per_year,
                   1.62824170038242295e-03 * days_per_year,
                   -9.51592254519715870e-05 * days_per_year,
                   5.15138902046611451e-05 * solar_mass)

    sun = Body(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, solar_mass)

    bodies = [sun, jupiter, saturn, uranus, neptune]

    init_sun(bodies)
#    @printf("%.9f\n", energy(bodies))
    for i = 1:N
        advance(bodies, 0.01)
    end
#    @printf("%.9f\n", energy(bodies))
end

end # module

ext: jl

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
"""
    cachedir([file=Base.PROGRAM_FILE])

Return the cache directory.
"""
function cachedir(file = Base.PROGRAM_FILE)
    name, _ = splitext(basename(file))
    dir = joinpath(dirname(file), "." * name * ".poptart")
    isabspath(file) || return joinpath(pwd(), dir)
    return dir
end

"""
    cache_file([file=Base.PROGRAM_FILE])

Return the files that will be cached:

1. `cmd.jl`: a julia script that contains the generated CLI code.
2. `checksum`: a checksum file for checking if the generated CLI code matches the original file.
"""
function cachefile(file = Base.PROGRAM_FILE)
    dir = cachedir(file)
    return joinpath(dir, "poptart_gui.jl"), joinpath(dir, "checksum")
end

"""
    create_cache(cmd[, file=Base.PROGRAM_FILE])

Create cache for given command `cmd` at file location `file`.
"""
function create_cache(cmd, file = Base.PROGRAM_FILE)
    isempty(file) && return
    dir = cachedir(file)
    if !ispath(dir)
        mkpath(dir)
    end

    cache_file, crc = cachefile(file)
    open(cache_file, "w") do io
        gui_write(io, cmd)
    end
    write(crc, string(checksum(file), base = 16))
    return
end

"""
    write([io], cmd::EntryCommand)

Write the generated CLI script into a Julia script file. Default is the [`cachefile`](@ref).
"""
function gui_write(io::IO, x::EntryCommand)
    println(io, "#= generated by ComoniconGUI =#")
    main_code = prettify(codegen(x))
    if main_code.head == :block
        for arg in main_code.args
            println(io, arg)
        end
    else
        println(io, main_code)
    end
end


text:
"""
    dijkstra(graph::Graph, source::Int, target::Int)::Tuple{AbstractFloat, Vector{Int}}

The shortest path distance from `source` to `target` using `Dijkstra's` algorithm
Also returns the shortest path
"""
function dijkstra(graph::Graph, source::Int, target::Int = -1)::Tuple{AbstractFloat, Vector{Int}}
    if target > 0
        condition = x -> target ∈ x
    else
        condition = x -> !isempty(x)
    end
    dist, pred = _dijkstra_core(graph, source, condition)
    return dist[target], path(pred, target, source)
end

"""
    dijkstra(graph::Graph, source::Int)::Tuple{Vector{AbstractFloat}, Vector{Int}}

The shortest path distances from `source` to all nodes using `Dijkstra's` algorithm
Also returns the list of predecessors
"""
function dijkstra(graph::Graph, source::Int)::Tuple{Vector{AbstractFloat}, Vector{Int}}
    condition = x -> !isempty(x)
    dist, pred =  _dijkstra_core(graph, source, condition)
    return dist, pred
end

"""
    _dijkstra_core(graph::Graph, source::Int, condition::Function)::Tuple{Vector{AbstractFloat}, Vector{Int}}

The function that contains the actual implementation of the `Dijkstra's` algorithm
Please run the `dijkstra` function and not this
"""
function _dijkstra_core(graph::Graph, source::Int, condition::Function)::Tuple{Vector{AbstractFloat}, Vector{Int}}
    node_size = size(graph, 1)
    dist = ones(Float64, node_size) * Inf
    dist[source] = 0.0
    Q = Set{Int64}()  # visited nodes
    T = Set{Int64}(1:node_size)  # unvisited nodes
    pred = ones(Int64, node_size) * -1
    while condition(T)
        # node selection
        untraversed_nodes = [(d, k) for (k, d) in enumerate(dist) if k in T]
        if minimum(untraversed_nodes)[1] == Inf
            break # Break if remaining nodes are disconnected
        end
        node_ind = untraversed_nodes[argmin(untraversed_nodes)][2]
        push!(Q, node_ind)
        delete!(T, node_ind)
        # distance update
        curr_node = graph.nodes[node_ind]
        for (neigh, edge) in adj(graph, curr_node)
            t_ind = neigh.index
            weight = edge.cost
            if dist[t_ind] > dist[node_ind] + weight
                dist[t_ind] = dist[node_ind] + weight
                pred[t_ind] = node_ind
            end
        end
    end
    return dist, pred
end


text:
using GeneralUtils
using MAT
using JLD

include("constrained_parabolic_minimization.jl")

"""
pdict = wallwrap(bdict, pdict)
Given bdict, a dictionary of symbols to [minval, maxval] vectors, and pdict, a dictionary of symbols
to values (or, alternatively, an Array of (Symbol, value) tuples], goes through each of the symbols in
bdict and modifies the corresponding value in pdict putting it through a tanh so the final output lies
within the limits in bdict.  Returns the new pdict.  Makes a copy of pdict so as not to modify the original.
"""
function wallwrap(bdict, epdict)
    local pdict = two_level_copy(epdict)  # Must be very careful here! I got bit by the bug of forgetting that without
    # an explicit copy() call, Julia does not make copies of the contents of arrays or dictionaries, making it
    # easy to inadvertently modify something one did not intend to perturb.  Note the two_level_copy() call,
    # necessary to make sure we don't mess up the content of the caller's dictionary.

    if typeof(pdict)<:Array
        pdict = Dict(pdict)
    end

    allkeys = keys(bdict)

    for k in allkeys
        local bbox = bdict[k]
        d = 0.5*(bbox[2] - bbox[1])
        m = 0.5*(bbox[2] + bbox[1])

        pdict[k] = bbox[1] + d*(tanh((pdict[k]-m)/d)+1)
    end
    return pdict
end

"""
params = vector_wrap(bbox, args, eparams)
Given bdict, a dictionary of symbols to [minval, maxval] vectors, args, an array of strings representing
symbols, and params, an array of values corresponding to the args list, puts each param that has an entry
in bdict through the tanh-walling mechanism, and returns the result. Does not modify the contents of the
original params vector (or bdict or args).
"""
function vector_wrap(bbox, args, eparams)
    local params = two_level_copy(eparams)
    pdict = wallwrap(bbox, make_dict(args, params))
    i=1; j=1
    for i=1:length(args)
        if typeof(args[i])<:Array
            params[j:j+args[i][2]-1] = pdict[Symbol(args[i][1])]
            j += args[i][2]-1
        else
            params[j] = pdict[Symbol(args[i])]
        end
    j = j+1
    end
    return params
end

"""
params = inverse_wall(bdict, args, wparams)
Given bdict, a dictionary of symbols to [minval, maxval] vectors, args, an array of strings representing
symbols, and wparams, an array of values corresponding to the args list where each param that has an entry
in bdict has alreadt been through the tanh-walling mechanism, UNwalls the ones that have a bdict entry and
returns the result. Does not modify the contents of the original params vector (or bdict or args).
"""
function inverse_wall(bdict, args, wparams)
    local params = two_level_copy(wparams)
    pdict = inverse_wall(bdict, make_dict(args, params))
    i=1; j=1
    for i=1:length(args)
        if typeof(args[i])<:Array
            params[j:j+args[i][2]-1] = pdict[Symbol(args[i][1])]
            j += args[i][2]-1
        else
            params[j] = pdict[Symbol(args[i])]
        end
        j = j+1
    end
    return params
end

"""
pdict = inverse_wall(bdict, wdict)
Given bdict, a dictionary of symbols to [minval, maxval] vectors, and wdict, a dictionary of symbols to values
(or vectors of values)  UNwalls the ones that have a bdict entry and
returns the result. Does not modify the contents of any dictionaries.
"""
function inverse_wall(bdict, wdict)
    local pdict = two_level_copy(wdict)

    allkeys = keys(bdict)
    for k in allkeys
        local bbox = bdict[k]
        d = 0.5*(bbox[2] - bbox[1])
        m = 0.5*(bbox[2] + bbox[1])

        pdict[k] = m + d*0.5*log((pdict[k]-bbox[1])./(2*d - pdict[k] + bbox[1]))
    end
    return(pdict)
end

######################################################
#                                                    #
#         BBOX_HESSIAN_KEYWORD_MINIMIZATION          #
#                                                    #
######################################################

function bbox_Hessian_keyword_minimization(seed, args, bbox, func; start_eta=0.1, tol=1e-6, maxiter=400,
    verbose=false, verbose_level=1, verbose_every=1, softbox=true, hardbox=false, report_file="")

    # --- check that saving will be done to a .jld file ---
    if length(report_file)>0 && splitext(report_file)[2] != ".jld"
        if splitext(report_file)[2] == ""
            report_file = resport_file * ".jld"
        else
            error("Sorry, report_file can only write to JLD files, the extension has to be .jld")
        end
    end

    # --------- Initializing the trajectory trace and function wrapper--------

    traj_increment = 100
    params = 0  # Make sure to have this here so that params stays defined beyond the try/catch
    if ( !(typeof(bbox)<:Dict) ); error("Currently only supporting softbox=true, bbox must be a Dict"); end;
    try
        params = copy(inverse_wall(bbox, args, seed))
    catch
        error("Were all initial param values within the indicated walls?")
    end
    eta = start_eta
    trajectory = zeros(2+length(params), traj_increment); cpm_traj = zeros(2, traj_increment)

    ftraj = Array{Any}(3,0)  # will hold gradient, hessian, and further_out,  per iteration

    further_out =[];  # We define this variable here so it will be available for stashing further outputs from func

    # Now we define a wrapper around func() to do three things: (a) wallwrap parameters using the softbox method;
    # (b) return as the desired scalar the first output of func; (c) stash in further_out any further outputs of func
    internal_func = (;pars...) -> begin
        fresults = func(;wallwrap(bbox, pars)...)   # note use of bbox external to this begin...end
        if typeof(fresults)<:Tuple
            answer = fresults[1]
            further_out = fresults[2:end]
        else
            answer = fresults
        end
        return answer  # we assume that the first output of func() will always be a scalar, and that's what we return for ForwardDiff
    end

    # --------- END Initializing the trajectory trace --------
    println("-------- end of Initialization --------")

    if verbose
        @printf "%d: eta=%g ps=" 0 eta
        print_vector(vector_wrap(bbox, args, params))
        @printf "\n"
    end

    if softbox
        if !(typeof(bbox)<:Dict); error("bhm: If softbox=true, then bbox must eb a Dict"); end
        cost, grad, hess = keyword_vgh(internal_func, args, params)  # further_out will be mutated
    elseif hardbox
        error("Sorry, no longer supporting hardbox=true")
    else
        error("Sorry, no longer supporting softbox=false")
    end

    chessdelta = zeros(size(params))

    i=0  # here so variable i is available outside the loop
    for i in [1:maxiter;]
        if i > size(trajectory, 2)
            trajectory = [trajectory zeros(2+length(params), traj_increment)]
            cpm_traj   = [cpm_traj   zeros(2, traj_increment)]
        end
        trajectory[1:2, i]   = [eta;cost]
        trajectory[3:end, i] = vector_wrap(bbox, args, params)
        ftraj = [ftraj [grad, hess, further_out]]

        if length(report_file)>0
            save(report_file, Dict("traj"=>trajectory[:,1:i], "cpm_traj"=>cpm_traj[:,1:i], "ftraj"=>ftraj))
        end

        # println(hess);
        hessdelta  = - pinv(hess)*grad
        try
            if verbose && verbose_level >= 2
                @printf("bhm: about to try cpm with grad : "); print_vector_g(grad); print("\n")
                @printf("bhm:   hess :"); print_vector_g(hess[:]); print("\n");
            end
            if verbose && verbose_level >= 2
                cpm_out = constrained_parabolic_minimization(hess, grad'', eta,
                    maxiter=500, tol=1e-20, do_plot=true, verbose=true)
            else
                cpm_out = constrained_parabolic_minimization(hess, grad'', eta, maxiter=500, tol=1e-20)
            end
            chessdelta = cpm_out[1]; cpm_traj[1,i] = cpm_out[5]; cpm_traj[2,i] = cpm_out[6]
            jumptype = "not failed"
        catch y
            jumptype = "failed"
            if verbose
                @printf "Constrained parabolic minimization failed with error %s\n" y
                @printf "\n"
                @printf "eta was %g\n" eta
                @printf "grad was\n"
                print_vector(grad)
                @printf "\n\nhess was\n"
                for k in [1:length(grad);]
                    print_vector(hess[k,:])
                    @printf "\n"
                end
                @printf "\n"
                matwrite("error_report.mat", Dict("grad"=>grad, "hess"=>hess, "eta"=>eta))
            end
            break
        end

        if norm(hessdelta) <= eta
            new_params = params + hessdelta
            jumptype = "Newton"
        elseif jumptype != "failed"
            new_params = params + chessdelta
            jumptype  = "constrained"
        end

        if jumptype != "failed"
            new_cost, new_grad, new_hess = keyword_vgh(internal_func, args, new_params)   # further_out may mutate
            if verbose && verbose_level >=2
                @printf("bhm: had new_params = : "); print_vector_g(vector_wrap(bbox, args, params)); print("\n");
                @printf("bhm: and my bbox was : "); print(bbox); print("\n")
                @printf("bhm: and my wallwrap output was : "); print(wallwrap(bbox, make_dict(args, new_params))); print("\n")
                @printf("bhm: and this produced new_grad : "); print_vector_g(new_grad); print("\n")
                @printf("bhm:   new_hess :"); print_vector_g(new_hess[:]); print("\n");
            end

            if abs(new_cost - cost) < tol || eta < tol
                if verbose
                    @printf("About to break -- tol=%g, new_cost-cost=%g, eta=%g\n", tol, new_cost-cost, eta)
                end
                break
            end
        end

        if jumptype == "failed" || new_cost >= cost
            if verbose
                @printf("eta going down: new_cost-cost=%g and jumptype='%s'\n", new_cost-cost, jumptype)
                if verbose_level >= 2
                    nwp = vector_wrap(bbox, args, new_params); wp = vector_wrap(bbox, args, params)
                    @printf("   vvv: proposed new params were : "); print_vector_g(nwp); print("\n")
                    @printf("   vvv: proposed delta params was : "); print_vector_g(nwp-wp); print("\n")
                    @printf("   vvv: grad was : "); print_vector_g(grad); print("\n")
                    costheta = dot(new_params-params, grad)/(norm(new_params-params)*norm(grad))
                    @printf("   vvv: costheta of proposed jump was %g\n", costheta)
                end
            end
            eta = eta/2
            costheta = NaN
            if eta < tol
                if verbose
                    @printf("About to break -- tol=%g, new_cost-cost=%g, eta=%g\n", tol, new_cost-cost, eta)
                end
                break
            end
        else
            eta = eta*1.1
            costheta = dot(new_params-params, grad)/(norm(new_params-params)*norm(grad))

            params = new_params
            cost = new_cost
            grad = new_grad
            hess = new_hess
        end

        if verbose
            if rem(i, verbose_every)==0
                @printf "%d: eta=%g cost=%g jtype=%s costheta=%.3f ps=" i eta cost jumptype costheta
                print_vector_g(vector_wrap(bbox, args, params))
                @printf "\n"
                if verbose_level >= 3
                    @printf "    At this point, grad is ="
                    print_vector_g(grad)
                    @printf "\n"
                end
            end
        end
    end

    trajectory = trajectory[:,1:i]; cpm_traj = cpm_traj[:,1:i]
    if length(report_file)>0
        save(report_file, Dict("traj"=>trajectory, "cpm_traj"=>cpm_traj, "ftraj"=>ftraj))
    end

    return vector_wrap(bbox, args, params), trajectory, cost, cpm_traj, ftraj, hess
end

function error_ellipse(x_bf, hessian, idx)

# x_bf : best-fit parameter for the offset of ellipse
# hessian : Hessian matrix for covariance matrix
# idx : id of parameters to get the covrariance error ellipse

    covariance = inv(hessian[idx,idx]);

    eigenval, eigenvec = eig(covariance);
    println(eigenvec)

    # Get the largest eigenvalue
    # Get the index of the largest eigenvector
    largest_eigenval, largest_eigenvec_ind_c = findmax(eigenval);
    largest_eigenvec = eigenvec[:, largest_eigenvec_ind_c];

    # Get the smallest eigenvector and eigenvalue
    if largest_eigenvec_ind_c == 1
        smallest_eigenval = eigenval[2]
        smallest_eigenvec = eigenvec[:,2];
    else
        smallest_eigenval = eigenval[1]
        smallest_eigenvec = eigenvec[1,:];
    end

    # Calculate the angle between the x-axis and the largest eigenvector
    angle = atan2(largest_eigenvec[2], largest_eigenvec[1]);

    # This angle is between -pi and pi.
    # Let's shift it such that the angle is between 0 and 2pi
    if angle < 0
        angle = angle + 2*pi;
    end

    # % Get the 95% confidence interval error ellipse
    chisquare_val_95 = 2.1459;

    theta_grid = linspace(0,2*pi,200);
    phi = angle;

    # % x0,y0 ellipse centre coordinates
    X0=x_bf[idx[1]];
    Y0=x_bf[idx[2]];
    a=sqrt(largest_eigenval);
    b=sqrt(smallest_eigenval);

    # % the ellipse in x and y coordinates
    ellipse_x_r  = chisquare_val_95*a*cos.( theta_grid );
    ellipse_y_r  = chisquare_val_95*b*sin.( theta_grid );

    # %Define a rotation matrix
    R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];

    # let's rotate the ellipse to some angle phi
    r_ellipse = [ellipse_x_r ellipse_y_r] * R;

    Xs = r_ellipse[:,1] + X0
    Ys = r_ellipse[:,2] + Y0
    return Xs,Ys
end

function error_ellipse_cov(x_bf, covariance, idx)

# x_bf : best-fit parameter for the offset of ellipse
# hessian : Hessian matrix for covariance matrix
# idx : id of parameters to get the covrariance error ellipse

    # covariance = inv(hessian[idx,idx]);

    eigenval, eigenvec = eig(covariance);

    # Get the largest eigenvalue
    # Get the index of the largest eigenvector
    largest_eigenval, largest_eigenvec_ind_c = findmax(eigenval);
    largest_eigenvec = eigenvec[:, largest_eigenvec_ind_c];

    # Get the smallest eigenvector and eigenvalue
    if largest_eigenvec_ind_c == 1
        smallest_eigenval = eigenval[2]
        smallest_eigenvec = eigenvec[:,2];
    else
        smallest_eigenval = eigenval[1]
        smallest_eigenvec = eigenvec[1,:];
    end

    # Calculate the angle between the x-axis and the largest eigenvector
    angle = atan2(largest_eigenvec[2], largest_eigenvec[1]);

    # This angle is between -pi and pi.
    # Let's shift it such that the angle is between 0 and 2pi
    if angle < 0
        angle = angle + 2*pi;
    end

    # % Get the 95% confidence interval error ellipse
    chisquare_val_95 = 2.1459;

    theta_grid = linspace(0,2*pi,200);
    phi = angle;

    # % x0,y0 ellipse centre coordinates
    X0=x_bf[idx[1]];
    Y0=x_bf[idx[2]];
    a=sqrt(largest_eigenval);
    b=sqrt(smallest_eigenval);

    # % the ellipse in x and y coordinates
    ellipse_x_r  = chisquare_val_95*a*cos.( theta_grid );
    ellipse_y_r  = chisquare_val_95*b*sin.( theta_grid );

    # %Define a rotation matrix
    R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];

    # let's rotate the ellipse to some angle phi
    r_ellipse = [ellipse_x_r ellipse_y_r] * R;

    Xs = r_ellipse[:,1] + X0
    Ys = r_ellipse[:,2] + Y0
    return Xs,Ys
end


text:
"""
    mutable struct Affine <: Block

Applies a linear transformation to the incoming data: y = w * x
"""
mutable struct Affine <: Block
    w::VarOrNil # input to hidden weights
    function Affine(inputSize::Int, hiddenSize::Int; type::Type=Array{Float32})
        T = eltype(type)
        A = sqrt(T(1/hiddenSize))
        w = randn(T, hiddenSize, inputSize) .* A
        new(Variable{type}(w,true,true,true))
    end
    function Affine()
        new(nothing)
    end
end

function clone(this::Affine; type::Type=Array{Float32})
    cloned = Affine()
    cloned.w = clone(this.w, type=type)
    return cloned
end

# pretty show
function Base.show(io::IO, m::Affine)
    SIZE = size(m.w)
    TYPE = typeof(m.w.value)
    print(io, "Affine($(SIZE[2]), $(SIZE[1]); type=$TYPE)")
end

"""
    unbiasedof(m::Affine)

unbiased weights of Affine block
"""
function unbiasedof(m::Affine)
    weights = Vector(undef, 1)
    weights[1] = m.w.value
    return weights
end

function weightsof(m::Affine)
    weights = Vector(undef, 1)
    weights[1] = m.w.value
    return weights
end

function gradsof(m::Affine)
    grads = Vector(undef, 1)
    grads[1] = m.w.delta
    return grads
end

function zerograds!(m::Affine)
    for v in gradsof(m)
        v .= 0.0
    end
end

function paramsof(m::Affine)
    params = Vector{Variable}(undef,1)
    params[1] = m.w
    return params
end

function xparamsof(m::Affine)
    xparams = Vector{XVariable}(undef,1)
    xparams[1] = ('w', m.w)
    return xparams
end

function nparamsof(m::Affine)
    return length(m.w)
end

function bytesof(model::Affine, unit::String="MB")
    n = nparamsof(model) * elsizeof(model.w)
    return blocksize(n, uppercase(unit))
end

function forward(m::Affine, x::Variable)
    w = m.w
    return w * x
end

function predict(m::Affine, x)
    w = m.w.value
    return w * x
end

function to(type::Type, m::Affine)
    m.w = to(type, m.w)
    return m
end

function to!(type::Type, m::Affine)
    m.w = to(type, m.w)
    return nothing
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
"""
Locator for a single data entry, uniquely defined by parameter names partitioned into file system path and group path.

**Fields**
- `root` -- the root directory of the data entry.
- `folder_entries`  -- the ordered partition of parameters names which define the file system path of the data entry.
- `group_entries`   -- the ordered partition of parameters names which define the group path of the data entry.
"""
struct DataEntry
    root::String
    folder_entries::Array{Array{String,1},1}
    group_entries::Array{Array{String,1},1}
    DataEntry(r, f, g) = new(r, [sort(i) for i in f], [sort(i) for i in g])
end

"""
    print_params(format_dict, value_dict, key)

Print out the path of parameter set given by `key`. `format_dict` and `value_dict` are look-up tables for the parameter formatters and values respectively. `key` should be an valid partition of parameter names.
"""
function print_params(format_dict, value_dict, key::Array{String,1})
    res = ""
    for name in key
        if haskey(value_dict, name)
            value = value_dict[name]
        elseif haskey(value_dict, Symbol(name))
            value = value_dict[Symbol(name)]
        else
            throw(KeyError(name))
        end
        data_string = cfmt(format_dict[name], value)
        res = res * name * "=" * data_string * "_"
    end
    try
        chop(res, tail=1)
    catch
        ""
    end
end

function print_params(format_dict, value_dict, key::Array{Array{String,1},1})
    res = ""
    for i in key
        param_str = print_params(format_dict, value_dict, i)
        res = res * param_str * "/"
    end
    try
        chop(res, tail=1)
    catch
        ""
    end
end

"""
    get_folder_path(d::DataEntry, values)

Print out the file system path given data entry `d` and the corresponding `values`.
"""
function get_folder_path(d::DataEntry, values::Dict)
    joinpath(d.root, print_params(get_param_set(), values, d.folder_entries))
end

function get_folder_path(d::DataEntry, values)
    v = Dict(values...)
    joinpath(d.root, print_params(get_param_set(), v, d.folder_entries))
end

"""
    get_group_path(d::DataEntry, values)

Print out the group path given data entry `d` and the corresponding `values`
"""
function get_group_path(d::DataEntry, values::Dict)
    print_params(get_param_set(), values, d.group_entries)
end

function get_group_path(d::DataEntry, values)
    v = Dict(values...)
    print_params(get_param_set(), v, d.group_entries)
end


text:
using DataFrames
using CSV
using Gurobi
using JuMP

distance_data = CSV.read("distance_data.csv",header=[:i, :j, :dist, :route], missing=:none)

N = length(unique(distance_data.i))
const distance_mx = zeros(Float64,(N,N))
for r in 1:nrow(distance_data)
   distance_mx[distance_data.i[r],distance_data.j[r]] = distance_data.dist[r]
end

#http://opensourc.es/blog/mip-tsp
N = 50 # this time we use 50 points

m = Model(solver=GurobiSolver());
@variable(m, x[f=1:N,t=1:N], Bin)
@objective(m, Min, sum( x[i,j]*distance_mx[i,j] for i=1:N,j=1:N))
@constraint(m, notself[i=1:N], x[i,i] == 0)
@constraint(m, oneout[i=1:N], sum(x[i,1:N]) == 1)
@constraint(m, onein[j=1:N], sum(x[1:N,j]) == 1)
for f=1:N, t=1:N
    @constraint(m, x[f,t]+x[t,f] <= 1)
end

function getcycle(m, N)
    x_val = getvalue(x)
    cycle_idx = Vector{Int}()
    push!(cycle_idx, 1)
    while true
        v, idx = findmax(x_val[cycle_idx[end], 1:N])
        if idx == cycle_idx[1]
            break
        else
            push!(cycle_idx, idx)
        end
    end
    cycle_idx
end

function callbackhandle(cb)
    cycle_idx = getcycle(cb, N)
    println("Callback! N= $N cycle_idx: ", cycle_idx)
    println("Length: ", length(cycle_idx))
    if length(cycle_idx) < N
        @lazyconstraint(cb, sum(x[cycle_idx,cycle_idx]) <= length(cycle_idx)-1)
        println("added a lazy constraint")
    end
end

addlazycallback(m, callbackhandle)
solve(m)

sbws_la = CSV.read("Subway_LV.csv", missing=:none)

using PyPlot
cycle_idx = getcycle(m, N)
ids = vcat(cycle_idx, cycle_idx[1])
PyPlot.cla()
PyPlot.ioff()
p = PyPlot.plot(sbws_la.long[ids], sbws_la.latt[ids])
PyPlot.display_figs()
PyPlot.show()


text:
@testset "constraints.jl" begin
    @testset "MOI.add_constraint(::NotEqualSet)" begin
        opt = SeaPearl.Optimizer()

        x = MOI.add_variable(opt)
        y = MOI.add_variable(opt)

        ci = MOI.add_constraint(opt, MOI.VectorOfVariables([x, y]), SeaPearl.NotEqualSet())
        @test ci == MOI.ConstraintIndex{MOI.VectorOfVariables, SeaPearl.NotEqualSet}(1)
        @test length(opt.moimodel.constraints) == 1
        @test opt.moimodel.constraints[1] == SeaPearl.MOIConstraint(SeaPearl.NotEqualSet, (x, y), ci)
    end
    @testset "create_CPConstraint(::NotEqualSet)" begin
        opt = SeaPearl.Optimizer()

        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(2)))
        x = MOI.VariableIndex(1)
        y = MOI.VariableIndex(2)

        ci = MOI.add_constraint(opt, MOI.VectorOfVariables([x, y]), SeaPearl.NotEqualSet())

        SeaPearl.bridge_variables!(opt)

        cp_constraint = SeaPearl.create_CPConstraint(opt.moimodel.constraints[1], opt)
        @test isa(cp_constraint, SeaPearl.NotEqual)
        @test SeaPearl.variablesArray(cp_constraint) == [SeaPearl.get_cp_variable(opt, x), SeaPearl.get_cp_variable(opt, y)]
    end
    @testset "MOI.add_constraint(::LessThan)" begin
        opt = SeaPearl.Optimizer()

        x = MOI.add_variable(opt)

        set = MOI.LessThan(5.)
        ci = MOI.add_constraint(opt, x, set)
        @test opt.moimodel.variables[1].max == 5
        @test isnothing(opt.moimodel.variables[1].min)
    end
    @testset "MOI.add_constraint(::GreaterThan)" begin
        opt = SeaPearl.Optimizer()

        x = MOI.add_variable(opt)

        set = MOI.GreaterThan(5.)
        ci = MOI.add_constraint(opt, x, set)
        @test opt.moimodel.variables[1].min == 5
        @test isnothing(opt.moimodel.variables[1].max)
    end
    @testset "MOI.add_constraint(::MOI.ScalarAffineFunction, ::MOI.EqualTo)" begin
        opt = SeaPearl.Optimizer()

        x = MOI.add_variable(opt)
        y = MOI.add_variable(opt)
        aff = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5., -3.], [x, y]), 3.)

        ci = MOI.add_constraint(opt, aff, MOI.EqualTo(5.))
        @test aff.constant == -2
        @test length(opt.moimodel.constraints) == 1
        @test opt.moimodel.constraints[1] == SeaPearl.MOIConstraint(MOI.EqualTo, (SeaPearl.AffineIndex(1),), ci)
        @test length(opt.moimodel.affines) == 1
        @test opt.moimodel.affines[1].content == aff
        @test isnothing(opt.moimodel.affines[1].cp_identifier)
    end
    @testset "create_CPConstraint(::MOI.EqualTo)" begin
        opt = SeaPearl.Optimizer()

        x = MOI.add_variable(opt)
        MOI.add_constraint(opt, x, MOI.LessThan(2.))
        MOI.add_constraint(opt, x, MOI.GreaterThan(1.))
        y = MOI.add_variable(opt)
        MOI.add_constraint(opt, y, MOI.LessThan(5.))
        MOI.add_constraint(opt, y, MOI.GreaterThan(0.))
        aff = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5., -3.], [x, y]), 3.)

        ci = MOI.add_constraint(opt, aff, MOI.EqualTo(5.))

        SeaPearl.bridge_variables!(opt)
        SeaPearl.bridge_affines!(opt)

        cp_constraint = SeaPearl.create_CPConstraint(opt.moimodel.constraints[1], opt)

        @test SeaPearl.assignedValue(opt.cpmodel.variables["7"]) == 0

        # Sometimes clearly infeasible
        opt = SeaPearl.Optimizer()

        x = MOI.add_variable(opt)
        MOI.add_constraint(opt, x, MOI.LessThan(2.))
        MOI.add_constraint(opt, x, MOI.GreaterThan(1.))
        y = MOI.add_variable(opt)
        MOI.add_constraint(opt, y, MOI.LessThan(5.))
        MOI.add_constraint(opt, y, MOI.GreaterThan(4.))
        aff = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5., -3.], [x, y]), 3.)

        ci = MOI.add_constraint(opt, aff, MOI.EqualTo(5.))

        SeaPearl.bridge_variables!(opt)
        SeaPearl.bridge_affines!(opt)

        @test_throws AssertionError SeaPearl.create_CPConstraint(opt.moimodel.constraints[1], opt)
    end
end


text:
using Navigation
using Test

@testset "Navigation.jl" begin
    tests = ["utility", "distances", "bearings", "points", "wind", "airspace"]
    # tests = ["airspace"]

    for t in tests
        include("$(t).jl")
    end
end


text:
function Mstep(d::LogNormal, y)
    mu = mean(log, y)
    sigma = sqrt(mean(x -> (log(x)-mu)^2, y))
    return LogNormal(mu, sigma)
end

function Mstep(d::GeneralizedGamma, y, lr=1)
    shp, scl, pwr = params(d)
    v = [log(shp), log(pwr)]
    g = ForwardDiff.gradient(v -> sum(y-> -logpdf(GeneralizedGamma(exp(v[1]),scl,exp(v[2])),y)), v)
    H = ForwardDiff.hessian(x -> sum(y-> -logpdf(GeneralizedGamma(exp(v[1]),scl,exp(v[2])),y)), v)
    v -= lr*(H \ g)
    b = mean(x -> x^exp(v[2]), y)^exp(-v[2])
    return GeneralizedGamma(exp(v[1]), b, exp(v[2]))
end

function Mstep(d::Gamma, y, lr=1)
    shp, scl = params(d)
    rho = log(shp)
    meanlogy = mean(log,y)
    #f(x,shp) = -shp*log(scl) + (shp-1)*log(x) - (x/scl) -loggamma(shp)
    g = shp*log(scl) - shp*meanlogy + digamma(shp)*shp
    H = shp*log(scl) - shp*meanlogy + trigamma(shp)*shp^2 + digamma(shp)*shp
    rho -= lr*(g/H)
    b = mean(y)/exp(rho)
    return Gamma(exp(rho), b)
end

function Mstep(d::Weibull, y, lr=1)
    shp, scl = params(d)
    rho = log(shp)
    logy = sum(log, y)
    n = length(y)
    A = ((y/scl).^shp)'*(log.(y/scl))
    B = ((y/scl).^shp)'*(log.(y/scl).^2)
    #f(x,shp) = log(shp)-shp*log(scl) + (shp-1)*log(x) - (x/scl)^shp
    g = n*(shp*log(scl)-1) - shp*logy + shp*A
    H = n*(shp*log(scl)-1) - shp*logy + shp*A + (shp^2)*B
    rho -= lr*(g/H)
    b = mean(x -> x^exp(rho), y)^exp(-rho)
    return Weibull(exp(rho), b)
end

# function Mstep(d::Gamma, y)
#     mlog = log(mean(y))
#     logm = mean(log, y)
#     a = shape(d)
#     ia = inv(a)
#     a = inv(ia + (-mlog + logm + log(a) - digamma(a)) / ((a^2) * (ia - trigamma(a))))
#     b = mean(y) / a
#     return Gamma(a,b)
# end

# function Mstep(d::Weibull, y)
#     a, b = params(d)
#     n = length(y)
#     powa = (y.^a)'
#     logy = log.(y)

#     fa = n*(inv(a) - inv(b)) + sum(logy)- n*log(b) - powa*(logy-log(b))/(b^a)
#     dfa = -n*inv(a^2) - (powa*(logy-log(b)).^2)/(b^a)

#     Delta = fx / dfx
#     a += Delta
#     b = mean(x -> x^a, ytilde)^inv(a)
#     return Weibull(a,b)
# end

function Mstep(d::Exponential, y)
    b = mean(y)
    return Exponential(b)
end

#######
#interval censored
function MCEMic(rng, dist, iter, EL, ER, S, lr=1)
    lp = zeros(iter)
    pars = params(dist)
    for it in 1:iter
    ytilde = [rand(rng, truncated(dist,S[i]-ER[i],S[i]-EL[i])) for i in eachindex(S)]
    dist = Mstep(dist, ytilde, lr)
    lp[it] = mean(x-> -logpdf(dist,x), ytilde)
    end
    return dist,lp
end

########
#interval censored with right truncated
function Estep_icrt(rng, dist, EL, ER, S, Tmax)
    N = length(S)
    ys = zeros(N)
    yb = []
    for i in 1:N
        ys[i] = rand(rng, truncated(dist,S[i]-ER[i],S[i]-EL[i]))
        u = rand(rng)
        E = EL[i] + (ER[i]-EL[i])*u
        q = cdf(dist, Tmax - E)
        if 0.0 < q < 1.0
            B = rand(rng, Geometric(q))
            if B > 0
                ysub = [rand(rng, truncated(dist,Tmax-E,Inf)) for i in 1:B]
                yb = [yb;ysub]
            end
        end
    end
    return [ys ; yb]
end

function MCEMicrt(rng, dist, iter, EL, ER, S, Tmax, lr=1)
    lp = zeros(iter)
    pars = params(dist)
    for it in 1:iter
    ytilde = Estep_icrt(rng, dist, EL, ER, S, Tmax)
    dist = Mstep(dist, ytilde, lr)
    lp[it] = mean(x-> -logpdf(dist,x), ytilde)
    end
    return dist,lp
end

######
#doubly interval censored
# function Estep_dic(rng, dist, EL, ER, SL, SR)
#     N = length(EL)
#     ys = zeros(N)
#     for i in 1:N
#         u = rand(rng)
#         E = EL[i] + (ER[i]-EL[i])*u
#         ys[i] = rand(rng, truncated(dist, SL[i]-E, SR[i]-E))
#     end
#     return ys
# end

function Estep_dic(rng, dist, EL, ER, SL, SR)
    N = length(EL)
    ys = zeros(N)
    for i in 1:N
        u = rand(rng)
        S = SL[i] + (SR[i]-SL[i])*u
        ys[i] = rand(rng, truncated(dist, S-ER[i], S-EL[i]))
    end
    return ys
end

function MCEMdic(rng, dist, iter, EL, ER, SL, SR, lr=1)
    lp = zeros(iter)
    pars = params(dist)
    for it in 1:iter
    ytilde = Estep_dic(rng, dist, EL, ER, SL, SR)
    dist = Mstep(dist, ytilde, lr)
    lp[it] = mean(x-> -logpdf(dist, x), ytilde)
    end
    return dist,lp
end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


text:
export @union, @default, @align, @struct

macro union(T, TT)
    TT.args[1] == :Union || throw(ArgumentError("2nd argument must be a `Union{T1,T2...}` type"))
    return quote
        typealias $(esc(T)) $(esc(TT))
        $(esc(:(FlatBuffers.typeorder))){TT}(::Type{$(esc(T))}, ::Type{TT}) = $(Dict{DataType,Int}([(eval(current_module(), typ), i - 1) for (i, typ) in enumerate(TT.args[2:end])]))[TT]
        $(esc(:(FlatBuffers.typeorder)))(::Type{$(esc(T))}, i::Integer) = $(Dict{Int,DataType}([(i - 1, eval(current_module(), typ)) for (i, typ) in enumerate(TT.args[2:end])]))[i]
    end
end

macro align(T, sz)
    return quote
        $(esc(:(FlatBuffers.alignment)))(::Type{$(esc(T))}) = $sz
    end
end

macro enumtype(T, typ)
    return quote
        $(esc(:(FlatBuffers.enumtype)))(::Type{$(esc(T))}) = $typ
    end
end

# recursively finds largest field of a struct
fbsizeof{T<:Enum}(::Type{T}) = sizeof(enumtype(T))
fbsizeof{T}(::Type{T}) = sizeof(T)

maxsizeof{T<:Enum}(::Type{T}) = sizeof(enumtype(T))
maxsizeof{T}(::Type{T}) = isbitstype(T) ? sizeof(T) : maximum(map(x->maxsizeof(x), T.types))

nextsizeof{T}(::Type{T}) = isbitstype(T) ? sizeof(T) : nextsizeof(T.types[1])

function fieldlayout(typ, exprs...)
    fields = Expr[]
    values = []
    largest_field = maximum(map(x->maxsizeof(eval(current_module(), x.args[2])), exprs))
    sz = cur_sz = 0
    x = 0
    for (i,expr) in enumerate(exprs)
        T = eval(current_module(), expr.args[2])
        if !isbitstype(T)
            exprs2 = [Expr(:(::), nm, typ) for (nm,typ) in zip(fieldnames(T),T.types)]
            fields2, values2 = fieldlayout(T, exprs2...)
            append!(fields, map(x->Expr(:(::), Symbol(string(expr.args[1],'_',T,'_',x.args[1])), x.args[2]), fields2))
            append!(values, map(x->x == 0 ? 0 : Expr(:call, :getfield, expr.args[1], QuoteNode(x)), values2))
        else
            push!(fields, expr)
            push!(values, expr.args[1])
        end
        sz += cur_sz = fbsizeof(T)
        if sz % largest_field == 0
            sz = cur_sz = 0
            continue
        end
        nextsz = i == length(exprs) ? 0 : nextsizeof(eval(current_module(), exprs[i+1].args[2]))
        if i == length(exprs) || cur_sz < nextsz || (sz + nextsz) > largest_field
            # this is the last field and we're not `sz % largest_field`
            # potential diffs = 7, 6, 5, 4, 3, 2, 1
            sym = expr.args[1]
            diff = cur_sz < nextsz ? nextsz - cur_sz : largest_field - sz
            if diff == 7
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
                push!(values, 0); push!(values, 0); push!(values, 0)
            elseif diff == 6
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
                push!(values, 0); push!(values, 0)
            elseif diff == 5
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
                push!(values, 0); push!(values, 0)
            elseif diff == 4
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
                push!(values, 0)
            elseif diff == 3
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
                push!(values, 0); push!(values, 0)
            elseif diff == 2
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
                push!(values, 0)
            elseif diff == 1
                push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
                push!(values, 0)
            end
            sz = (sz + diff) % largest_field == 0 ? 0 : (cur_sz < nextsz ? sz + diff : 0)
            cur_sz = 0
        end
    end
    return fields, values
end

macro struct(expr)
    !expr.args[1] || throw(ArgumentError("@struct is only applicable for immutable types"))
    exprs = filter(x->x.head !== :line, expr.args[3].args)
    fields, values = FlatBuffers.fieldlayout(expr.args[2], exprs...)
    expr.args[3].args = fields
    # generate convenience outer constructors if necessary
     # if there are nested structs or padding:
        # build an outer constructor that takes all direct, original fields
        # recursively flatten/splat all nested structs into one big args tuple
        # adding zeros for padded arguments
        # pass big, flat, args tuple to inner constructor
    T = expr.args[2]
    if any(x->!FlatBuffers.isbitstype(eval(current_module(), x.args[2])), exprs) ||
       length(fields) > length(exprs)
       exprs2 = map(x->FlatBuffers.isbitstype(eval(current_module(), x.args[2])) ? x.args[1] : x, exprs)
       sig = Expr(:call, T, exprs2...)
       body = Expr(:call, T, values...)
       outer = Expr(:function, sig, body)
    else
        outer = :(nothing)
    end
    return quote
        $(esc(expr))
        $(esc(outer))
    end
end

macro default(T, kwargs...)
    if eval(current_module(), T) <: Enum
        return quote
            # default{T<:Enum}(::Type{T}) = enumtype(T)(T(0))
            $(esc(:(FlatBuffers.default)))(::Type{$(esc(T))}) = $(esc(:(FlatBuffers.enumtype)))($(esc(T)))($(esc(kwargs[1])))
        end
    else
        ifblock = quote end
        for kw in kwargs
            push!(ifblock.args, :(if sym == $(QuoteNode(kw.args[1]))
                                      return $(kw.args[2])
                                  end))
        end
        return quote
            function $(esc(:(FlatBuffers.default)))(::Type{$(esc(T))}, TT, sym)
                $ifblock
                return FlatBuffers.default(TT)
            end
        end
    end
end

#TODO:
# handle default values
# handle unions
# handle id?
# handle deprecated
# nested_flatbuffer
macro table(expr)

end


text:
#------------------------------------------------------------------------------
"""
    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::BitArray(T)`: 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)

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

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

  vvM = falses(T,xNum)
  for i = 1:xNum                             #loop over inputs
    vvM[:,i] = any(isnan.(x[i]) .| ismissing.(x[i]),dims=dims)
  end

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

  return vvb

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

#------------------------------------------------------------------------------
"""
    FindNNDictPs(D...)

Find rows that have no missing values in any of the arrays in the dictionaries D...,
assuming all arrays have the same number of rows.

# Input
- `D::Dict`:    One (D) or several (D1,D2) dictionaries  with N different Tx? arrays.

"""
function FindNNDictPs(D...)

  N    = length(D)

  TM = hcat( [hcat([size(x,1) for x in values(D[i])]...) for i=1:N]... )  #size(x,1)
  all(TM .== TM[1]) ? T = TM[1] : error("different number of rows")

  vv2 = falses(T,N)
  for i = 1:N
    vv       = hcat([FindNNPs(x) for x in values(D[i])]...)
    vv2[:,i] = vec(all(vv,dims=2))
  end

  vvb = vec(all(vv2,dims=2))

  return vvb

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

#------------------------------------------------------------------------------
"""
    DictSelectRowsPs(D,vv)

Create new dictionary where we pick rows of D[i] to vv ()

"""
function DictSelectRowsPs(D,vv)
  Dnew = Dict()
  for (k,v) in D
    ndims(v) == 1 ? Dnew[k] = v[vv] : Dnew[k] = v[vv,:]
  end
  return Dnew
end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    FindNoNaNPs(Keepdim,x...)

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

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

# Output
- `vvb::BitArray(T)`: 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).

Paul.Soderlind@unisg.ch

"""
function FindNoNaNPs(Keepdim,x...)

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

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

  vvM = falses(T,xNum)
  for i = 1:xNum                             #loop over inputs
    vvM[:,i] = any(isnan.(x[i]) .| ismissing.(x[i]),dims=dims)
  end

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

  return vvb

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

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
using Plots
include("funcs.jl")

function rupture_plot(dir_name, begin_cycle, end_cycle)

    cd(string("../../erickson/output_files/", dir_name))
    @printf "found directory...\n"
    slip_file = collect(eachline(open("slip.dat":, "r")))
    y = get_line(slip_file, 2)[3:end]
    # slip data without first two lines
    slip_data = @view slip_file[3:end, :]

    get_cycle_indices(slip_data)
    # get indices in file where breaks are
    break_number, break_indices = get_break_indices(slip_data)
    # convert first cycle number to break number
    begin_break = 2*begin_cycle - 1
    # convert end cycle to break number
    # 2 * end_cycle + 1 is the interseismic phase of end_cycle + 1
    end_break = 2*end_cycle + 1

    # index of begin_cycle interseismic
    b_inter_index = break_indices[begin_break]
    # index of begin_cycle coseismic
    b_co_index = break_indices[begin_break + 1]
    # get last index of end_cycle coseismic
    end_index = break_indices[end_break] - 2

    # truncate slip_data around those indices
    slip_data = @view slip_data[b_inter_index:end_index, :]

end


text:
export monte_carlo_integrate, monte_carlo_𝒯integrate

function monte_carlo_integrate(func, a::Array{Float64,1}, b::Array{Float64,1}, ncomp::Int; args=(), method::Symbol=:vegas, kwargs...)
    integrand = generate_integrand(func, a, b, args)
    ndim = length(a)
    @eval begin
        $(method)($integrand,$ndim,$ncomp; $kwargs...)
    end
end

function generate_integrand(func, a, b, args)
    (x,f)-> f .= func(a + (b-a).*x, args...)
end

function monte_carlo_𝒯integrate(func, t, ndim::Int, ncomp::Int; args=(), method::Symbol=:vegas)
    integrand = generate_𝒯integrand_cuba(func, t, args=args)
    @eval begin
        $(method)($integrand,$ndim,$ncomp)
    end
end

function generate_𝒯integrand_cuba(func, t; args=())
    function integrand(x, f)
        y = 1.0
        x[1] = x[1]*t
        for i in 2:length(x)
            x[i] = x[i]*x[i-1]
            y = y*x[i-1]
        end
        f.=func(x, args...)*y*t
    end
end


text:
const gridnum =200
Tn = Array{Float64,2}(undef, gridnum, gridnum)

deriveM=[2*(gridnum-i) for i in 1:gridnum]

function chderive(xmin, xmax, input, ngrid)
    cder = similar(input)
    cder[ngrid] = Float64(0.0)
    cder[ngrid-1] = 2 * (ngrid - 1) * input[ngrid]
    for j = ngrid-2:-1:1
        cder[j] = cder[j+2] + 2 * j * input[j+1]
    end
    con = (Float64(2.0)) / (xmax - xmin)
    return con * cder
end

#
# coeffa = Tnt * u0
#
#
#
# coeffb = chder(Float64(0.0), ρmax, coeffa, gridnum)

for i = 1:gridnum
    for j = 1:gridnum
        Tn[j, i] = cos(π * (i - 1) * (j - Float64(0.5)) / gridnum)
    end
end

Tnt = (2 / gridnum * Tn)' |> collect
Tn[:, 1] .= Float64(0.5) * Tn[:, 1]

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection


text:
function demoSumProd()
    #DEMOSUMPROD Sum-Product algorithm test :
    # Variable order is arbitrary
    variables=1:5
    a, b, c, d, e=variables
    nstates=round.(Int64,3*rand(1,5)+2) # random number of states for each variable

    pot=Array{PotArray}(5)
    pot[1]=PotArray([a b],rand(nstates[[a b]]...))
    pot[2]=PotArray([b c d],rand(nstates[[b c d]]...))
    pot[3]=PotArray([c],rand(nstates[[c]]...))
    pot[4]=PotArray([e d],rand(nstates[[e d]]...))
    pot[5]=PotArray([d],rand(nstates[[d]]...))

    A = FactorGraph(pot)
    marg, mess, normconst=sumprod(pot,A)

    # check if this is correct:
    jointpot =prod(pot); V=length(pot)
    for i=1:V
        mtable=[condpot(marg[i]).content condpot(jointpot,i).content]
        println("marginal of variable $i:\n Factor Graph\t\t Raw summation");
        println(mtable)
    end
    tmp=sum(prod(pot))
    println("Normalisation constant:\nFactorGraph\t=$(normconst.content)\nRaw summation\t=$(tmp.content)")

end


text:
using RCall
using Base.Test,DataArrays,DataFrames

# write your own tests here
attenu = dataset(:attenu)
@test typeof(attenu) == DataFrame
dist = attenu[:dist]
@test typeof(dist) == DataVector{Float64}


text:
using jInv.Mesh
using Multigrid
using Multigrid.DomainDecomposition
using SparseArrays
using LinearAlgebra

speye = (n) -> SparseMatrixCSC(1.0I,n,n);

domain = [0.0, 1.0, 0.0, 1.0];
n = [64,64];
Mr = getRegularMesh(domain,n)
numCores = 4;
innerKatzIter = 5;
inner = 5;
nrhs = 2;

println("************************************************* Hybrid Keczmarz for DivSigGrad 2D ******************************************************");

m  = exp.(randn(prod(n)));
Ar = getNodalDivSigGradMatrix(Mr,vec(m));
Ar = Ar + 2e-1*opnorm(Ar,1)*speye(size(Ar,2));
Ar = convert(SparseMatrixCSC{Float64,Int64},Ar)

N = size(Ar,2);
b = Ar*rand(Float64,N,nrhs); b = b/norm(b);
x = zeros(Float64,N,nrhs);

HKparam = getHybridKaczmarz(Float64,Int64, Ar,Mr,[4,4], getNodalIndicesOfCell,0.8,numCores,innerKatzIter);
prec = getHybridKaczmarzPrecond(HKparam,Ar,nrhs);
x = FGMRES_relaxation(Ar,copy(b),x,inner,prec,1e-5*norm(b),true,numCores)[1]
println("Norm is: ",norm(Ar*x - b))

println("************************************************* Hybrid Keczmarz for Elasticity 2D ******************************************************");

mu = 2.0*ones(prod(Mr.n));
lambda = mu;
Ar = GetLinearElasticityOperator(Mr,mu,lambda);
Ar = Ar + 2e-1*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
N = size(Ar,2);
b = Ar*rand(N,nrhs); b = b/norm(b);
x = zeros(N,nrhs);

HKparam = getHybridKaczmarz(Float64,Int64, Ar,Mr,[4,4], getFacesStaggeredIndicesOfCellNoPressure,0.8,numCores,innerKatzIter);
prec = getHybridKaczmarzPrecond(HKparam,Ar,nrhs);
x = FGMRES_relaxation(Ar,copy(b),x,inner,prec,1e-5*norm(b),true,numCores)[1]
println("Norm is: ",norm(Ar*x - b))


text:
module Hrgnwake

#=
Copyright (c) 2021 p1scescom

This software is released under the MIT License.
http://opensource.org/licenses/mit-license.php
=#

export main

using Images, CUDA, Flux, ETLCDBReader, FixedPointNumbers, Colors, Statistics
using Flux.Data
using Flux.Losses: crossentropy
using Flux: onecold, onehot, onehotbatch

allhrgn = ['あ', 'い', 'う', 'え', 'お', 'か', 'が', 'き', 'ぎ', 'く', 'ぐ', 'け', 'げ', 'こ', 'ご', 'さ', 'ざ', 'し', 'じ', 'す', 'ず', 'せ', 'ぜ', 'そ', 'ぞ', 'た', 'だ', 'ち', 'ぢ', 'つ', 'づ', 'て', 'で', 'と', 'ど', 'な', 'に', 'ぬ', 'ね', 'の', 'は', 'ば', 'ぱ', 'ひ', 'び', 'ぴ', 'ふ', 'ぶ', 'ぷ', 'へ', 'べ', 'ぺ', 'ほ', 'ぼ', 'ぽ', 'ま', 'み', 'む', 'め', 'も', 'や', 'ゆ', 'よ', 'ら', 'り', 'る', 'れ', 'ろ', 'わ', 'を', 'ん']
hrgnl = length(allhrgn)

imgsize = (63,64)
batchsize = 30

# ETLCDB http://etlcdb.db.aist.go.jp/specification-of-etl-9?lang=ja
function gethrgndata(;dir="$(ENV["HOME"])/Downloads/ETL9B/", testlen = 20)
    etls = ETLCDBReader.getetl9b(dir)
    hrgnetls = filter(x -> x.charcode in allhrgn, etls)
    cuetlsdata = map(x -> Float32.(x.data), hrgnetls)

    train_x = cuetlsdata[1:end-hrgnl*testlen]
    train_y = onehotbatch(map(x -> x.charcode, hrgnetls)[1:end-hrgnl*testlen], allhrgn)

    test_x = cat(permutedims(cuetlsdata[end-hrgnl*testlen+1:end])..., dims=4)
    test_y = onehotbatch(map(x -> x.charcode, hrgnetls)[end-hrgnl*testlen+1:end], allhrgn)

    return train_x,train_y, test_x,test_y
end

function generatemodel()
    model = Chain(Conv((3,3), 1 => 8, relu, pad=1, stride=1),
                      flatten,
                      Dense(imgsize[1]*imgsize[2]*8, hrgnl),
                      softmax)
    return model
end

function train(model, datas, labels, batchsize = 20)
    loss(x,y) = begin
            crossentropy(model(x), y)
        end
    accuracy(x,y) = mean(onecold(model(x)) .== onecold(y))
    opt = ADAM()

    for _ in 1:1
        trainloader = DataLoader((datas, labels), batchsize=batchsize, shuffle=true)
        for (data,label) in trainloader
            cdata = cat(permutedims(data)..., dims=4)
            d = [(cdata, label)]
            Flux.train!(loss, params(model),d , opt)
            cdata = nothing
            d = nothing
        end
    end
end

# init function
function main(model = gpu(generatemodel()))
    train_x,train_y, test_x,test_y = gpu.(gethrgndata())
    train(model, train_x,train_y)
    return model
end

end # module


text:
function open(filename::AbstractString, new::Bool=false)
    Base.warn_once("ZipFile.open is deprecated, use ZipFile.Reader or ZipFile.Writer.")
    new ? Writer(filename) : Reader(filename)
end

const Dir = Reader
const WritableDir = Writer
const File = ReadableFile


text:
#!/usr/bin/env julia
using Dierckx
using Test
using Random: seed!

# Answers 'ans' are from scipy.interpolate,
# generated with genanswers.py script.

# -----------------------------------------------------------------------------
# Spline1D

x = [1., 2., 3.]
y = [0., 2., 4.]
spl = Spline1D(x, y; k=1, s=length(x))

yi = evaluate(spl, [1.0, 1.5, 2.0])
@test yi ≈ [0.0, 1.0, 2.0]
@test evaluate(spl, 1.5) ≈ 1.0
@test get_knots(spl) ≈ [1., 3.]
@test get_coeffs(spl) ≈ [0., 4.]
@test isapprox(get_residual(spl), 0.0, atol=1.e-30)

@test spl([1.0, 1.5, 2.0]) ≈ [0.0, 1.0, 2.0]
@test spl(1.5) ≈ 1.0

# test that a copy is returned by get_knots()
knots = get_knots(spl)
knots[1] = 1000.
@test get_knots(spl) ≈ [1., 3.]

# test ported from scipy.interpolate testing this bug:
# http://mail.scipy.org/pipermail/scipy-dev/2008-March/008507.html
x = [-1., -0.65016502, -0.58856235, -0.26903553, -0.17370892,
     -0.10011001, 0., 0.10011001, 0.17370892, 0.26903553, 0.58856235,
     0.65016502, 1.]
y = [1.,0.62928599, 0.5797223, 0.39965815, 0.36322694, 0.3508061,
     0.35214793, 0.3508061, 0.36322694, 0.39965815, 0.5797223,
     0.62928599, 1.]
w = [1.00000000e+12, 6.88875973e+02, 4.89314737e+02, 4.26864807e+02,
     6.07746770e+02, 4.51341444e+02, 3.17480210e+02, 4.51341444e+02,
     6.07746770e+02, 4.26864807e+02, 4.89314737e+02, 6.88875973e+02,
     1.00000000e+12]
spl = Spline1D(x, y; w=w, s=Float64(length(x)))
desired = [0.35100374, 0.51715855, 0.87789547, 0.98719344]
actual = evaluate(spl, [0.1, 0.5, 0.9, 0.99])
@test isapprox(actual, desired, atol=5e-4)

# test periodic
x = [1., 2., 3., 4., 5.]
y = [4., 1., 4., 1., 4.]
spl = Spline1D(x, y, periodic=true)

@test derivative(spl, 1) ≈ derivative(spl, 5)
@test derivative(spl, 1, nu=2) ≈ derivative(spl, 5, nu=2)

# tests for out-of-range
x = [0.0:4.0;]
y = x.^3

xp = range(-8.0, stop=13.0, length=100)
xp_zeros = Float64[(0. <= xi <= 4.) ? xi : 0.0 for xi in xp]
xp_clip = Float64[(0. <= xi <= 4.) ? xi : (xi < 0.0) ? 0.0 : 4. for xi in xp]

spl = Spline1D(x, y)
t = get_knots(spl)[2: end-1]  # knots, excluding those at endpoints
spl2 = Spline1D(x, y, t)

@test evaluate(spl, xp) ≈ xp_clip.^3
@test evaluate(spl2, xp) ≈ xp_clip.^3

# test other bc's
spl = Spline1D(x, y; bc="extrapolate")
@test evaluate(spl, xp) ≈ xp.^3
spl = Spline1D(x, y; bc="zero")
@test evaluate(spl, xp) ≈ xp_zeros.^3
spl = Spline1D(x, y; bc="error")
@test_throws ErrorException evaluate(spl, xp)

# test unknown bc
@test_throws ErrorException Spline1D(x, y; bc="unknown")

# test derivative
x = range(0, stop=1, length = 70)
y = x.^3
spl = Spline1D(x, y)
xt = [0.3, 0.4, 0.5]
@test derivative(spl, xt) ≈ 3xt.^2

# test integral
x = range(0, stop=10, length = 70)
y = x.^2
spl = Spline1D(x, y)
@test integrate(spl, 1.0, 5.0) ≈ 5.0^3/3 - 1/3

# test roots
x = range(0, stop=10, length = 70)
y = (x .- 4).^2 .- 1
spl = Spline1D(x, y)
@test roots(spl) ≈ [3, 5]

# test that show works.
io = IOBuffer()
show(io, spl)
seek(io, 0)
s = read(io, String)
@test s[1:9] == "Spline1D("

# test equality
seed!(0)
x = sort(rand(10))
y = rand(10)
sp1 = Spline1D(x, y)
sp2 = Spline1D(x, y)
sp3 = Spline1D(x.+1, y)
sp4 = Spline1D(x, y.+1)
@test sp1 == sp2
@test allunique([sp1,sp3,sp4])

# -----------------------------------------------------------------------------
# ParametricSpline

u = [1., 2., 3.]
x = [1. 2. 3.; 0. 2. 4.]
spl = ParametricSpline(u, x, k=1, s=size(x, 2))

xi = evaluate(spl, [1.0, 1.5, 2.0])
@test xi ≈ [1.0 1.5 2.0; 0.0 1.0 2.0]
@test evaluate(spl, 1.5) ≈ [1.5, 1.0]
@test get_knots(spl) ≈ [1., 3.]
@test get_coeffs(spl) ≈ [1.0 3.0; 0.0 4.0]
@test isapprox(get_residual(spl), 0.0, atol=1.e-30)

@test spl([1.0, 1.5, 2.0]) ≈ [1.0 1.5 2.0; 0.0 1.0 2.0]
@test spl(1.5) ≈ [1.5, 1.0]

# test that a copy is returned by get_knots()
knots = get_knots(spl)
knots[1] = 1000.
@test get_knots(spl) ≈ [1., 3.]

# test periodic
x = [23. 24. 25. 25. 24. 23.;
     13. 12. 12. 13. 13. 13.]
spl = ParametricSpline(x, periodic=true)
@test evaluate(spl, 0) ≈ evaluate(spl, 1)
@test derivative(spl, 0) ≈ derivative(spl, 1)
@test derivative(spl, 0, nu=2) ≈ derivative(spl, 1, nu=2)

# tests for out-of-range
u = 0.0:4.0
x = [u'.^2; u'.^3]

up = range(-8.0, stop=13.0, length = 100)
up_zeros = Float64[(0. <= ui <= 4.) ? ui : 0.0 for ui in up]
up_clip = Float64[(0. <= ui <= 4.) ? ui : (ui < 0.0) ? 0.0 : 4. for ui in up]

spl = ParametricSpline(u, x)
t = get_knots(spl)[2: end-1]  # knots, excluding those at endpoints
spl2 = ParametricSpline(u, x, t)

@test evaluate(spl, up) ≈ [up_clip'.^2; up_clip'.^3]
@test evaluate(spl2, up) ≈ [up_clip'.^2; up_clip'.^3]

# test other bc's
spl = ParametricSpline(u, x; bc="extrapolate")
@test evaluate(spl, up) ≈ [up'.^2; up'.^3]
spl = ParametricSpline(u, x; bc="zero")
@test evaluate(spl, up) ≈ [up_zeros'.^2; up_zeros'.^3]
spl = ParametricSpline(u, x; bc="error")
@test_throws ErrorException evaluate(spl, up)

# test unknown bc
@test_throws ErrorException ParametricSpline(u, x; bc="unknown")

# test derivative
u = range(0, stop=1, length = 70)
x = [u'.^2; u'.^3]
spl = ParametricSpline(u, x)
ut = [0.3, 0.4, 0.5]
@test derivative(spl, 0.3) ≈ [2*0.3, 3*0.3^2]
@test derivative(spl, ut) ≈ [2*ut'; 3*ut'.^2]
@test derivative(spl, 0.3, nu=2) ≈ [2.0, 6*0.3]
@test derivative(spl, ut, nu=2) ≈ [2*ones(3)'; 6*ut']

# test integral
u = range(0, stop=10, length = 70)
x = [u'.^2; u'.^3]
spl = ParametricSpline(u, x)
@test integrate(spl, 1.0, 5.0) ≈ [5.0^3/3 - 1/3, 5.0^4/4 - 1/4]

# test that show works.
io = IOBuffer()
show(io, spl)
seek(io, 0)
s = read(io, String)
@test s[1:17] == "ParametricSpline("

# test equality
seed!(0)
x = sort(rand(10))
y = rand(3, 10)
sp1 = ParametricSpline(x, y)
sp2 = ParametricSpline(x, y)
sp3 = ParametricSpline(x.+1, y)
sp4 = ParametricSpline(x, y.+1)
@test sp1 == sp2
@test allunique([sp1,sp3,sp4])

# test too many roots warning
x = (0:100)
y = (-1).^(0:100)
sp = Spline1D(x,y)
@test_logs (:warn,Regex("number of zeros exceeded")) roots(sp)

# -----------------------------------------------------------------------------
# Spline2D

# test linear
x = [1., 1., 1., 2., 2., 2., 3., 3., 3.]
y = [1., 2., 3., 1., 2., 3., 1., 2., 3.]
z = [0., 0., 0., 2., 2., 2., 4., 4., 4.]
spl = Spline2D(x, y, z; kx=1, ky=1, s=length(x))
tx, ty = get_knots(spl)
@test tx ≈ [1., 3.]
@test ty ≈ [1., 3.]
@test isapprox(get_residual(spl), 0.0, atol=1e-16)
@test evaluate(spl, 2.0, 1.5) ≈ 2.0
@test evalgrid(spl, [1.,1.5,2.], [1.,1.5]) ≈ [0. 0.; 1. 1.; 2. 2.]

# test 1-d grid arrays
@test evalgrid(spl, [2.0], [1.5])[1, 1] ≈ 2.0

# In this setting, lwrk2 is too small in the default run.
x = range(-2, stop=2, length = 80)
y = range(-2, stop=2, length = 80)
z = x .+ y
spl = Spline2D(x, y, z; s=length(x))
@test evaluate(spl, 1.0, 1.0) ≈ 2.0

# In this setting lwrk2 is too small multiple times!
# Eventually an error about s being too small is thrown.
seed!(0)
x = rand(100)
y = rand(100)
z = sin.(x) .* sin.(y)
@test_throws ErrorException Spline2D(x, y, z; kx=1, ky=1, s=0.0)

# test grid input creation
x = [0.5, 2., 3., 4., 5.5, 8.]
y = [0.5, 2., 3., 4.]
z = [1. 2. 1. 2.;  # shape is (nx, ny)
     1. 2. 1. 2.;
     1. 2. 3. 2.;
     1. 2. 2. 2.;
     1. 2. 1. 2.;
     1. 2. 3. 1.]
spl = Spline2D(x, y, z)

# element-wise output
xi = [1., 1.5, 2.3, 4.5, 3.3, 3.2, 3.]
yi = [1., 2.3, 5.3, 0.5, 3.3, 1.2, 3.]
ans = [2.94429906542,
       1.25537598131,
       2.00063588785,
       1.0,
       2.93952664,
       1.06482509358,
       3.0]
zi = evaluate(spl, xi, yi)
@test zi ≈ ans

zi = spl(xi, yi)
@test zi ≈ ans

# grid output
xi = [1., 1.5, 2.3, 4.5]
yi = [1., 2.3, 5.3]
ans = [2.94429906542  1.16946130841  1.99831775701;
       2.80393858478  1.25537598131  1.99873831776;
       1.67143209613  1.94853338542  2.00063588785;
       1.89392523364  1.8126946729  2.01042056075]
zi = evalgrid(spl, xi, yi)
@test zi ≈ ans

# Test 2-d integration
test2d_1(x, y) = 1 - x^2 -y^2
test2d_2(x, y) = cos(x) + sin(y)
test2d_3(x, y) = x*exp(x-y)
for (f, domain, exact) in [(test2d_1, (0.0, 1.0, 0.0, 1.0), 1.0/3.0),
                           (test2d_2, (0.0, pi, 0.0, pi), 2.0*pi),
                           (test2d_3, (0.0, 1.0, 0.0, 1.0), (ℯ-1.0)/ℯ)]
    (x0, x1, y0, y1) = domain

    # define grids for x and y dimensions:
    npoints = 50
    xgrid = range(x0, stop=x1, length = npoints)
    ygrid = range(y0, stop=y1, length = npoints)

    fxygrid = zeros(npoints, npoints)
    for (j, y) in enumerate(ygrid)
        local j, y
        for (i, x) in enumerate(xgrid)
            local i, x
            fxygrid[i,j] = f(x, y)
        end
    end

    spl1 = Spline2D(xgrid, ygrid, fxygrid)
    @test isapprox(integrate(spl1, x0, x1, y0, y1), exact, atol=1e-6)
end

# test equality
seed!(0)
x = sort(rand(10))
y = sort(rand(10))
z = rand(10,10)
sp1 = Spline2D(x, y, z)
sp2 = Spline2D(x, y, z)
sp3 = Spline2D(x.+1, y, z)
sp4 = Spline2D(x, y.+1, z)
sp5 = Spline2D(x, y, z.+1)
@test sp1 == sp2
@test allunique([sp1, sp3, sp4, sp5])

println("All tests passed.")


text:
module GenInvGaussian

export GeneralizedInverseGaussian

import Base.mean
import Base.rand
using Distributions
import Distributions: mean, var, mode, mgf, cf, rand, params, minimum, maximum
# import Distributions: @check_args, @distr_support

immutable GeneralizedInverseGaussian <: ContinuousUnivariateDistribution
    p::Float64
    a::Float64
    b::Float64

    function GeneralizedInverseGaussian(p::Real, a::Real, b::Real)
#        @check_args(GeneralizedInverseGaussian, a > zero(a) && b > zero(b))
        new(p, a, b)
    end
end

minimum(::GeneralizedInverseGaussian) = 0.0
maximum(::GeneralizedInverseGaussian) = Inf

#### Parameters

params(d::GeneralizedInverseGaussian) = (d.p, d.a, d.b)

#### Statistics

function mean(d::GeneralizedInverseGaussian)
    (p, a, b) = params(d)
    ω = sqrt(a*b)
    sqrt(b/a) * besselk(p + 1, ω) / besselk(p, ω)
end

function var(d::GeneralizedInverseGaussian)
    (p, a, b) = params(d)
    ω = sqrt(a*b)
    num = besselk(p, ω)*besselk(p+2, ω) - besselk(p+1, ω)*besselk(p+1, ω)
    denom = besselk(p, ω) * besselk(p, ω)
    b/a * num / denom
end

function mode(d::GeneralizedInverseGaussian)
    (p, a, b) = params(d)
    m = (p - 1) + sqrt((p - 1)*(p - 1) + a*b)
    if m < 0
        error("Distribution has no mode!")
    end
    m / a
end

function mgf(d::GeneralizedInverseGaussian, t::Real)
    (p, a, b) = params(d)
    at = a - 2t
    ω = sqrt(a*b)
    ωt = sqrt(at*b)
    (a/at)^(p/2) * besselk(p, ωt)/besselk(p, ω)
end

function cf(d::GeneralizedInverseGaussian, t::Real)
    (p, a, b) = params(d)
    at = a - 2im * t
    ω = sqrt(d.a*d.b)
    ωt = sqrt(at*b)
    (a/at)^(p/2) * besselk(p, ωt)/besselk(p, ω)
end

#### Evaluation & Sampling

# Simple rejection sampling using a gamma distribution as a proposal distribution.
# Probably not the best solution out there. However, for p >= 1, the expected number of
# rejection is lower than sqrt(2).
function rand(d::GeneralizedInverseGaussian)
    (p, a, b) = params(d)

    # Handle the case p < 0
    if p < 0
        p = -p
        a, b = b, a
        invert = true
    else
        invert = false
    end

    ω = sqrt(a*b)

    r = p/2 + sqrt(p*p + a*b)/2
    dproposal = Gamma(r, 2/a)
    daccept = Uniform()

    xm = b / ( r - p) / 2

    acceptrate(x::Real) = (x./xm)^(p - r) * exp(-b/2*(1./x - 1./xm))

    while true
        sproposal = rand(dproposal)
        saccept = rand(daccept)

        if saccept < acceptrate(sproposal)
            return invert ? 1/sproposal : sproposal
        end
    end
end

end # module


text:
using Test

abstract type SolutionIO end

Base.close(sio::SolutionIO) = error("Base.close(::SolutionIO) not implemented for ", typeof(sio))

save_attributes(sio::SolutionIO, ::Solution) = error("save_attributes(::SolutionIO, ::Solution) not implemented for ", typeof(sio))
save_solution(sio::SolutionIO, ::Solution) = error("save_solution(::SolutionIO, ::Solution) not implemented for ", typeof(sio))
load_solution(sio::SolutionIO, ::Solution) = error("load_solution(::SolutionIO, ::Solution) not implemented for ", typeof(sio))

function test_interface(sio::SolutionIO, sol::Solution)

    @test_nowarn save_attributes(sio, sol)
    @test_nowarn save_solution(sio, sol)
    @test_nowarn load_solution(sio, sol)
    @test_nowarn close(sio)

end


text:
# using MLDatasets: MNIST
using Base.Iterators
import Pkg; Pkg.add("Augmentor"); using Augmentor
using DataFrames

function test_lenet()

    augm = CropSize(28,28)
    trn, vld = mk_image_minibatch("../data/flowers",
                4; split=true, fr=0.2,
                balanced=false, shuffle=true,
                train=true,
                aug_pipl=augm, pre_proc=preproc_imagenet)

    lenet = Classifier(Conv(5,5,3,20),
                    Pool(),
                    BatchNorm(),
                    Conv(5,5,20,50),
                    BatchNorm(trainable=true, channels = 50),
                    Pool(),
                    BatchNorm(trainable=true),
                    Flat(),
                    Dense(800,512),
                    Linear(512, 512, actf=relu),
                    Dense(512,3,actf=identity)
            )

    mdl = tb_train!(lenet, Adam, trn, vld, epochs=1,
            acc_fun=accuracy,
            eval_size=0.25, eval_freq=2, mb_loss_freq=100,
            tb_name="test_run", tb_text="NNHelferlein example")

    acc = accuracy(mdl, data=vld)

    tst = mk_image_minibatch("../data/flowers",
                4; split=false, fr=0.2,
                balanced=false, shuffle=true,
                train=false,
                aug_pipl=augm, pre_proc=nothing)
    p = predict_imagenet(mdl, tst, top_n=2)
    return acc isa Real && acc <= 1.0 &&
           isdir("logs") &&
           size(p) == (3,8)
end

function acc_fun(mdl; data=data)
    a = 0.0
    for (x,y) in data
        a += mean(abs2, mdl(x) .- y)
    end
    return a/length(data)
end

function test_mlp()

        trn = DataFrame(x1=randn(16), x2=randn(16),
                        x3=randn(16), x4=randn(16),
                        x5=randn(16), x6=randn(16),
                        x7=randn(16), x8=randn(16),
                        y=collect(range(0, 1, length=16)))

        mb = dataframe_minibatches(trn, size=4)

        mlp = Regressor(Dense(8,8, actf=relu),
                         Dense(8,8),
                         Dense(8,1, actf=identity))

        mlp = tb_train!(mlp, Adam, mb, epochs=10, acc_fun=nothing,
                lr=0.001, lr_decay=0.0001, lrd_freq=5)

        mlp = tb_train!(mlp, Adam, mb, epochs=10, acc_fun=acc_fun,
                lr=0.001, lr_decay=0.0001, lrd_freq=5)

        acc = NNHelferlein.calc_acc(mlp, acc_fun, data=mb)
        return acc isa Real
end

function test_signatures()

        trn = DataFrame(x1=randn(16), x2=randn(16),
                        x3=randn(16), x4=randn(16),
                        x5=randn(16), x6=randn(16),
                        x7=randn(16), x8=randn(16),
                        y=collect(range(0, 1, length=16)))

        mb = dataframe_minibatches(trn, size=4)

        mlp = Regressor(Dense(8,8, actf=relu),
                         Dense(8,8),
                         Dense(8,1, actf=identity))
        mlp = tb_train!(mlp, Adam, mb, epochs=1, acc_fun=nothing)

        y = mlp(rand(Float32, 8,4))
        test_sign = size(y) == (1,4)

        loss = mlp(rand(Float32, 8,4), rand(Float32, 1,4))
        test_sign = test_sign && typeof(loss) <: Real

        loss = mlp(mb)
        test_sign = test_sign && typeof(loss) <: Real

        loss = mlp(first(mb))
        test_sign = test_sign && typeof(loss) <: Real

        x,y = first(mb)
        loss = mlp(x,y)
        test_sign = test_sign && typeof(loss) <: Real

        return test_sign
end

function test_decay_cp()

        trn = DataFrame(x1=randn(16), x2=randn(16),
                        x3=randn(16), x4=randn(16),
                        x5=randn(16), x6=randn(16),
                        x7=randn(16), x8=randn(16),
                        y=collect(range(0, 1, length=16)))

        mb = dataframe_minibatches(trn, size=4)

        chain = Chain(Dense(8,8, actf=relu),
                      Dense(8,8))
        mlp = Regressor(chain,
                        Dense(8,1, actf=identity))

        mlp = tb_train!(mlp, Adam, mb, epochs=2, acc_fun=nothing,
                cp_freq=1, lr=0.01, lr_decay=0.99, l2=1e-6)
        acc = NNHelferlein.calc_acc(mlp, acc_fun, data=mb)
        return acc isa Real
end

function test_hamming()

        p = [0  0  0  0  0
             4  4  3  4  5
             2  2  2  2  5
             1  3  1  3  3
             0  0  0  0  0]
        t = [0  0  0  0  0
             4  4  5  3  4
             3  1  2  4  4
             2  5  5  4  2
             0  0  0  0  0]

        acc = hamming_dist(p,t)
        hd = hamming_dist(p,t, accuracy=true)

        return acc ≈ 2.4 && hd ≈ 0.2
end


text:
using Statistics

function portwine_data()
    A1 = [
    7 5 7 5 5 6 5 6;
    0 6 2 7 7 8 4 6;
    5 6 5 7 6 6 10 6;
    8 3 5 4 4 1 3 5
    ]
    A2 = [
    4 3 3 1 2 1 0 2;
    0 6 3 6 5 5 4 6;
    5 5 7 3 5 4 2 4
    ]
    A3 = [
    7 2 6 2 5 3 2 4;
    4 0 3 0 1 0 0 0;
    2 6 4 6 5 5 4 4;
    6 6 7 4 6 5 3 5
    ]
    A4 = [
    9 8 10 7 8 8 6 8;
    7 6 6 7 7 8 5 9;
    9 7 7 8 8 10 10 10
    ]
    m = 4
    A = [A1', A2', A3', A4']
    A = map(M -> M .- mean(M, dims=1), A)
    maxdiff_optim = [419.65133740484237, 542.3275506383522, 568.19290485317]
    maxbet_optim  = [605.6934523038526, 779.7572614382912, 829.7087289316665]
    A, maxdiff_optim, maxbet_optim
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.58980225336104 <= q <= 16.332536463331582)

                     add_NL_constraint(m, :(sigmoid(0.1097699546178803 + 0.5720816873107069*sigmoid(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + 0.11808627755560508*sigmoid(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.12730937872071246*sigmoid(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.6166642606025916*sigmoid(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + 0.45996993866072433*sigmoid(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + sigmoid(-0.8288330741680832 + -0.5580861219922526*sigmoid(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.8786358511296881*sigmoid(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.6915444540383735*sigmoid(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.5554450106282389*sigmoid(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.20937059223327026*sigmoid(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + sigmoid(0.13615042538554878 + 0.08502572267265318*sigmoid(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.9980047609363689*sigmoid(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.5678523898140919*sigmoid(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.29618367003382495*sigmoid(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.7256587322261439*sigmoid(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + sigmoid(-0.16139551949919406 + -0.5886735107938796*sigmoid(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + -0.02630736114123211*sigmoid(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.03441572289730743*sigmoid(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + 0.3935924617025486*sigmoid(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + 0.3425056157482871*sigmoid(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) + sigmoid(-0.5329713534130147 + -0.3183971019485696*sigmoid(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3])) + 0.5847814320797045*sigmoid(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3])) + -0.7330877632484647*sigmoid(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3])) + -0.37731199331975374*sigmoid(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3])) + -0.08817423127621415*sigmoid(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module Optimization

end # module


text:
@testset "Atmosphere structure" begin
    forced_date = Dates.DateTime("2021-09-15T16:24:00.929")

    # Testing Atmosphere with some random values:
    @test Atmosphere(date = forced_date, T = 25, Wind = 5, Rh = 0.3) == Atmosphere{Float64,Dates.DateTime,Float64}(Dates.DateTime("2021-09-15T16:24:00.929"), 1.0, 25.0,
5.0, 101.325, 0.3, 400.0, 0.9540587244435038, 3.180195748145013, 2.2261370237015092, 1.1838896840018194, 2.441875e6, 0.06757907523556121, 0.5455578187331258, 0.19009500927530176, 9999.9, 9999.9, 9999.9, 9999.9, 9999.9, 9999.9)

    # Testing Rh with values given in %:
    @test_logs (:warn, "Rh should be 0 < Rh < 1, assuming it is given in % and dividing by 100") Atmosphere(T = 25, Wind = 5, Rh = 30)
    # Testing Rh with values given with wrong value:
    @test_logs (:error, "Rh should be 0 < Rh < 1, and its value is 300") Atmosphere(T = 25, Wind = 5, Rh = 300)

    test_met = Atmosphere(T = 25, Wind = 5, Rh = 30)
    @test test_met.Rh == 0.3
end;


text:
filter((x) -> length(x) >= row, trial)

for x in 1:length(quadrantdf[:, 1])
    if quadrantdf[x, 3] == "quad 1"
        push!(quad1, quadrantdf[x, :])
    elseif quadrantdf[x, 3] == "quad 2"
        push!(quad2, quadrantdf[x, :])
    elseif quadrantdf[x, 3] == "quad 3"
        push!(quad3, quadrantdf[x, :])
    elseif quadrantdf[x, 3] == "quad 4"
        push!(quad4, quadrantdf[x, :])
    else
    end
end


text:
using DataFrames
using Statistics
using CubicSplines
using DSP
using Dates

export standardize
export spline_fill
export digital_filter
export fill_na
export date_col_name

"return boolean vector of whether df columns are numeric"
numeric_cols(df::DataFrame) = [(Float64 <: eltype(df[!, i])) | (Int64 <: eltype(df[!, i])) for i in 1:ncol(df)]

"""
    standardize all numeric columns in a dataframe. _newxᵢ = (xᵢ - μ) / σ²_
    parameters:
        df : DataFrame
            dataframe to standardize
    returns: DataFrame
        dataframe with standardized columns
"""
function standardize(df::DataFrame)
    tmp = copy(df)
    is_numeric = numeric_cols(tmp)
    for i in 1:ncol(tmp)
        if is_numeric[i]
            tmp[!, i] = (tmp[!, i] .- mean(skipmissing(tmp[!, i]))) / std(skipmissing(tmp[!, i]))
        end
    end
    return tmp
end

"""
    replace intermediate missings in a series with cubic spline
    parameters:
        series : Array
            series to fill missings
    returns: Array{Float64}
        series with intermediate missings filled in
"""
function spline_fill(series::Array)
    tmp = copy(series) |> x->
        [ismissing(a) ? missing : Float64(a) for a in x]
    xs = findall(.!ismissing.(tmp))
    ys = tmp[xs]
    spline = CubicSpline(xs, ys)
    tmp[xs[1]:xs[end]] = spline[xs[1]:xs[end]]
    return tmp
end

"""
    fill in missing tails of series using 1-D digital filter
    parameters:
        series : Array
            series to fill leading/trailing missings
        k : Int
            rational transfer function argument's numerator
    returns: Array
        series with filled missings
"""
function digital_filter(series::Array)
    tmp = copy(series) |> x->
        [ismissing(a) ? missing : Float64(a) for a in x]
    missing_indices = findall(ismissing.(tmp))
    tmp[missing_indices] .= median(skipmissing(tmp))
    tmp = tmp .|> Float64

    responsetype = Lowpass(0.9)
    designmethod = Butterworth(1)
    ma = filt(digitalfilter(responsetype, designmethod), tmp)

    tmp[missing_indices] = ma[missing_indices]
    return tmp
end

"""
    remove missing/NAs in a dataframe
    parameters:
        df : DataFrame
            dataframe to fill NAs
    returns: Dict
        Dict(:output => df with NAs filled, :na_indices => df of Bool is na)
"""
function fill_na(df::DataFrame)
    tmp = copy(df)
    # adding a false at end for the index column
    is_numeric = [numeric_cols(tmp); false]
    na_indices = ismissing.(tmp)
    tmp[!, :index] = 1:nrow(tmp)

    # if more than 80% of columns are missing in the beginning, drop the row
    threshold = 0.8
    mask = ([(x |> Array .|> ismissing |> sum) / ncol(tmp) for x in eachrow(tmp)] .<= threshold)
    mask = (mask .== 1) .| ([mask[1:end-1]; 1] .== 1)
    tmp = tmp[mask, :]

    # if all entries are missing at the end, drop the row
    mask = [!all(ismissing, Array(x)) for x in eachrow(tmp[!, is_numeric])]
    tmp = tmp[mask, :]

    # fill missing values between with cubic spline, fill head/tail with digital filter
    for i in 1:ncol(tmp)
        if is_numeric[i]
            tmp[!, i] = spline_fill(tmp[!, i])
            tmp[!, i] = digital_filter(tmp[!, i])
        end
    end

    return Dict(:output => tmp[!, Not(:index)], :na_indices => na_indices[tmp.index,:])
end

"""
    generate matrix of monthly and quarterly variables
    parameters:
        dates : Array{Dates.Date}
            array of corresponding dates for the matrix
        df : DataFrame
            dataframe to get monthly and quarterly variables from
    returns: Array{Bool, 1}
        array of 1's (is quarterly) and 0's (is not quarterly, i.e. monthly)
"""
function gen_monthly_quarterly(dates::Array{Dates.Date}, df::DataFrame)::Array{Bool, 1}
    is_quarterly = []
    for i in 1:ncol(df)
        quarterly = ([Month(x).value for x in dates[.!ismissing.(df[!, i])]] |> Set) |> x->
            isequal(x, Set([3,6,9,12]))
        push!(is_quarterly, quarterly)
    end
    return is_quarterly
end

"""
    return name of first date column of a dataframe
    parameters:
        df : DataFrame
            dataframe containing a date column
    returns: Symbol
        name of the first date column
"""
function date_col_name(df::DataFrame)
    try
        return findall(x->x==Dates.Date, eltype.(eachcol(df)))[1] |> i->
            Symbol(names(df)[i])
    catch e
        println("No column of type Dates.Date")
    end
end

"""
    given a dataframe, create a lagged dataframe where lag months back will be set to missing
    parameters:
        x : DataFrame
            dataframe to be artificially lagged
        lag : Int
            number of months to lag dataframe by
    returns: DataFrame
        lagged dataframe
"""
function create_lag(x, lag)
    new = copy(x)
    for i in 2:ncol(x)
        last_data_index = findall(!ismissing, x[!, i])[end]
        new[(last_data_index-lag+1):last_data_index, i] .= missing
    end
    return new
end


text:
# testing the autodiff and rootfinding packages in Julia
# please see the notebook in this folder

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
# Copyright (c) 2021 California Institute of Technology (“Caltech”).
# U.S. Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# - Neither the name of Caltech nor its operating division, the Jet Propulsion
#   Laboratory, nor the names of its contributors may be used to endorse or
#   promote products derived from this software without specific prior written
#   permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

############### Angular cryofunks ########################

function get_angular_cryofunks(mask; G0=nothing, nside_hires=npix2nside(length(mask)))
    if G0 == nothing
        @time G0 = get_green(npix2nside(length(mask)), nside_hires=nside_hires)
    end
    G = Hermitian(collect(G0[mask.>0,mask.>0]))

    #fsky = sum(mask .> 0) / length(mask)
    #G = Hermitian(G + I / (4*π*fsky*length(mask)))

    G = monopolize_green(G)

    B = get_B(mask)
    @time ℓ, Z = get_cryobasis(G)

    ℓ[1] = 0  # monopole
    Z[:,1] .= mean(Z[:,1])  # monopole

    ℓ, Z = sort_cryobasis(ℓ, Z)
    Y = B \ Z
    Yinv = Z' * B

    return ℓ, Y, Yinv
end

function get_cryobasis(G)
    f = eigen(G)
    λ = f.values  # already sorted smallest to largest
    Z = f.vectors

    @show λ[1:4] λ[end-4:end]
    epsmax = eps(maximum(abs.(λ)))

    if maximum(abs.(imag(λ))) <= epsmax
        # assume this is numerical noise, and imag really vanishes, as for all
        # Hermitian matrices.
        λ = real(λ)
    end

    if abs(λ[1]) <= epsmax
        # assume this is numerical noise, and it should really vanish, as for
        # the monopole.
        λ[1] = 0
    end

    @show λ[1] λ[end]
    #λ = complex(λ)

    ℓ = @. sqrt(1 / λ + 1/4) - 1/2
    #sel = @. imag(ℓ) != 0
    #ℓ[sel] .= abs.(ℓ[sel]) .+ maximum(real(ℓ))
    #ℓ = real(ℓ)

    return ℓ, Z
end

# simple form to fill in the Green's matrix
function get_green1(mask; kwargs...)  # 'kwargs' are ignored
    bmask = collect(mask .> 0)

    nside = npix2nside(length(bmask))
    npix = nside2npix(nside)
    Ωᵢ = 4 * π / npix

    Gii_same = Ωᵢ / (4 * π) * (1 - log(Ωᵢ/(2*π)))
    θ = fill(NaN, npix)
    ϕ = fill(NaN, npix)
    for i=1:npix
        θ[i], ϕ[i] = pix2angRing(Resolution(nside), i)
    end
    sinθ = sin.(θ)

    @show npix^2 * 8/1024^3
    Gij = fill(NaN, npix, npix)
    for j=1:npix, i=j:npix
        (bmask[i] && bmask[j]) || continue
        if i == j
            Gij[i,j] = Gii_same
        else
            Δθ = θ[j] - θ[i]
            Δϕ = ϕ[j] - ϕ[i]
            s²ρ2 = sin(Δθ/2)^2 + sinθ[i]*sinθ[j]*sin(Δϕ/2)^2
            Gij[i,j] = - log(2 * s²ρ2) * Ωᵢ / (4 * π)
            Gij[j,i] = Gij[i,j]'
        end
    end

    return Hermitian(Gij)
end

# integrate over each pixel
function get_green3(mask; nside_hires=1*npix2nside(length(mask)))
    T = Float64
    nside = npix2nside(length(mask))
    npix = nside2npix(nside)
    npix_hires = nside2npix(nside_hires)

    bmask = collect(mask .> 0)

    Ωᵢhires = 4 * π / npix_hires
    Gii_same = Ωᵢhires / (4 * π) * (1 - log(Ωᵢhires/(2*π)))
    θ = fill(NaN, npix_hires)
    ϕ = fill(NaN, npix_hires)
    for i=1:npix_hires
        θ[i], ϕ[i] = pix2angRing(Resolution(nside_hires), i)
    end
    sinθ = sin.(θ)
    cosθ = cos.(θ)

    pixmap = [ang2pix(mask, θ[i], ϕ[i]) for i=1:npix_hires]

    @show nside nside_hires
    @show npix npix_hires
    @show npix^2 * 8/1024^3
    @show npix_hires^2 * 8/1024^3
    Gmn = fill(T(0), npix, npix)
    for j=1:npix_hires
        n = pixmap[j]
        bmask[n] || continue
        for i=1:npix_hires
            m = pixmap[i]
            bmask[m] || continue
            if i == j
                Gmn[m,n] += Gii_same
            else
                Δθ = θ[j] - θ[i]
                Δϕ = ϕ[j] - ϕ[i]

                # arcsin formula
                s²ρ2 = sin(Δθ/2)^2 + sinθ[i]*sinθ[j]*sin(Δϕ/2)^2
                Gij = - log(2 * s²ρ2) * Ωᵢhires / (4 * π)

                # arctan formula
                r² = (sinθ[j] * sin(Δϕ))^2 + (sinθ[i] * cosθ[j] - cosθ[i] * sinθ[j] * cos(Δϕ))^2
                d = cosθ[i] * cosθ[j] + sinθ[i] * sinθ[j] * cos(Δϕ)
                ρ = atan(sqrt(r²), d)  # Don't do a naive atan! atan(x) looses info.
                Gij2 = - Ωᵢhires / (4*π) * log(2 * sin(ρ/2)^2)

                #@show Gij Gij2 θ[i],ϕ[j] Δθ,Δϕ 2asin(√s²ρ2)-π/2 ρ+π/2
                #@show tan(2asin(√s²ρ2)) tan(ρ)
                #tanρ = sqrt(r²) / d
                #@show s²ρ2 tanρ^2 / (1 + sqrt(1 + tanρ^2)) / sqrt(1 + tanρ^2) / 2
                false && if abs(Gij - Gij2) > 1e-10 * (abs(Gij) + abs(Gij2)) / 2
                    θ1 = big(θ[i])
                    θ2 = big(θ[j])
                    ϕ1 = big(ϕ[i])
                    ϕ2 = big(ϕ[j])
                    Δθ = θ2 - θ1
                    Δϕ = ϕ2 - ϕ1

                    s²ρ2 = sin(Δθ/2)^2 + sin(θ1)*sin(θ2)*sin(Δϕ/2)^2
                    Gij = - log(2 * s²ρ2) * Ωᵢhires / (4 * π)

                    # arctan formula
                    r² = (sin(θ2) * sin(Δϕ))^2 + (sin(θ1) * cos(θ2) - cos(θ1) * sin(θ2) * cos(Δϕ))^2
                    d = cos(θ1) * cos(θ2) + sin(θ1) * sin(θ2) * cos(Δϕ)
                    ρ = atan(sqrt(r²) / d)
                    @error "Mismatch!" θ1 ϕ1 θ2 ϕ2 Δθ Δϕ atan(tan(2asin(√s²ρ2))) atan(tan(ρ))
                    error("Wrong!")
                end

                Gmn[m,n] += Gij2
            end
        end
    end
    Ωᵢ = 4 * π / npix
    Gmn .*= Ωᵢhires / Ωᵢ
    return Hermitian(Gmn)
end

function calc_green(mask; nside_hires=npix2nside(length(mask)))
    bmask = collect(mask .> 0)
    mask = HealpixMap{eltype(mask), Healpix.RingOrder}(mask)
    return Hermitian(collect(get_green3(mask, nside_hires=nside_hires)[bmask,bmask]))
end

function get_green(nside::Int, fname=nothing; nside_hires=nside)
    if fname == nothing
        fname = "$cachedir_greensmatrix/Gang_nside$(nside)_$(nside_hires).bin"
    end
    npix = nside2npix(nside)
    if !isfile(fname)
        @timeit "Angular Green's matrix (nside,nside_hires)=($nside,$nside_hires)" G = calc_green(fill(1.0, npix), nside_hires=nside_hires)
        write(fname, G)
    end
    return Hermitian(read!(fname, fill(NaN, npix, npix)))
end

function monopolize_green(G)
    u = fill(1.0, size(G,1))
    u = u / sqrt(u'*u)
    P = Hermitian(I - u*u')
    G′ = P * G * P'
    if typeof(G) <: Hermitian
        # Note: Here we set the type, because that will significantly speed up
        # the eigenvector decomposition. It really should be considered a bug
        # in Julia that P*G*P' does not automatically result in the best type.
        G′ = Hermitian(G′)
        G′.data .= G′
    end
    @show u'*u typeof(u*u')
    @show typeof(G′)
    return G′
end

function sort_cryobasis(ℓ, Z)
    perm = sortperm(real(ℓ))
    ℓ .= ℓ[perm]
    Z[:,:] .= Z[:,perm]
    return collect(ℓ), collect(Z)
end

function get_B(mask)
    nside = npix2nside(length(mask))
    npix = size(mask,1)
    Ωᵢ = 4 * π / npix
    B = √Ωᵢ* I
    return B
end

######################### Angular Functions, integrate directly ######

function Healpix.pix2angRing(nside::Integer, idxs::AbstractArray)
    reso = Resolution(nside)
    θ = fill(NaN, length(idxs))
    ϕ = fill(NaN, length(idxs))
    for i=1:length(idxs)
        θ[i], ϕ[i] = pix2angRing(reso, idxs[i])
    end
    return θ, ϕ
end

function get_angular_continuous_cryofunks(mask; nside_hires=nside_hires)
    npix = length(mask)
    npix_hires = nside2npix(nside_hires)
    nside = npix2nside(npix)
    Ω = 4*π / npix

    ℓ, Y = calc_pixwin(nside, nside_hires=nside_hires)

    Yinv = Y' * Ω  # one of the eigenvalues vanishes, but this is all we need.

    return ℓ, Y, Yinv
end

# vim: set sw=4 et sts=4 :