text:
################################################################################
# Convert the csv files gotten from ContactNets into the jld2 format
################################################################################

function convert_csv()
	for k = 0:569
		csv_path = joinpath(module_dir(), "examples", "real2sim", "data", "tosses_csv", "$(k).csv")
		jld2_path = joinpath(module_dir(), "examples", "real2sim", "data", "tosses_jld2", "$(k).jld2")
		toss = split.(readlines(csv_path), ",")
		toss = [parse.(Float64, t) for t in toss]
		jldsave(jld2_path; toss=toss)
	end
end

# convert_csv()

################################################################################
# Convert the JLD2 files into trajectories and store them into a dataset
################################################################################

function toss2z(toss, timestep; s=1)
	N = length(toss)
	M = Int((N-1 - (N-1)%S) / S)
	z = []
	for i = 1:M
		vec1 = toss[1+S*(i-1)]
		vec2 = toss[1+S*i]

		x1 = vec1[1:3]
		q1 = vec1[4:7]

		x2 = vec2[1:3]
		q2 = vec2[4:7]

		v15 = (x2 - x1) / timestep
		ω15 = angular_velocity(Quaternion(q1...),
			Quaternion(q2...), timestep)

		z2 = [x2; v15; q2; ω15]
		push!(z, z2)
	end
	return z
end

function build_pairs(z)
    pairs = []
	N = length(z)
    for t = 1:N-1
        z1 = z[t]
        z2 = z[t+1]
        pair = [z1, z2]
        push!(pairs, pair)
    end
    return pairs
end

function generate_hardware_dataset(;N::Int=10,
		sleep_ratio = 0.0,
		show_contact = true,
		s=1,
		)
	H = 1.00
	timestep= 1/148 * S
	gscaled = -9.81*20

    mechanism = get_mechanism(:block, timestep=timestep, gravity=gravityscaled);
    trajs = []
	pairs = []
    for i = 1:N
		file = jldopen(joinpath(module_dir(), "examples", "real2sim", "data", "tosses_jld2", "$(i).jld2"))
		toss = file["toss"]
		# z = toss2z(toss[1:7], timestep)
		# z = toss2z(toss[end-10:end], timestep)
		z = toss2z(toss, timestep, S=S)
		storage = generate_storage(mech, z)
		push!(pairs, build_pairs(z)...)
        push!(trajs, storage)
        visualize(mechanism, storage, vis=vis, show_contact=show_contact)
		sleep(H*sleep_ratio)
    end
	data = [
		0.2, 0,0,0, +1, +1, -1,
		0.2, 0,0,0, +1, -1, -1,
		0.2, 0,0,0, -1, +1, -1,
		0.2, 0,0,0, -1, -1, -1,
		0.2, 0,0,0, +1, +1, +1,
		0.2, 0,0,0, +1, -1, +1,
		0.2, 0,0,0, -1, +1, +1,
		0.2, 0,0,0, -1, -1, +1]
    params = Dict(:N => N, :H => H, :timestep => timestep, :g => gscaled, :data => data)
    jldsave(joinpath(@__DIR__, "dataset", datafilename(:hardwarebox; N = N, S = S));
        params=params, trajs=trajs, pairs=pairs)
    return nothing
end


text:
name = "Jacob"
num = 42

if true
    println("True!")
elseif false
    println("False!")
else
    println("WTF?!")
end

for x in 1:10
    y = x^2
    println("$(x) squared is $(y)")
end

function area(width, height)
    return width * height
end

println(area(12, 34))

ext: jl

max_stars_repo_name: jgphilpott/babel

max_issues_repo_name: jgphilpott/babel

max_forks_repo_name: jgphilpott/babel


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)
function fva_preprocess(S, b, lb, ub, rxns;
        check_obj = nothing,
        verbose = true, eps = 0.0,
        batchlen = 50,
        check_obj_atol::Real = 1e-4,
        ignore = [], # rxns skip preprocess
        protect = [], # rxns skip blocking
    )

    # FVA
    _bidx = trues(size(S, 2))
    _bidx[ignore] .= false
    non_ignored = findall(_bidx)

    fvalb, fvaub = (lb, ub) .|> copy
    fvalb[non_ignored], fvaub[non_ignored] = fva(
        S, b, lb, ub, non_ignored;
        check_obj_atol, check_obj,
        batchlen, verbose
    )

    return del_blocked(S, b, fvalb, fvaub, rxns; eps, protect)
end

function fva_preprocess(metnet::MetNet;
            check_obj = nothing,
            verbose = true, eps = 0.0,
            return_blocked = false,
            batchlen = 50,
            check_obj_atol::Real = 1e-4,
            ignore = [], # rxns skip preprocess
            protect = [] # rxns skip blocking
        )
    ignore = map((r) -> rxnindex(metnet, r), ignore)
    protect = map((r) -> rxnindex(metnet, r), protect)
    check_obj = isnothing(check_obj) ? nothing : rxnindex(metnet, check_obj)

    model_fields = _extract_dense(metnet, [:S, :b, :lb, :ub])
    rxns = eachindex(metnet.rxns)
    S, b, lb, ub, rxnis, blocked = fva_preprocess(
        model_fields..., rxns;
        check_obj, verbose, eps,
        ignore, protect, batchlen,
        check_obj_atol
    );

    rxns = metnet.rxns[rxnis]
    metnet = MetNet(metnet; S, b, lb, ub, rxns)
    return return_blocked ? (metnet, blocked) : metnet
end


text:
module Fortran

using PyCall

const fparser = PyNULL()
const fparser_api = PyNULL()

function __init__()
    copy!(fparser, pyimport_conda("fparser", "fparser", "conda-forge"))
    copy!(fparser_api, pyimport_conda("fparser.api", "fparser", "conda-forge"))
end

fparse(source::String; kwargs...) = fparser_api.parse(source; kwargs...)
fparse(io::IO; kwargs...) = fparse(String(read(io)); kwargs...)

export fparse

end


text:
#!/usr/bin/env julia
using Dierckx
using Compat
using Base.Test

# 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_approx_eq(yi, [0.0, 1.0, 2.0])
@test_approx_eq(evaluate(spl, 1.5), 1.0)
@test_approx_eq(get_knots(spl), [1., 3.])
@test_approx_eq(get_coeffs(spl), [0., 4.])
@test_approx_eq_eps(get_residual(spl), 0.0, 1.e-30)

# test call on v0.4+
if VERSION >= v"0.4.0"
    @test_approx_eq(spl([1.0, 1.5, 2.0]), [0.0, 1.0, 2.0])
    @test_approx_eq(spl(1.5), 1.0)
end

# test that a copy is returned by get_knots()
knots = get_knots(spl)
knots[1] = 1000.
@test_approx_eq(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=@compat(Float64(length(x))))
desired = [0.35100374, 0.51715855, 0.87789547, 0.98719344]
actual = evaluate(spl, [0.1, 0.5, 0.9, 0.99])
@test_approx_eq_eps(actual, desired, 5e-4)

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

xp = linspace(-8.0, 13.0, 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_approx_eq(evaluate(spl, xp), xp_clip.^3)
@test_approx_eq(evaluate(spl2, xp), xp_clip.^3)

# test other bc's
spl = Spline1D(x, y; bc="extrapolate")
@test_approx_eq(evaluate(spl, xp), xp.^3)
spl = Spline1D(x, y; bc="zero")
@test_approx_eq(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 = linspace(0, 1, 70)
y = x.^3
spl = Spline1D(x, y)
xt = [0.3, 0.4, 0.5]
@test_approx_eq(derivative(spl, xt), 3xt.^2)

# test integral
x = linspace(0, 10, 70)
y = x.^2
spl = Spline1D(x, y)
@test_approx_eq(integrate(spl, 1.0, 5.0), 5.^3/3 - 1/3)

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

# -----------------------------------------------------------------------------
# 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_approx_eq(tx, [1., 3.])
@test_approx_eq(ty, [1., 3.])
@test_approx_eq_eps(get_residual(spl), 0.0, 1e-16)
@test_approx_eq(evaluate(spl, 2.0, 1.5), 2.0)
@test_approx_eq(evalgrid(spl, [1.,1.5,2.], [1.,1.5]), [0. 0.; 1. 1.; 2. 2.])

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

# In this setting, lwrk2 is too small in the default run.
x = linspace(-2, 2, 80)
y = linspace(-2, 2, 80)
z = x .+ y
spl = Spline2D(x, y, z; s=length(x))
@test_approx_eq(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.
srand(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_approx_eq(zi, ans)

# test call on v0.4+
if VERSION >= v"0.4.0"
    zi = spl(xi, yi)
    @test_approx_eq(zi, ans)
end

# 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_approx_eq(zi, ans)

println("All tests passed.")


text:
module BasisSetModule
export BasisSet, readBasisSetTX93
using ..BaseModule
using ..AtomModule

immutable PrimitiveGaussianBasisFunctionDefinition
    exponent::Float64
    prefactor::Float64
end

immutable ContractedGaussianBasisFunctionDefinition
    lQuantumNumber::LQuantumNumber
    primitives::Array{PrimitiveGaussianBasisFunctionDefinition,1}
end

immutable BasisSet
    definitions::Dict{Element,Array{ContractedGaussianBasisFunctionDefinition,1}}
end

function readBasisSetTX93(filename::AbstractString)
  basSet = BasisSet(Dict())                                 # initialize return value

  elem = Element("H")                                       # functionglobal elem
  lqn = "S"                                         # functionglobal lqn
  contractedDefinition = ContractedGaussianBasisFunctionDefinition(LQuantumNumber("S"),[])  # functionglobal contractedDefinition

  fd = open(filename)
  for line in eachline(fd)
    # comments
    if ismatch(r"^(!|\s*$)",line) continue end              # skip

    # new element
    if ismatch(r"^FOR",line)
      if(contractedDefinition.primitives != [])                 #   and not the first
    append!(basSet.definitions[elem],[contractedDefinition])                    #   finish contractedDefinition
    contractedDefinition = ContractedGaussianBasisFunctionDefinition(LQuantumNumber(lqn),[])    #   start new contractedDefinition
      end
      elem = Element(match(r"^FOR\s*(\w*)",line).captures[1])           #   update elem
      basSet.definitions[elem] = ContractedGaussianBasisFunctionDefinition[]    #   initialize elem-definitions empty

    # definition line
    else
      regmatch = match(Regex(raw"(?P<lqn>\w*)?\s*(?P<exp>"*floatregex*raw")\s*(?P<lin>"*floatregex*raw")"),line)
      # new contracted
      if regmatch[:lqn] != ""                   # if new contracted
    lqn = regmatch[:lqn]
    if(contractedDefinition.primitives != [])               #    not a new element
      append!(basSet.definitions[elem],[contractedDefinition])                  #   finish contractedDefinition
    end
      contractedDefinition = ContractedGaussianBasisFunctionDefinition(LQuantumNumber(lqn),[])  #   start new contractedDefinition
      end
      # always (add primitive to current contracted)
      exponent = float(regmatch[:exp])
      linfactor = float(regmatch[:lin])
      append!(contractedDefinition.primitives,[PrimitiveGaussianBasisFunctionDefinition(exponent,linfactor)])
    end
  end
  if(contractedDefinition.primitives != [])                 #   if still contracteds left
    append!(basSet.definitions[elem],[contractedDefinition])        #   finish contractedDefinition
  end
  close(fd)
  return basSet
end

end # module


text:
include("helper.jl")

##############################
##### Initialize Figure ######
##############################
@info("Initializing Figure...")

noto_sans = "./resources/NotoSans-Regular.ttf"
noto_sans_bold=  "./resources/NotoSans-Bold.ttf"

tickfontsize    = 28
labelfontsize   = 32

f = Figure(font=noto_sans, resolution=(1600,400))

####################################
###### Plot Differential Form ######
####################################
@info("Plotting Differential Form...")

vanderpol_arrow(u) = vanderpol(u,1,1)
ax1 = Axis(f[1,1], aspect=1,
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x₍₁₎", ylabel = "x₍₂₎",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize)
u_grid = -4:.333:4

us = zeros(length(u_grid)^2)
uv = zeros(length(u_grid)^2)
i=1
for y in u_grid
    for x in u_grid
        du = vanderpol_arrow([x,y])
        du /= sqrt(sum(abs2,du))
        us[i], uv[i] = du
        global i += 1
    end
end

us = reshape(us,length(u_grid),length(u_grid))/4
uv = reshape(uv,length(u_grid),length(u_grid))/4

arrows!(    ax1,
            u_grid,
            u_grid,
            us,
            uv)

xlims!(ax1, (-4,4))
ylims!(ax1, (-4,4))

###############################
###### Plot Trajectories ######
###############################

ax2 = Axis( f[1,2], aspect=1,
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x₍₁₎", ylabel = "x₍₂₎",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize)

# Four different initial conditions
u0 = zeros(2,4)
u0[2,:] = [-.1,.1,-4,4]

for i in 1:4
    μ = 1.0
    tspan = (0.0,20.0)

    prob = ODEProblem(vanderpol!,u0[:,i],tspan,μ)

    sol = solve(prob,saveat=0.05)

    lines!(ax2,sol, color=:black)
end
xlims!(ax2, (-4,4))
ylims!(ax2, (-4,4))

##########################################
############ Phase Plot ##################
##########################################
@info("Generating Phase Plot...")

# Determine period of limit cycle

let
  u0 = [1,1]
  μ = 1.0
  tspan = (0.0,20.0)

  prob = ODEProblem(vanderpol!,u0,tspan,μ)

  sol = solve(prob,saveat=0.01)

  # Compute the distances from the final point and get one cycle
  distance = zeros(2001-100)
  for i in 1:2001-100
      distance[i] = sum(abs2,sol[end] - sol[i])
  end
  ~,ind = findmin(distance)
  period = 2001-ind
  global positions = sol.u[end-period+1:end]
  global positions = cat(positions...,dims=2)
end

u_grid = -4:.01:4

results = zeros(length(u_grid)^2)

Threads.@threads for i in 1:length(u_grid)
  for j in 1:length(u_grid)
    y = u_grid[j]
    x = u_grid[i]
    results[j + (i-1)*length(u_grid)] = phasemap([x,y],positions)
  end
end

results = reshape(results,length(u_grid),length(u_grid))

ax3 = Axis( f[1,4],
            aspect=1,
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x₍₁₎", ylabel = "x₍₂₎",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize)

heatmap!(ax3, u_grid, u_grid, results, colormap=:phase)
contour!(ax3, u_grid, u_grid, results, levels=31, color=:black)
xlims!(ax3,(-4,4))
ylims!(ax3,(-4,4))

#####################################
###### Generate Arrow Plot ##########
#####################################
@info("Plotting Informative Measurements...")

u_grid = -4:.3333:4
out_u = zeros(length(u_grid),length(u_grid))
out_v = zeros(length(u_grid),length(u_grid))

Threads.@threads for i in 1:length(u_grid)
    for j in 1:length(u_grid)
        u₀ = [u_grid[i], u_grid[j]]
        dx = ForwardDiff.jacobian(simsystem,u₀)
        ~, Σ, V = svd(dx)
        # Force consistent arrow direction
        if u_grid[j] > 2.5 && u_grid[i] > 2.5
            test_vec = [1, 1]
        elseif u_grid[j] < -2.5 && u_grid[i] < -2.5
            test_vec = [-1, -1]
        else
            test_vec = [u_grid[j], -u_grid[i]]
        end
        if dot(test_vec,V[:,1]) < 0
            out_u[i,j] = V[1,1]
            out_v[i,j] = V[2,1]
        else
            out_u[i,j] = -V[1,1]
            out_v[i,j] = -V[2,1]
        end
    end
end

ax4 = Axis(f[1,3], aspect=1,
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="x₍₁₎", ylabel = "x₍₂₎",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize)

arrows!(    ax4,
            u_grid,
            u_grid,
            out_u./4,
            out_v./4)

xlims!(ax4, (-4,4))
ylims!(ax4, (-4,4))

##########################################
########### Panel Labels  ################
##########################################

Label(f[1,1,TopLeft()], "A", textsize = 34, padding=(0,5,5,0), halign=:right, font = noto_sans_bold)
Label(f[1,2,TopLeft()], "B", textsize = 34, padding=(0,5,5,0), halign=:right, font = noto_sans_bold)
Label(f[1,3,TopLeft()], "C", textsize = 34, padding=(0,5,5,0), halign=:right, font = noto_sans_bold)
Label(f[1,4,TopLeft()], "D", textsize = 34, padding=(0,5,5,0), halign=:right, font = noto_sans_bold)

colgap!(f.layout, 0)
rowgap!(f.layout, 0)

##########################################
############### Save Plot ################
##########################################
@info("Saving Figure...")
save("out/demo.pdf",f);


text:
using Ditherings
using Base.Test

@test true


text:
function hex_neighbors(hex::Tuple{Int,Int})
    i,j = hex
    [(i+1,j),(i,j+1),(i-1,j+1),(i-1,j),(i,j-1),(i+1,j-1)]
end

struct HexWorldMDP
    # Problem has |hexes| + 1 states, where last state is consuming.
    hexes::Vector{Tuple{Int,Int}}

    # The exact same problem as a DiscreteMDP
    mdp::DiscreteMDP

    # The special hex rewards used to construct the MDP
    special_hex_rewards::Dict{Tuple{Int,Int}, Float64}

    function HexWorldMDP(
        hexes::Vector{Tuple{Int,Int}},
        r_bump_border::Float64,
        p_intended::Float64,
        special_hex_rewards::Dict{Tuple{Int,Int}, Float64},
        γ::Float64,
        )

        nS = length(hexes) + 1 # Hexes plus one terminal state
        nA = 6 # Six directions. 1 is east, 2 is north east, 3 is north west, etc.
               # As enumerated in hex_neighbors.

        s_absorbing = nS

        T = zeros(Float64, nS, nA, nS)
        R = zeros(Float64, nS, nA)

        p_veer = (1.0 - p_intended)/2 # Odds of veering left or right.

        for s in 1 : length(hexes)
            hex = hexes[s]
            if !haskey(special_hex_rewards, hex)
                # Action taken from a normal tile
                neighbors = hex_neighbors(hex)
                for (a,neigh) in enumerate(neighbors)
                    # Indended transition.
                    s′ = findfirst(h -> h == neigh, hexes)
                    if s′ == nothing
                        # Off the map!
                        s′ = s
                        R[s,a] += r_bump_border*p_intended
                    end
                    T[s,a,s′] += p_intended

                    # Unintended veer left.
                    a_left = mod1(a+1, nA)
                    neigh_left = neighbors[a_left]
                    s′ = findfirst(h -> h == neigh_left, hexes)
                    if s′ == nothing
                        # Off the map!
                        s′ = s
                        R[s,a] += r_bump_border*p_veer
                    end
                    T[s,a,s′] += p_veer

                    # Unintended veer right.
                    a_right = mod1(a-1, nA)
                    neigh_right = neighbors[a_right]
                    s′ = findfirst(h -> h == neigh_right, hexes)
                    if s′ == nothing
                        # Off the map!
                        s′ = s
                        R[s,a] += r_bump_border*p_veer
                    end
                    T[s,a,s′] += p_veer
                end
            else
                # Action taken from an absorbing hex
                # In absorbing hex, your action automatically takes you to the absorbing state and you get the reward.
                for a in 1 : nA
                    T[s,a,s_absorbing] = 1.0
                    R[s,a] += special_hex_rewards[hex]
                end
            end
        end

        # Absorbing state stays where it is and gets no reward.
        for a in 1 : nA
            T[s_absorbing,a,s_absorbing] = 1.0
        end

        mdp = DiscreteMDP(T,R,γ)

        return new(hexes, mdp, special_hex_rewards)
    end
end

const HexWorldRBumpBorder = -1.0 # Reward for falling off hex map
const HexWorldPIntended = 0.7 # Probability of going intended direction
const HexWorldDiscountFactor = 0.9
# const HexWorld = HexWorldMDP(
#     [(0,0),(1,0),(2,0),(3,0),(0,1),(1,1),(2,1),(-1,2),
#      (0,2),(1,2),(2,2),(3,2),(4,2),(5,2),(6,2),(7,2),
#      (8,2),(4,1),(5,0),(6,0),(7,0),(7,1),(8,1),(9,0)],
#     HexWorldRBumpBorder,
#     HexWorldPIntended,
#     Dict{Tuple{Int,Int}, Float64}(
#         (0,1)=>  5.0, # left side reward
#         (2,0)=>-10.0, # left side hazard
#         (9,0)=> 10.0, # right side reward
#     ),
#     HexWorldDiscountFactor
# )

function HexWorld()
    HexWorld = HexWorldMDP(
        [(0,0),(1,0),(2,0),(3,0),(0,1),(1,1),(2,1),(-1,2),
         (0,2),(1,2),(2,2),(3,2),(4,2),(5,2),(6,2),(7,2),
         (8,2),(4,1),(5,0),(6,0),(7,0),(7,1),(8,1),(9,0)],
        HexWorldRBumpBorder,
        HexWorldPIntended,
        Dict{Tuple{Int,Int}, Float64}(
            (0,1)=>  5.0, # left side reward
            (2,0)=>-10.0, # left side hazard
            (9,0)=> 10.0, # right side reward
        ),
        HexWorldDiscountFactor
    )
    return HexWorld
end

# const StraightLineHexWorld = HexWorldMDP(
#     [(0,0),(1,0),(2,0),(3,0),(4,0),(5,0),(6,0)],
#     HexWorldRBumpBorder,
#     HexWorldPIntended,
#     Dict{Tuple{Int,Int}, Float64}(
#         (6,0)=> 10.0, # right side reward
#     ),
#     HexWorldDiscountFactor
# )

function StraightLineHexWorld()
    StraightLineHexWorld = HexWorldMDP(
        [(0,0),(1,0),(2,0),(3,0),(0,1),(1,1),(2,1),(-1,2),
         (0,2),(1,2),(2,2),(3,2),(4,2),(5,2),(6,2),(7,2),
         (8,2),(4,1),(5,0),(6,0),(7,0),(7,1),(8,1),(9,0)],
        HexWorldRBumpBorder,
        HexWorldPIntended,
        Dict{Tuple{Int,Int}, Float64}(
            (0,1)=>  5.0, # left side reward
            (2,0)=>-10.0, # left side hazard
            (9,0)=> 10.0, # right side reward
        ),
        HexWorldDiscountFactor
    )
    return StraightLineHexWorld
end

n_states(mdp::HexWorldMDP) = n_states(mdp.mdp)
n_actions(mdp::HexWorldMDP) = n_actions(mdp.mdp)
discount(mdp::HexWorldMDP) = discount(mdp.mdp)
ordered_states(mdp::HexWorldMDP) = ordered_states(mdp.mdp)
ordered_actions(mdp::HexWorldMDP) = ordered_actions(mdp.mdp)
state_index(mdp::HexWorldMDP, s::Int) = s

transition(mdp::HexWorldMDP, s::Int, a::Int) = transition(mdp.mdp, s, a)
generate_s(mdp::HexWorldMDP, s::Int, a::Int) = generate_s(mdp.mdp, s, a)
reward(mdp::HexWorldMDP, s::Int, a::Int) = reward(mdp.mdp, s, a)
generate_sr(mdp::HexWorldMDP, s::Int, a::Int) = (generate_s(mdp, s, a), reward(mdp,s,a))

generate_start_state(mdp::HexWorldMDP) = rand(1:(n_states(mdp)-1)) # non-terminal state

function hex_distance(a::Tuple{Int,Int}, b::Tuple{Int,Int})
    az = -a[1] - a[2]
    bz = -b[1] - b[2]
    return max(abs(a[1] - b[1]), abs(a[2] - b[2]), abs(az - bz))
end

function DiscreteMDP(mdp::HexWorldMDP)
    return mdp.mdp
end
function MDP(mdp::HexWorldMDP)
    return MDP(mdp.mdp)
end


text:
using TDMSReader, Gtk, Plots, StatsPlots

struct SingleRead
    TestParams
    ImpedanceData
    GainData
    CurrentTimeDomain
    VoltageTimeDomain
end

struct FreqSweepRead
    FreqVec
    ZMagVec
    ZPhaseVec
    GainData
    TestParams
end

struct StatsRead
    ImpedanceData
    ZPhaseArray
    ZMagArray
    GainData
    TestParams
end

function UnpackSingleRead()
    UnpackSingleRead(open_dialog("Pick a file"))
end

function UnpackSingleRead(FilePath)
    AllData = readtdms(FilePath)
    ImpedanceData = AllData.groups["Impedance"].props
    ImpedanceData["Ls"]

    GainData = AllData.groups["Calib_and_Gain"].props
    CurrentTimeDomain = AllData.groups["SampleTimeDomain"].channels["Current"].data
    VoltageTimeDomain = AllData.groups["SampleTimeDomain"].channels["Voltage"].data

    TestParams =  AllData.groups["Test_Parameters"].props

    return SingleRead(TestParams,ImpedanceData,GainData,CurrentTimeDomain,VoltageTimeDomain)

end

function UnpackSweepRead()
    UnpackSweepRead(open_dialog("Pick a file"))
end
function UnpackSweepRead(FilePath)
    AllData = readtdms(FilePath)
    FreqVec = AllData.groups["Freq_Sweep_Data"].channels["Freq_Vec"].data
    ZMagVec = AllData.groups["Freq_Sweep_Data"].channels["ZMag_Vec"].data
    ZPhaseVec = AllData.groups["Freq_Sweep_Data"].channels["ZPhase_Vec"].data

    GainData = AllData.groups["Calib_and_Gain"].props

    TestParams =  AllData.groups["Test_Parameters"].props

    return FreqSweepRead(FreqVec,ZMagVec,ZPhaseVec,GainData,TestParams)

end

function UnpackStatsRead()
    UnpackStatsRead(open_dialog("Pick a file"))
end

function UnpackStatsRead(FilePath)
    AllData = readtdms(FilePath)
    ImpedanceData = AllData.groups["Impedance"].props
    ZPhaseArray = AllData.groups["Data Arrays"].channels["ZPhase_Array"].data
    ZMagArray = AllData.groups["Data Arrays"].channels["ZMag_Array"].data
    GainData = AllData.groups["Calib_and_Gain"].props

    TestParams =  AllData.groups["Test_Parameters"].props

    return StatsRead(ImpedanceData,ZPhaseArray,ZMagArray,GainData,TestParams)

end

function PlotFreqSweep(DataIn::FreqSweepRead;XScaling= :identity, YScaling = :log10)

    XVec = [DataIn.FreqVec , DataIn.FreqVec]
    YVec = [DataIn.ZMagVec , DataIn.ZPhaseVec]
    plot(XVec,YVec,layout = (2,1),xscale = XScaling,label = ["|Z|" "∠Z"],legend = :best,lc=[:blue :green],size = (800,400),xlabel = "Frequency[Hz]",ylabel = ["|Z| [Ω]" "∠Z [deg.]"])
    # , yscale = [YScaling , :identity],
end

function PolarDeg2Cart(R,Θ)
    RealPart = real.(R.*exp.(im*π/180 .*Θ))
    ImagPart = imag.(R.*exp.(im*π/180 .*Θ))
    return RealPart, ImagPart
end
function PolarDeg2Complex(R,Θ)

    return (R.*exp.(im*π/180 .*Θ))
end

function TmpUnpackStatsRead()
    AllData = readtdms((open_dialog("Pick a file")))
    ImpedanceData = AllData.groups["Impedance"].props
    println(ImpedanceData)

end


text:
# define common PETSc constants, independant from PETSC_ARCH
# -> directly copied from https://github.com/JuliaParallel/PETSc.jl.git

const PetscErrorCode = Cint

const PETSC_DEFAULT = -2
const PETSC_DECIDE  = -1

"""
    Macro to automatically export all items of an enum
"""
macro exported_enum(name, args...)
    esc(quote
        @enum($name, $(args...))
        export $name
        $([:(export $arg) for arg in args]...)
        end)
end

@enum PetscBool PETSC_FALSE PETSC_TRUE

@enum PetscDataType begin
    PETSC_DATATYPE_UNKNOWN  = 0
    PETSC_DOUBLE            = 1
    PETSC_COMPLEX           = 2
    PETSC_LONG              = 3
    PETSC_SHORT             = 4
    PETSC_FLOAT             = 5
    PETSC_CHAR              = 6
    PETSC_BIT_LOGICAL       = 7
    PETSC_ENUM              = 8
    PETSC_BOOL              = 9
    PETSC_FLOAT128          = 10
    PETSC_OBJECT            = 11
    PETSC_FUNCTION          = 12
    PETSC_STRING            = 13
    PETSC___FP16            = 14
    PETSC_STRUCT            = 15
    PETSC_INT               = 16
    PETSC_INT64             = 17
end

const Petsc64bitInt  = Int64
const PetscLogDouble = Cdouble

@enum InsertMode begin
    NOT_SET_VALUE
    INSERT_VALUES
    ADD_VALUES
    MAX_VALUES
    MIN_VALUES
    INSERT_ALL_VALUES
    ADD_ALL_VALUES
    INSERT_BC_VALUES
    ADD_BC_VALUES
end

@enum NormType begin
    NORM_1 = 0
    NORM_2 = 1
    NORM_FROBENIUS = 2
    NORM_INFINITY = 3
    NORM_1_AND_2 = 4
end

@enum MatAssemblyType begin
    MAT_FLUSH_ASSEMBLY = 1
    MAT_FINAL_ASSEMBLY = 0
end

@enum MatFactorType begin
    MAT_FACTOR_NONE     = 0
    MAT_FACTOR_LU       = 1
    MAT_FACTOR_CHOLESKY = 2
    MAT_FACTOR_ILU      = 3
    MAT_FACTOR_ICC      = 4
    MAT_FACTOR_ILUDT    = 5
end

@enum PetscViewerFormat begin
    PETSC_VIEWER_DEFAULT
    PETSC_VIEWER_ASCII_MATLAB
    PETSC_VIEWER_ASCII_MATHEMATICA
    PETSC_VIEWER_ASCII_IMPL
    PETSC_VIEWER_ASCII_INFO
    PETSC_VIEWER_ASCII_INFO_DETAIL
    PETSC_VIEWER_ASCII_COMMON
    PETSC_VIEWER_ASCII_SYMMODU
    PETSC_VIEWER_ASCII_INDEX
    PETSC_VIEWER_ASCII_DENSE
    PETSC_VIEWER_ASCII_MATRIXMARKET
    PETSC_VIEWER_ASCII_VTK
    PETSC_VIEWER_ASCII_VTK_CELL
    PETSC_VIEWER_ASCII_VTK_COORDS
    PETSC_VIEWER_ASCII_PCICE
    PETSC_VIEWER_ASCII_PYTHON
    PETSC_VIEWER_ASCII_FACTOR_INFO
    PETSC_VIEWER_ASCII_LATEX
    PETSC_VIEWER_ASCII_XML
    PETSC_VIEWER_ASCII_GLVIS
    PETSC_VIEWER_ASCII_CSV
    PETSC_VIEWER_DRAW_BASIC
    PETSC_VIEWER_DRAW_LG
    PETSC_VIEWER_DRAW_LG_XRANGE
    PETSC_VIEWER_DRAW_CONTOUR
    PETSC_VIEWER_DRAW_PORTS
    PETSC_VIEWER_VTK_VTS
    PETSC_VIEWER_VTK_VTR
    PETSC_VIEWER_VTK_VTU
    PETSC_VIEWER_BINARY_MATLAB
    PETSC_VIEWER_NATIVE
    PETSC_VIEWER_HDF5_PETSC
    PETSC_VIEWER_HDF5_VIZ
    PETSC_VIEWER_HDF5_XDMF
    PETSC_VIEWER_HDF5_MAT
    PETSC_VIEWER_NOFORMAT
    PETSC_VIEWER_LOAD_BALANCE
end

@enum MatOperation begin
    MATOP_SET_VALUES=0
    MATOP_GET_ROW=1
    MATOP_RESTORE_ROW=2
    MATOP_MULT=3
    MATOP_MULT_ADD=4
    MATOP_MULT_TRANSPOSE=5
    MATOP_MULT_TRANSPOSE_ADD=6
    MATOP_SOLVE=7
    MATOP_SOLVE_ADD=8
    MATOP_SOLVE_TRANSPOSE=9
    MATOP_SOLVE_TRANSPOSE_ADD=10
end


text:
#
# Test the specific vector trnasport along implementations that do iterative transport,
# also the Schild and pole special cases
#
using ManifoldsBase, Test

struct NonDefaultEuclidean <: Manifold{ManifoldsBase.ℝ} end
ManifoldsBase.log!(::NonDefaultEuclidean, v, x, y) = (v .= y .- x)
ManifoldsBase.exp!(::NonDefaultEuclidean, y, x, v) = (y .= x .+ v)
function ManifoldsBase.vector_transport_to!(
    ::NonDefaultEuclidean,
    vto,
    x,
    v,
    y,
    ::ParallelTransport,
)
    return copyto!(vto, v)
end

@testset "vector_transport_along" begin
    M = NonDefaultEuclidean()
    types = [Vector{Float64}, Vector{Float32}]
    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])
            v2 = log(M, pts[1], pts[3])
            c = [0.5 * (pts[1] + pts[2]), pts[2], 0.5 * (pts[2] + pts[3]), pts[3]]
            @test vector_transport_along(M, pts[1], v2, c, SchildsLadderTransport()) == v2
            @test vector_transport_along(M, pts[1], v2, c, PoleLadderTransport()) == v2
            @test vector_transport_along(M, pts[1], v2, c, ParallelTransport()) == v2
            @test vector_transport_along(M, pts[1], v2, [], SchildsLadderTransport()) == v2
            @test vector_transport_along(M, pts[1], v2, [], PoleLadderTransport()) == v2
            @test vector_transport_along(M, pts[1], v2, [], ParallelTransport()) == v2
            # check mutating ones with defaults
            p = allocate(pts[1])
            ManifoldsBase.pole_ladder!(M, p, pts[1], pts[2], pts[3])
            # -log_p3 p == log_p1 p2
            @test isapprox(M, -log(M, pts[3], p), log(M, pts[1], pts[2]))
            ManifoldsBase.schilds_ladder!(M, p, pts[1], pts[2], pts[3])
            @test isapprox(M, log(M, pts[3], p), log(M, pts[1], pts[2]))
        end
    end
end


text:
include("../../configs.jl")
include("../core/metrics.jl")

function encode_gt_and_get_indices(gt, priors, losses, pos_thold, neg_thold)
    variances = size(priors, 2) == 102300 ? [1, 1] : [0.2, 0.1]
    iou_vals = iou(gt[1:4,:], _to_min_max_form(priors))
    # selecting positive prior boxes
    pos_pairs = findall(iou_vals .>= pos_thold)
    pos_gt = zeros(15, length(pos_pairs))
    gt_idx = getindex.(pos_pairs, [1 2])[:, 2]
    prior_idx = getindex.(pos_pairs, [1 2])[:, 1]

    # enlarging gt matrix to match selected anchors
    pos_gt .= gt[:,gt_idx]
    num_poses = length(prior_idx)
    if num_poses == 0
        # if no positive anchor boxes are found, then no loss will be calculated
        return nothing, nothing, nothing
    end

    #selecting negative prior boxes
    max_prior_vals, max_prior_idx = findmax(iou_vals; dims=2)
    neg_indices = getindex.(findall(max_prior_vals .<= neg_thold), [1 2])[:,1]
    neg_indices = neg_indices[sortperm(losses[neg_indices])]

    neg_cnt = ohem_ratio * num_poses
    neg_cnt = neg_cnt < 1 ? 1 : neg_cnt
    neg_indices = neg_indices[1:neg_cnt]

    # gt bbox conversion
    selected_priors = priors[:,prior_idx]

    pos_gt[1:4,:] = _to_center_length_form(pos_gt[1:4,:])
    pos_gt[3:4,:] = log.(pos_gt[3:4,:] ./ selected_priors[3:4,:]) ./ variances[1]

    pos_gt[1:2,:] = (pos_gt[1:2,:] .- selected_priors[1:2,:]) ./ (variances[2] .* selected_priors[3:4,:])
    pos_gt[5:6,:] = (pos_gt[5:6,:] .- selected_priors[1:2,:]) ./ (variances[2] .* selected_priors[3:4,:])
    pos_gt[7:8,:] = (pos_gt[7:8,:] .- selected_priors[1:2,:]) ./ (variances[2] .* selected_priors[3:4,:])
    pos_gt[9:10,:] = (pos_gt[9:10,:] .- selected_priors[1:2,:]) ./ (variances[2] .* selected_priors[3:4,:])
    pos_gt[11:12,:] = (pos_gt[11:12,:] .- selected_priors[1:2,:]) ./ (variances[2] .* selected_priors[3:4,:])
    pos_gt[13:14,:] = (pos_gt[13:14,:] .- selected_priors[1:2,:]) ./ (variances[2] .* selected_priors[3:4,:])

    # print("Pos_indices: ", prior_idx, "&& Neg_indices: ", neg_indices, "\n")
    return pos_gt, prior_idx, neg_indices
end

function nms(scores, points; thold=0.4)
    x1 = points[1,:]; y1 = points[2,:]; x2 = points[3,:]; y2 = points[4,:];

    areas = (x2 - x1) .* (y2 - y1)
    order = sortperm(scores, rev=true)
    keep = []

    for i in order
        if keep == []
            push!(keep, i)
        else
            iou_vals = iou(points[:,i:i], points[:,keep])
            val, idx = findmax(iou_vals, dims=1)
            if val[1] <= thold
                push!(keep, i)
            end
        end
    end
    return keep
end

"""
Conversion from:
(center_x, center_y, width, height) --> (min_x, min_y, max_x, max_y)
"""
function _to_min_max_form(boxes)
    init_shape = size(boxes)
    boxes = reshape(boxes, (size(boxes)[1], prod(size(boxes)[2:end]))) # converting to 2D
    half_lens = boxes[3:4,:] ./ 2
    return reshape(cat(boxes[1:2,:] .- half_lens, boxes[1:2,:] .+ half_lens, dims=1), init_shape)
end

"""
Conversion from:
(min_x, min_y, max_x, max_y) --> (center_x, center_y, width, height)
"""
function _to_center_length_form(boxes)
    init_shape = size(boxes)
    boxes = reshape(boxes, (size(boxes)[1], prod(size(boxes)[2:end]))) # converting to 2D
    lengths = boxes[3:4,:] .- boxes[1:2,:]
    return reshape(cat(boxes[1:2,:] .+ (lengths ./ 2), lengths, dims=1), init_shape)
end

"""
Returns the anchor boxes with their center_x, center_y, width, height information.
"""
function _get_priorboxes(num_anchors, anchor_info, img_size)
    feature_maps = [Int(ceil(img_size / scale["stride"])) for scale in anchor_info]
    num_proposals = num_anchors * sum([i*i for i in feature_maps])
    anchors = zeros(4, num_proposals)

    counter = 1
    for (idx, f) in enumerate(feature_maps)
        scaler = anchor_info[idx]["stride"]
        for h in 1:f
            cy = (h - 0.5) * scaler
            for w in 1:f
                cx = (w - 0.5) * scaler
                for s in anchor_info[idx]["anchors"]
                    anchors[:,counter] = [cx cy s s]
                    counter += 1
                end
            end
        end
    end
    return anchors
end

"""
Decoder functions:
The main motivation is that the network does not directly predict the proposals but
the combination of prediction and prior anchor box style and rescaling gives the actual
bounding box and landmark coordinations.
"""
function decode_points(bboxes, landmarks, priors)
    if length(size(priors)) == 2
        priors = reshape(priors, (size(priors)..., 1))
    end
    decoded_bboxes = _decode_bboxes(bboxes, priors)
    decoded_landmarks = _decode_landmarks(landmarks, priors)
    return decoded_bboxes, decoded_landmarks
end

function _decode_bboxes(bbox, priors)
    variances = size(priors, 2) > 100000 ? [1, 1] : [0.2, 0.1]
    if length(size(priors)) == 2
        priors = reshape(priors, (size(priors)..., 1))
    end
    centers = priors[1:2,:,:] .+ bbox[1:2,:,:] .* variances[2] .* priors[3:4,:,:]
    lengths = exp.(bbox[3:4,:,:] .* variances[1]) .* priors[3:4,:,:]
    return cat(centers, lengths, dims=1)
end

function _decode_landmarks(landmarks, priors)
    variances = size(priors, 2) > 100000 ? [1, 1] : [0.2, 0.1]
    if length(size(priors)) == 2
        priors = reshape(priors, (size(priors)..., 1))
    end
    lm1 = priors[1:2,:,:] .+ landmarks[1:2,:,:]  .* variances[2] .* priors[3:4,:,:]
    lm2 = priors[1:2,:,:] .+ landmarks[3:4,:,:]  .* variances[2] .* priors[3:4,:,:]
    lm3 = priors[1:2,:,:] .+ landmarks[5:6,:,:]  .* variances[2] .* priors[3:4,:,:]
    lm4 = priors[1:2,:,:] .+ landmarks[7:8,:,:]  .* variances[2] .* priors[3:4,:,:]
    lm5 = priors[1:2,:,:] .+ landmarks[9:10,:,:] .* variances[2] .* priors[3:4,:,:]
    return cat(lm1, lm2, lm3, lm4, lm5, dims=1)
end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
#Simple parallel reduction for random eigenvalue computation

N =1024
t1=100
include("randeig.jl")
nps = Float64[]
gflops = Float64[]
for np in [1,2,4,8,16,32,64]#1:64#Sys.CPU_CORES
    ps = addprocs(np==1 ? 1 : np-nworkers())
    @sync for p in ps
        @spawnat p include("randeig.jl")
    end
    t = @elapsed randeig(t1*np, N)
    gflop = t1*np*4/3*N^3*1e-9/t
    println(nworkers(), '\t', t, '\t', gflop)
    push!(nps, nworkers())
    push!(gflops, gflop)
end

using Plots
plot(nps, gflops, c = :blue, label = "")
scatter!(nps, gflops, c = :blue, label="")
plot!([extrema(nps)...], [gflops[1], nps[end]*gflops[1], c = :grey, label="Ideal"])
png("../plots/parfor")


text:
using Plots

L=2
function φ(x,y;V=10)
    L=2
    s=0
    for n=1:2:100000
        s=s.+map(z  -> sin(n *π *z/L) ,x).*map( z -> exp(-n *π *z/L)./n,y)
    end
    4s*V/ π
end

Space = zeros(100,100)
x=0:2/(99):2;y=x

length(x)
for i ∈ 1:100  , j ∈ 1:100
    Space[i,j]=φ(x[i],y[j])
end
plot(x,y,Space,linetype=:contourf,title = "Analitica")

max_stars_repo_name: munozariasjm/Laplace-and-Julia

max_issues_repo_name: munozariasjm/Laplace-and-Julia

max_forks_repo_name: munozariasjm/Laplace-and-Julia


text:
# Stochastic Chebyshev Polynomial method
# Based on APPROXIMATING SPECTRAL SUMS OF LARGE-SCALE MATRICES USING STOCHASTIC CHEBYSHEV APPROXIMATIONS By Han, Malioutov, Avron and Shin

export ChebyHutchSpace, chebyhutch, chebydiagonal, lczeigen

using LinearAlgebra
using Parameters
using Statistics

# Lanczos iteration for finding extremal eigenvalues of the matrix
function lczeigen(A, fn, dfn)
    # Estimate eigmax and eigmin for Chebyshev bounds
    mval = Int64(ceil(log(ϵ/(1.648 * sqrt(size(A, 1))))/(-2 * sqrt(ξ))))
    w = SLQWorkspace(A, fn = fn, dfn = dfn, m = mval)
    dfn(w.v)
    w.v .= w.v ./ norm(w.v)
    lcz(w)
    λₘ = eigmax(w.T)
    λ₁ = eigmin(w.T)
    return λ₁, λₘ
end

# Chebyshev polynomial method helper functions
𝓍(k, n) = cos((π * (k + 0.5))/(n+1))

function T(j, x, Tvals)
    if j == 0
        return 1
    end
    if j == 1
        return x
    end
    if haskey(Tvals, j)
        return Tvals[j]
    else
        Tvals[j] = (2 * x * T(j-1, x, Tvals)) - T(j-2, x, Tvals)
        return Tvals[j]
    end
end

function coeff(j, n, a, b, fn)
    fs = zero(eltype(a))
    for k in 0:n
        x = 𝓍(k, n)
        Tvals = Dict{Int, Float64}()
        fs = fs + fn((((b-a)/2) * x) + (b+a)/2) * T(j, x, Tvals)
    end
    if j == 0
        return (1/(n+1)) * fs
    end
    return (2/(n+1)) * fs
end

mutable struct ChebyHutchSpace{elt, TM, FN<:Function, FN2<:Function, TA<:AbstractArray{elt, 1}, TV<:AbstractVecOrMat{elt}, I<:Int64}
    A::TM
    a::elt
    b::elt
    fn::FN
    dfn::FN2
    C::TA
    w₀::TV
    w₁::TV
    w₂::TV
    v::TV
    u::TV
    m::I
    n::I
end
"""
    ChebyHutchSpace(A::AbstractMatrix, a::Number, b::Number; fn::Function=invfun, dfn::Function=rademacherDistribution!, m = 4, n = 6)

Create Chebyshev-Hutchinson Workspace for chebyhutch or chebydiagonal methods.

# Arguments
- `A` : Symmetric Positive Definite Matrix
- `a` : Bound for minimum eigenvalue of A
- `b` : Bound for maximum eigenvalue of A
- `fn` : Function to appy. By default uses inverse function.
- `dfn` : Distribution function that returns a vector v. By default uses rademacherDistribution!
- `m` : Iteration number, increase this for precision. By default m = 4
- `n` : Polynomial degree, increase this for accuracy. By default n = 6
"""
function ChebyHutchSpace(A::AbstractMatrix, a::Number, b::Number; fn::Function=invfun, dfn::Function=rademacherDistribution!, m = 4, n = 6)
    elt = eltype(A)
    s = size(A, 1)
    C = elt[]
    v = Matrix{elt}(undef, s, m)
    w₀ = similar(v)
    w₁ = similar(v)
    w₂ = similar(v)
    u = similar(v)
    return ChebyHutchSpace(A, a, b, fn, dfn, C, w₀, w₁, w₂, v, u, m, n)
end

function chebypm(w::ChebyHutchSpace)
    @unpack A, a, b, C, fn, dfn, v, u, w₀, w₁, w₂, m, n = w
    tr = zero(eltype(A))
    for j in 0:n
        push!(C, coeff(j, n, a, b, fn))
    end
    dfn(v)
    w₀ .= v
    mul!(w₁, A, v)
    #rmul!(w₁, 2/(b-a))
    w₁ .= 2/(b-a) .* w₁ .- (((b+a)/(b-a)) .* v)
    u .= (C[1] .* w₀) .+ (C[2] .* w₁)
    for j in 2:n
        mul!(w₂, A, w₁)
        #rmul!(w₂, 4/(b-a))
        w₂ .= 4/(b-a) .* w₂ .- ((2(b+a)/(b-a)) .* w₁) .- w₀
        u .= u .+ (C[j+1] .* w₂)
        w₀ .= w₁
        w₁ .= w₂
    end
    return v, u, m
end
"""
    chebyhutch(w::ChebyHutchSpace)
    chebyhutch(A::AbstractMatrix, m::Integer, n::Integer; a::Float64, b::Float64, fn::Function=invfun, dfn::Function=rademacherDistribution!)
    chebyhutch(A::AbstractMatrix; fn::Function=invfun, dfn::Function=rademacherDistribution!)

Chebyshev-Hutchinson to estimate tr(fn(A)), for given matrix A and an analytic function fn.

# Arguments
- `A` : Symmetric Positive Definite Matrix
- `m` : Iteration number, increase this for precision. By default m = 4
- `n` : Polynomial degree, increase this for accuracy. By default n = 6
- `fn` : Function to appy. By default uses inverse function.
- `dfn` : Distribution function that returns a vector v. By default uses rademacherDistribution!
"""
function chebyhutch(w::ChebyHutchSpace)
    v, u, m = chebypm(w)
    return dot(v, u) / m
end

function chebyhutch(A::AbstractMatrix, m::Integer, n::Integer; a::Float64=-1.0, b::Float64=-1.0, fn::Function=invfun, dfn::Function=rademacherDistribution!)
    # calculate extremal eigenvals
    if a == -1.0 || b == -1.0
        λ₁, λₘ = lczeigen(A, fn, dfn)
    end
    if a != -1.0
        λ₁ = a
    end
    if b != -1.0
        λₘ = b
    end

    wx = ChebyHutchSpace(A, λₘ, λ₁, fn=fn, dfn=dfn, m = m, n = n)
    return chebyhutch(wx)
end

function chebyhutch(A::AbstractMatrix; m = -1, n = -1, fn::Function=invfun, dfn::Function=rademacherDistribution!)
    # calculate extremal eigenvals
    λ₁, λₘ = lczeigen(A, fn, dfn)

    # calculate values of m and n
    # these bounds are for theoretical purposes only
    κ = λₘ/λ₁
    ρ = sqrt(2 * κ - 1) - 1
    if m == -1
        m = Int64(ceil(54 * (ϵ)^(-2) * log(2/ξ)/16))
    end
    if n == -1
        nVal = Int64(ceil((log(8/ϵ)*ρ*κ)/(log((2/ρ) + 1))/16))
    end

    wx = ChebyHutchSpace(A, λₘ, λ₁, fn=fn, dfn=dfn, m = m, n = n)
    return chebyhutch(wx)
end
"""
    chebydiagonal(w::ChebyHutchSpace)
    chebydiagonal(A::AbstractMatrix, m::Integer, n::Integer; fn::Function=invfun, dfn::Function=rademacherDistribution!)

Chebyshev-Hutchinson to estimate diagonal elements of fn(A), for given matrix A and an analytic function fn.

# Arguments
- `A` : Symmetric Positive Definite Matrix
- `m` : Iteration number, increase this for precision. By default m = 4
- `n` : Polynomial degree, increase this for accuracy. By default n = 6
- `fn` : Function to appy. By default uses inverse function.
- `dfn` : Distribution function that returns a vector v. By default uses rademacherDistribution!
"""
function chebydiagonal(w::ChebyHutchSpace)
    v, u, m = chebypm(w)
    return vec(mean(v .* u, dims=2))
end

function chebydiagonal(A, m, n; fn::Function=invfun, dfn::Function=rademacherDistribution!)
    # calculate extremal eigenvals
    λ₁, λₘ = lczeigen(A, fn, dfn)
    wx = ChebyHutchSpace(A, λₘ, λ₁, fn=fn, dfn=dfn, m = m, n = n)
    return chebydiagonal(wx)
end


text:
## The package 'IntervalCensored' is Available on https://github.com/abikoushi/IntervalCensored.jl
using Distributions
using Random
using Plots
using StatsPlots
using QuadGK
using DataFrames
#using RCall
#using Optim
#using ForwardDiff
#using StatsFuns
#using SpecialFunctions
#using LinearAlgebra
#using ProgressMeter
]add "./projects/IntervalCensored"
#]add "https://github.com/abikoushi/IntervalCensored.jl"
using IntervalCensored

function calclp_dic_s(rng, d, EL, ER, SL, SR)
    ll = zero(EL[1])
    for i in eachindex(EL)
        if isfinite(EL[i])
            u = rand(rng)
            E = EL[i] + (ER[i]-EL[i])*u
            ys = rand(rng, truncated(d, SL[i]-E, SR[i]-E))
            ll += logpdf(d, ys)
        end
    end
    return -ll
end

function calclp_dic2(d, EL, ER, SL, SR)
    ll = zero(EL[1])
    logmu = log(mean(d))
    for i in eachindex(EL)
        if ER[i] < SL[i]
            ll += logmu+log(eqcdf(d,SR[i]-ER[i])-eqcdf(d,SL[i]-ER[i])-(eqcdf(d,SR[i]-EL[i])-eqcdf(d,SL[i]-EL[i]))) - log(ER[i] - EL[i]) - log(SR[i] - SL[i])
        elseif ER[i] < SR[i]
            ll += logmu+log(1-eqcdf(d,SL[i]-ER[i])-(eqcdf(d,SR[i]-EL[i])-eqcdf(d,SL[i]-EL[i]))) - log(ER[i] - EL[i])- log(SR[i] - SL[i])
        elseif SR[i] <= ER[i]
            ll += logmu+log(eqcdf(d,SR[i]-EL[i]) - eqcdf(d,SL[i]-EL[i]))- log(ER[i] - EL[i])- log(SR[i] - SL[i])
        end
    end
    return -ll
end

function calclp_ic_s(rng, d, EL, ER, S)
    ll = zero(EL[1])
    for i in eachindex(EL)
        if isfinite(EL[i])
            u = rand(rng)
            E = EL[i] + (ER[i]-EL[i])*u
            ll += logpdf(d,S[i]-E)
        end
    end
    return -ll
end

function simentropy(td, n, iter)
    rng = MersenneTwister()
    lp1 = zeros(iter)
    lp2 = zeros(iter)
    for i in 1:iter
        dat = make_dic(rng, td, n)
        lp1[i] = calclp_dic_s(rng, td, dat[1], dat[2], dat[3], dat[4])/n
        lp2[i] = calclp_dic2(td, dat[1], dat[2], dat[3], dat[4])/n
    end
    return lp1, lp2
end

function simentropy_ic(td, n, iter)
    rng = MersenneTwister()
    lp1 = zeros(iter)
    lp2 = zeros(iter)
    for i in 1:iter
        dat = make_ic(rng, td, n)
        lp1[i] = calclp_ic_s(rng, td, dat[1], dat[2], dat[3])/n
        lp2[i] = calclp_ic(td, dat[1], dat[2], dat[3])/n
    end
    return lp1, lp2
end

lp1, lp2 = simentropy_ic(Weibull(1.5,7), 100, 1000)
ms = [mean(lp1), mean(lp2)]
ss = [std(lp1), std(lp2)]
df = stack(DataFrame(H1=lp1, H2=lp2))
@df df violin(:variable, :value, fill="white", legend=false, tick_direction=:out, xtickfontsize=12, trim=false)
scatter!(["H1","H2"], ms[1:2], yerror = ss[1:2], color="black", ms=6)
Plots.abline!(0, entropy(Weibull(1.5, 7)), ls=:dash, color="black")
mean(lp1)-entropy(Weibull(2,7))
histogram(lp1,legend=false, color = "lightgrey")
vline!([entropy(Weibull(2,7))], linewidth=3, color = "royalblue")
vline!([mean(lp1)], linewidth=3, color = "tomato")

function sim_dic(td, md, N, iter, seed)
    rng = MersenneTwister()
    aic1 = zeros(iter)
    aic2 = zeros(iter)
    ge = zeros(iter)
    K = length(params(md))
    theta = zeros(iter,K)
    for i in 1:iter
        dat = make_dic(rng, td, N)
        fit = MCEMdic(rng, md, 10, dat[1], dat[2], dat[3], dat[4])
        ge[i] = quadgk(x -> -logpdf(fit[1],x)*pdf(td,x), 0, Inf)[1]
        aic1[i] = (calclp_dic(rng, fit[1], dat[1], dat[2], dat[3], dat[4]) + K)/N
        aic2[i] = fit[2][end]+K/N
        theta[i,:] .= params(fit[1])
    end
    return aic1, aic2, ge, theta
end

@time simout_dic = sim_dic(Weibull(0.8,7), Weibull(1,7), 100, 1000, 1234)
ms = [mean(simout_dic[1]-simout_dic[3]), mean(simout_dic[2]-simout_dic[3])]
ss = [std(simout_dic[1]-simout_dic[3]), std(simout_dic[2]-simout_dic[3])]
df = stack(DataFrame(AIC1=simout_dic[1]-simout_dic[3], AIC2=simout_dic[2]-simout_dic[3]))
@df df violin(:variable, :value, fill="white", legend=false, tick_direction=:out, xtickfontsize=12, trim=false)
scatter!(["AIC1","AIC2"], ms[1:2], yerror = ss[1:2], color="black", ms=6)
Plots.abline!(0, 0, ls=:dash, color="black")

function sim_ic(td, md, N, iter, seed)
    rng = MersenneTwister(seed)
    aic1 = zeros(iter)
    aic2 = zeros(iter)
    ge = zeros(iter)
    K = length(params(md))
    theta = zeros(iter,K)
    for i in 1:iter
        dat = make_ic(rng, td, N)
        fit = MCEMic(rng, md, 100, dat[1], dat[2], dat[3])
        ge[i] = quadgk(x -> -logpdf(fit[1],x)*pdf(td,x), 0, Inf)[1]
        aic1[i] = (calclp_ic(fit[1], dat[1], dat[2], dat[3]) + K)/N
        aic2[i] = fit[2][end]+K/N
        theta[i,:] .= params(fit[1])
    end
    return aic1, aic2, ge, theta
end

function sim_ic(td, md, N, p, iter, seed)
    rng = MersenneTwister(seed)
    aic1 = zeros(iter)
    aic2 = zeros(iter)
    ge = zeros(iter)
    K = length(params(md))
    theta = zeros(iter,K)
    for i in 1:iter
        dat = make_ic(rng, td, N, p)
        fit = MCEMic(rng, md, 100, dat[1], dat[2], dat[3])
        ge[i] = quadgk(x -> -logpdf(fit[1],x)*pdf(td, x), 0, Inf)[1]
        aic1[i] = (calclp_ic(fit[1], dat[1], dat[2], dat[3]) + K)/N
        aic2[i] = fit[2][end]+K/N
        theta[i,:] .= params(fit[1])
    end
    return aic1, aic2, ge, theta
end

@time simout_ic = sim_ic(Weibull(1.5, 3), Weibull(2, 3), 100, 500, 123)
ms = [mean(simout_ic[1]-simout_ic[3]), mean(simout_ic[2]-simout_ic[3])]
ss = [std(simout_ic[1]-simout_ic[3]), std(simout_ic[2]-simout_ic[3])]
df = stack(DataFrame(AIC1=simout_ic[1]-simout_ic[3], AIC2=simout_ic[2]-simout_ic[3]))

@df df violin(:variable, :value, fill="white", legend=false, tick_direction=:out, xtickfontsize=12, trim=false)
scatter!(["AIC1","AIC2"], ms[1:2], yerror = ss[1:2], color="black", ms=6)
Plots.abline!(0, 0, ls=:dash, color="black")

@time simout_icm = sim_ic(Gamma(1.5, 2.0), Gamma(1.0, 2.0), 100, 0.8, 500, 1111)
ms = [mean(simout_icm[1]-simout_icm[3]), mean(simout_icm[2]-simout_icm[3])]
ss = [std(simout_icm[1]-simout_icm[3]), std(simout_icm[2]-simout_icm[3])]
df2 = stack(DataFrame(AIC1=simout_icm[1]-simout_icm[3], AIC2=simout_icm[2]-simout_icm[3]))

@df df2 violin(:variable, :value, fill="white", legend=false, tick_direction=:out, xtickfontsize=12, trim=false)
scatter!(["AIC1","AIC2"], ms[1:2], yerror = ss[1:2], color="black", ms=6)
Plots.abline!(0, 0, ls=:dash, color="black")

########
function sim_icrt(td, md, N, tmax, iter, seed)
    rng = MersenneTwister(seed)
    aic1 = zeros(iter)
    aic2 = zeros(iter)
    ge = zeros(iter)
    K = length(params(md))
    theta = zeros(iter,K)
    for i in 1:iter
        dat = make_icrt(rng, td, tmax, N)
        fit = MCEMicrt(rng, md, 100, dat[1], dat[2], dat[3], tmax, 1)
        ge[i] = quadgk(x -> -logpdf(td,x)*pdf(td,x),0,Inf)[1]
        n = length(dat[1])
        aic1[i] = (calclp_icrt(fit[1] ,dat[1], dat[2], dat[3], tmax) + K)/n
        aic2[i] = fit[2][end]+K/n
        theta[i,:] .= params(fit[1])
    end
    return aic1, aic2, ge, theta
end
Tmax = 30
@time simout_icrt = sim_icrt(Weibull(1.5, 4), Weibull(1.5, 4), 100, Tmax, 500, 101)
rng = MersenneTwister()
dat = make_icrt(rng, Weibull(1.5, 4), Tmax, 100)
all(dat[2]-dat[1].>0)
plot(x -> (eqcdf(Gamma(1.5, 4),x-2)-eqcdf(Gamma(1.5, 4),x)),0,9,legend=false)
plot(x -> (ccdf(Gamma(1.5, 4),x-2)-ccdf(Gamma(1.5, 4),x-1)),0,20,legend=false)

f(t) = quadgk(x -> (ccdf(Gamma(1.5, 4),x-3)-cdf(Gamma(1.5, 4),x-1))/2,3,t)[1]
g(t) = mean(Gamma(1.5,4))*(-eqcdf(Gamma(1.5, 4),max(0,t-3))+eqcdf(Gamma(1.5, 4),max(t-1,0)))/2
f(5), g(5)
min(t-1,1)

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


text:
using Test
using LinearAlgebra
using Ferrite
using NonconvexIpopt

using TopOpt
using Arpack

fea_ins_dir = joinpath(@__DIR__, "instances", "fea_examples");
gm_ins_dir = joinpath(@__DIR__, "instances", "ground_meshes");

# @testset "Tim buckling log-barrier $problem_dim" for problem_dim in ["2d"] # , "3d"
#     file_name = "tim_$(problem_dim).json"
#     problem_file = joinpath(gm_ins_dir, file_name)

#     mats = TrussFEAMaterial(1.0, 0.3);
#     crossecs = TrussFEACrossSec(800.0);

#     node_points, elements, _, _ , fixities, load_cases = load_truss_json(problem_file)
#     ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
#     loads = load_cases["0"]

#     problem = TrussProblem(Val{:Linear}, node_points, elements, loads, fixities, mats, crossecs);

#     xmin = 0.0001 # minimum density
#     p = 4.0 # penalty
#     V = 0.5 # maximum volume fraction
#     x0 = fill(1.0, ncells) # initial design

#     solver = FEASolver(Direct, problem);

#     # # * Before optimization, check initial design stability
#     # solver.vars = x0
#     # solver()
#     # K, G = buckling(problem, solver.globalinfo, solver.elementinfo; u=solver.u);
#     # @test isfinite(logdet(cholesky(K+G)))
#     # sparse_eigvals, buckmodes = eigs(-G,K, nev=1, which=:LR)
#     # smallest_pos_eigval = 1/sparse_eigvals[1]
#     # @test smallest_pos_eigval >= 1.0

#     Nonconvex.NonconvexCore.show_residuals[] = true

#     comp = TopOpt.Compliance(problem, solver)
#     # TODO "manual" interior point loop, adjusting the c value every iter
#     for c in [0.1] # 10:-0.1:0.1
#         function obj(x)
#             solver.vars = x
#             # trigger assembly
#             solver()
#             K, G = buckling(problem, solver.globalinfo, solver.elementinfo);
#             # minimize compliance
#             return comp(x) - c*logdet(cholesky(Array(K+G)))
#         end
#         function constr(x)
#             # volume fraction constraint
#             return sum(x) / length(x) - V
#         end

#         m = Model(obj)
#         addvar!(m, zeros(length(x0)), ones(length(x0)))
#         Nonconvex.add_ineq_constraint!(m, constr)

#         options = MMAOptions(
#             maxiter=1000, tol = Tolerance(kkt = 1e-4, f = 1e-4),
#         )
#         TopOpt.setpenalty!(solver, p)
#         r = Nonconvex.optimize(
#             m, MMA87(dualoptimizer = ConjugateGradient()),
#             x0, options = options,
#         );
#     end

#     # check result stability
#     solver = FEASolver(Direct, problem; xmin=xmin);
#     solver.vars = r.minimizer;
#     solver()
#     K, G = buckling(problem, solver.globalinfo, solver.elementinfo; u=solver.u);
#     @test isfinite(logdet(cholesky(K+G)))
#     sparse_eigvals, buckmodes = eigs(-G,K, nev=1, which=:LR)
#     smallest_pos_eigval = 1/sparse_eigvals[1]
#     @test smallest_pos_eigval >= 1.0

#     # using Makie
#     # using TopOpt.TrussTopOptProblems.TrussVisualization: visualize
#     # fig = visualize(problem; topology=r.minimizer)
# end

@testset "Tim buckling SDP constraint $problem_dim" for problem_dim in ["2d"] # , "3d"
    file_name = "tim_$(problem_dim).json"
    problem_file = joinpath(gm_ins_dir, file_name)

    mats = TrussFEAMaterial(10.0, 0.3)
    crossecs = TrussFEACrossSec(800.0)

    node_points, elements, _, _, fixities, load_cases = load_truss_json(problem_file)
    ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
    loads = load_cases["0"]

    problem = TrussProblem(
        Val{:Linear}, node_points, elements, loads, fixities, mats, crossecs
    )

    xmin = 0.0001 # minimum density
    p = 4.0 # penalty
    V = 0.5 # maximum volume fraction
    x0 = fill(V, ncells) # initial design

    solver = FEASolver(Direct, problem)
    ch = problem.ch
    dh = problem.ch.dh

    comp = TopOpt.Compliance(problem, solver)
    dp = TopOpt.Displacement(solver)
    assemble_k = TopOpt.AssembleK(problem)
    element_k = ElementK(solver)
    truss_element_kσ = TrussElementKσ(problem, solver)

    # * comliance minimization objective
    obj = comp
    c = 1.0 # buckling load multiplier

    function buckling_matrix_constr(x)
        # * Array(K + c*Kσ) ⋟ 0, PSD
        # * solve for the displacement
        u = dp(x)

        # * x -> Kes, construct all the element stiffness matrices
        # a list of small matrices for each element (cell)
        Kes = element_k(x)

        # * Kes -> K (global linear stiffness matrix)
        K = assemble_k(Kes)
        K = apply_boundary_with_meandiag!(K, ch)

        # * u_e, x_e -> Ksigma_e
        Kσs = truss_element_kσ(u, x)

        # * Kσs -> Kσ
        Kσ = assemble_k(Kσs)
        Kσ = apply_boundary_with_zerodiag!(Kσ, ch)

        return Array(K + c * Kσ)
    end

    function vol_constr(x)
        # volume fraction constraint
        return sum(x) / length(x) - V
    end

    # * Before optimization, check initial design stability
    @test isfinite(logdet(cholesky(buckling_matrix_constr(x0))))
    @test vol_constr(x0) <= 0

    m = Model(obj)
    addvar!(m, zeros(length(x0)), ones(length(x0)))
    Nonconvex.add_ineq_constraint!(m, vol_constr)
    Nonconvex.add_sd_constraint!(m, buckling_matrix_constr)

    Nonconvex.NonconvexCore.show_residuals[] = false
    alg = SDPBarrierAlg(; sub_alg=IpoptAlg())
    options = SDPBarrierOptions(; sub_options=IpoptOptions(; max_iter=200))
    r = Nonconvex.optimize(m, alg, x0; options=options)
    # println("$(r.convstate)")

    # * check result stability
    @test isfinite(logdet(cholesky(buckling_matrix_constr(r.minimizer))))

    # using Makie
    # import GLMakie
    # using TopOpt.TrussTopOptProblems.TrussVisualization: visualize
    # fig = visualize(problem; topology=r.minimizer)
    # Makie.display(fig)
end


text:
"""
    OrthogonalPolynomialRatio(P,x)

is a is equivalent to the vector `P[x,:] ./ P[x,2:end]`` but built from the
recurrence coefficients of `P`.
"""
mutable struct OrthogonalPolynomialRatio{T, PP<:AbstractQuasiMatrix{T}} <: AbstractCachedVector{T}
    P::PP # OPs
    x::T
    data::Vector{T}
    datasize::Tuple{Int}

    function OrthogonalPolynomialRatio{T, PP}(P::PP, x::T) where {T,PP<:AbstractQuasiMatrix{T}}
        μ = inv(sqrt(sum(orthogonalityweight(P))))
        new{T, PP}(P, x, [Base.unsafe_getindex(P,x,1)/Base.unsafe_getindex(P,x,2)], (1,))
    end
end

OrthogonalPolynomialRatio(P::AbstractQuasiMatrix{T}, x) where T = OrthogonalPolynomialRatio{T,typeof(P)}(P, convert(T, x))

size(K::OrthogonalPolynomialRatio) = (ℵ₀,)

function cache_filldata!(R::OrthogonalPolynomialRatio, inds)
    A,B,C = recurrencecoefficients(R.P)
    x = R.x
    data = R.data
    @inbounds for n in inds
        data[n] = inv(A[n]*x + B[n] - C[n] * data[n-1])
    end

end


text:
# const CPX_INFBOUND = 1e20
# const CPX_STR_PARAM_MAX = 512

function get_param_type(env::Env, indx::Int)
  ptype = Vector{Cint}(1)
  stat = @cpx_ccall(getparamtype, Cint, (
                    Ptr{Void},
                    Cint,
                    Ptr{Cint}
                    ),
                    env.ptr, convert(Cint,indx), ptype)
  if stat != 0
    throw(CplexError(env, stat))
  end
  if ptype[1] == 0
    ret = :None
  elseif ptype[1] == 1
    ret = :Int
  elseif ptype[1] == 2
    ret = :Double
  elseif ptype[1] == 3
    ret = :String
  elseif ptype[1] == 4
    ret = :Long
  else
    error("Parameter type not recognized")
  end

  return ret
end

get_param_type(env::Env, name::String) = get_param_type(env, paramName2Indx[name])

function set_param!(env::Env, _pindx::Int, val, ptype::Symbol)
  pindx = convert(Cint, _pindx)
  if ptype == :Int
    stat = @cpx_ccall(setintparam, Cint, (Ptr{Void}, Cint, Cint), env.ptr, pindx, convert(Cint,val))
  elseif ptype == :Double
    stat = @cpx_ccall(setdblparam, Cint, (Ptr{Void}, Cint, Cdouble), env.ptr, pindx, float(val))
  elseif ptype == :String
    stat = @cpx_ccall(setstrparam, Cint, (Ptr{Void}, Cint, Cstring), env.ptr, pindx, String(val))
  elseif ptype == :Long
    stat = @cpx_ccall(setlongparam, Cint, (Ptr{Void}, Cint, Clonglong), env.ptr, pindx, convert(Clonglong, val))
  elseif ptype == :None
    warn("Trying to set a parameter of type None; doing nothing")
  else
    error("Unrecognized parameter type")
  end
  if stat != 0
    throw(CplexError(env, stat))
  end
end

set_param!(env::Env, pindx::Int, val) = set_param!(env, pindx, val, get_param_type(env, pindx))

set_param!(env::Env, pname::String, val) = set_param!(env, paramName2Indx[pname], val)

# set_params!(env::Env, args...)
#   for (name, v) in args
#     set_param!(prob, string(name), v)
#   end
# end

function get_param(env::Env, pindx::Int, ptype::Symbol)
  if ptype == :Int
    val_int = Vector{Cint}(1)
    stat = @cpx_ccall(getintparam, Cint, (Ptr{Void}, Cint, Ptr{Cint}), env.ptr, convert(Cint,pindx), val_int)
    if stat != 0
      throw(CplexError(env, stat))
    end
    return val_int[1]
  elseif ptype == :Double
    val_double = Vector{Cdouble}(1)
    stat = @cpx_ccall(getdblparam, Cint, (Ptr{Void}, Cint, Ptr{Cdouble}), env.ptr, convert(Cint,pindx), val_double)
    if stat != 0
      throw(CplexError(env, stat))
    end
    return val_double[1]
  elseif ptype == :String
    buf = Vector{Cchar}(CPX_STR_PARAM_MAX) # max str param length is 512 in Cplex 12.51
    stat = @cpx_ccall(getstrparam, Cint, (Ptr{Void}, Cint, Ptr{Cchar}), env.ptr, convert(Cint,pindx), buf)
    if stat != 0
      throw(CplexError(env, stat))
    end
    return bytestring(pointer(buf))
  elseif ptype == :Long
    val_long = Vector{Clonglong}(1)
    stat = @cpx_ccall(getlongparam, Cint, (Ptr{Void}, Cint, Ptr{Clonglong}), env.ptr, convert(Cint,pindx), val_long)
    if stat != 0
      throw(CplexError(env, stat))
    end
    return val_long[1]
  elseif ptype == :None
    warn("Trying to set a parameter of type None; doing nothing")
  else
    error("Unrecognized parameter type")
  end
  nothing
end

get_param(env::Env, pindx::Int) = get_param(env, pindx, get_param_type(env, pindx))

get_param(env::Env, pname::String) = get_param(env, paramName2Indx[pname])

tune_param(model::Model) = tune_param(model, Dict(), Dict(), Dict())

function tune_param(model::Model, intfixed::Dict, dblfixed::Dict, strfixed::Dict)
  intkeys = Cint[k for k in keys(intfixed)]
  dblkeys = Cint[k for k in keys(dblfixed)]
  strkeys = Cint[k for k in keys(strfixed)]
  tune_stat = Vector{Cint}(1)
  stat = @cpx_ccall(tuneparam, Cint, (Ptr{Void},
                         Ptr{Void},
                         Cint,
                         Ptr{Cint},
                         Ptr{Cint},
                         Cint,
                         Ptr{Cint},
                         Ptr{Cdouble},
                         Cint,
                         Ptr{Cint},
                         Ptr{Ptr{Cchar}},
                         Ptr{Cint}),
                        model.env,
                        model.lp,
                        convert(Cint, length(intkeys)),
                        intkeys,
                        Cint[intfixed[int(k)] for k in intkeys],
                        convert(Cint, length(dblkeys)),
                        dblkeys,
                        Cdouble[dblfixed[int(k)] for k in dblkeys],
                        convert(Cint, length(strkeys)),
                        strkeys,
                        [strkeys[int(k)] for k in strkeys],
                        tune_stat)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  for param in keys(paramName2Indx)
    print(param * ": ")
    val = get_param(model.env, param)
    println(val)
  end
  return tune_stat[1]
end


text:
import_arviz() = _import_dependency("arviz", "arviz"; channel="conda-forge")

arviz_version() = VersionNumber(arviz.__version__)

function check_needs_update(; update=true)
    if arviz_version() < _min_arviz_version
        @warn "ArviZ.jl only officially supports arviz version $(_min_arviz_version) or " *
              "greater but found version $(arviz_version())."
        if update
            if update_arviz()
                # yay, but we still already imported the old version
                msg = """
                Please rebuild ArviZ.jl with `using Pkg; Pkg.build("ArviZ")` and re-launch Julia
                to continue.
                """
            else
                msg = """
                Could not automatically update arviz. Please manually update arviz, rebuild
                ArviZ.jl with `using Pkg; Pkg.build("ArviZ")`, and then re-launch Julia to
                continue.
                """
            end
            @warn msg
        end
    end
    return nothing
end

function check_needs_rebuild()
    if arviz_version() != _precompile_arviz_version
        msg = """
        ArviZ.jl was built using arviz version $(_precompile_arviz_version) but loaded with
        version $(arviz_version()). Please recompile with `using Pkg; Pkg.build("ArviZ")`
        and re-launch Julia to continue.
        """
        @warn msg
    end
    return nothing
end

function initialize_arviz()
    ispynull(arviz) || return nothing
    copy!(arviz, import_arviz())
    check_needs_update(; update=true)
    check_needs_rebuild()

    append!(SUPPORTED_GROUPS, map(Symbol, arviz.data.inference_data.SUPPORTED_GROUPS))

    pytype_mapping(arviz.InferenceData, InferenceData)

    # pytypemap-ing RcParams produces a Dict
    copy!(_rcParams, py"$(arviz).rcparams.rcParams"o)

    # use 1-based indexing by default within arviz
    rcParams["data.index_origin"] = 1

    # handle Bokeh showing ourselves
    rcParams["plot.bokeh.show"] = false

    initialize_xarray()
    initialize_numpy()
    return nothing
end

function initialize_xarray()
    ispynull(xarray) || return nothing
    copy!(xarray, _import_dependency("xarray", "xarray"; channel="conda-forge"))
    _import_dependency("dask", "dask"; channel="conda-forge")
    pytype_mapping(xarray.Dataset, Dataset)
    return nothing
end

function initialize_numpy()
    # Trigger NumPy initialization, see https://github.com/JuliaPy/PyCall.jl/issues/744
    PyObject([true])
    return nothing
end

function initialize_pandas()
    ispynull(pandas) || return nothing
    copy!(pandas, _import_dependency("pandas", "pandas"; channel="conda-forge"))
    return nothing
end

function update_arviz()
    # updating arviz can change other packages, so we always ask for permission
    if _using_conda() && _isyes(Base.prompt("Try updating arviz using conda? [Y/n]"))
        # this syntax isn't officially supported, but it works (for now)
        try
            Conda.add("arviz>=$_min_arviz_version"; channel="conda-forge")
            return true
        catch e
            println(e.msg)
        end
    end
    if _has_pip() && _isyes(Base.prompt("Try updating arviz using pip? [Y/n]"))
        # can't specify version lower bound, so update to latest
        try
            run(PyCall.python_cmd(`-m pip install --upgrade -- arviz`))
            return true
        catch e
            println(e.msg)
        end
    end
    return false
end

function _import_dependency(modulename, pkgname=modulename; channel=nothing)
    try
        return if channel === nothing
            pyimport_conda(modulename, pkgname)
        else
            pyimport_conda(modulename, pkgname, channel)
        end
    catch e
        if _has_pip() && _isyes(Base.prompt("Try installing $pkgname using pip? [Y/n]"))
            # installing with pip is riskier, so we ask for permission
            run(PyCall.python_cmd(`-m pip install -- $pkgname`))
            return pyimport(modulename)
        end
        # PyCall has a nice error message
        throw(e)
    end
end

_isyes(s) = isempty(s) || lowercase(strip(s)) ∈ ("y", "yes")
_isyes(::Nothing) = true

_using_conda() = PyCall.conda

_has_pip() = _has_pymodule("pip")

_has_pymodule(modulename) = !ispynull(pyimport_e(modulename))


text:
using Dates, TimeZones
using DelimitedFiles
using CSV, DataFrames
using Flux, Flux.Tracker
using MicroLogging
using ProgressMeter

using Mise

function mkpath_resdir(res_dir::String, output_size::Integer)
    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        Base.Filesystem.mkpath(res_dir * "$(i_pad)/")
    end
end

function post_station(input_path::String, model_path::String, res_dir::String,
    sample_size::Integer, output_size::Integer,
    test_fdate::ZonedDateTime, test_tdate::ZonedDateTime,
    eltype::DataType)

    stn_codes = [111121, 111123, 111131, 111141, 111142,
        111151, 111152, 111161, 111171, 111181,
        111191, 111201, 111212, 111221, 111231,
        111241, 111251, 111261, 111262, 111273,
        111274, 111281, 111291, 111301, 111311]
    stn_names = ["중구", "종로구", "용산구", "광진구", "성동구",
        "중랑구", "동대문구", "성북구", "도봉구", "은평구",
        "서대문구", "마포구", "강서구", "구로구", "영등포구",
        "동작구", "관악구", "강남구", "서초구", "송파구",
        "강동구", "금천구", "강북구", "양천구", "노원구"]

    stn_stats_df = DataFrame(
        code = Int64[],
        name = String[],
        lat = Real[],
        lon = Real[],
        colname = String[],
        RMSE = Real[],
        MAE = Real[],
        MSPE = Real[],
        MAPE = Real[])

    #stn_codes = [ 111123, 111142, 111161, 111273]
    #stn_names = ["종로구", "성동구", "성북구", "송파구"]

    @info "    Test with other station"
    p = Progress(length(stn_codes), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)

    stations = (; zip(Symbol.(stn_names), stn_codes)...)
    df = load_data_DNN("/input/jongro_seoul.csv", stations)
    df = filter_raw_data(df, test_fdate, test_tdate)

    for (stn_code, stn_name) in Base.Iterators.zip(stn_codes, stn_names)
        stn_df = filter_station_DNN(df, stn_code)

        row_PM10 = test_station(model_path, stn_df, :PM10, stn_code, stn_name,
            sample_size, output_size, res_dir,
            "PM10_" * string(stn_code) * "_$(stn_name)", test_fdate, test_tdate,
            eltype, false)
        row_PM25 = test_station(model_path, stn_df, :PM25, stn_code, stn_name,
            sample_size, output_size, res_dir,
            "PM25_" * string(stn_code) * "_$(stn_name)", test_fdate, test_tdate,
            eltype, false)

        push!(stn_stats_df, row_PM10)
        push!(stn_stats_df, row_PM25)
        ProgressMeter.next!(p);
    end

    CSV.write(res_dir * "station_stats.csv", stn_stats_df)
end

function post_feature(input_path::String, model_path::String, res_dir::String,
    sample_size::Integer, output_size::Integer,
    test_fdate::ZonedDateTime, test_tdate::ZonedDateTime,
    eltype::DataType)

    fea_stats_df = DataFrame(
        code = Int64[],
        name = String[],
        lat = Real[],
        lon = Real[],
        colname = String[],
        RMSE = Real[],
        MAE = Real[],
        MSPE = Real[],
        MAPE = Real[],
        rm_fea = String[])

    stn_codes = [111123]
    stn_names = ["jongro"]
    stn_code = 111123
    stn_name = "jongro"

    stations = (; zip(Symbol.(stn_names), stn_codes)...)
    stn_df = load_data_DNN("/input/jongro_seoul.csv", stations)
    stn_df = filter_raw_data(stn_df, test_fdate, test_tdate)
    rm_features = [[:SO2], [:CO], [:O3], [:NO2],[:PM10], [:PM25],  [:temp], [:u, :v], [:pres], [:humid], [:prep, :snow]]
    rm_features_str = ["SO2", "CO", "O3", "NO2", "PM10", "PM25", "temp", "u_v", "pres", "humid", "prep_snow"]

    #rm_features = [[:PM10], [:PM25]]
    #rm_features_str = ["PM10", "PM25"]

    @info "    Test by removing some features"
    p = Progress(length(rm_features), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)

    features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
    norm_prefix = "norm_"
    _norm_feas = [Symbol(eval(norm_prefix * String(f))) for f in features]
    stn_df2 = copy(stn_df)
    #zscore!(stn_df2, features, _norm_feas)
    min_max_scaling!(stn_df2, features, _norm_feas)

    for (idx, r_fea) in enumerate(rm_features)
        r_fea_str = join(string.(r_fea))

        rmf_df = copy(stn_df)
        for r_f in r_fea
            rmf_df[:, r_f] .= 0.0
        end

        features = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
        norm_prefix = "norm_"
        _norm_feas = [Symbol(eval(norm_prefix * String(f))) for f in features]
        #zscore!(rmf_df, features, _norm_feas)
        min_max_scaling!(rmf_df, features, _norm_feas)

        for r_f in r_fea
            norm_r_fea = Symbol(norm_prefix, string(r_f))
            # replace NaN to zero at r_f column (necessary when testing removing some features)
            rmf_df[findall(x -> isnan(x), rmf_df[:, norm_r_fea]), norm_r_fea] .= 0.0
        end

        row_PM10 = test_station(model_path, rmf_df, stn_df2, :PM10, stn_code, stn_name,
            sample_size, output_size, res_dir, "PM10_rm_$(r_fea_str)", test_fdate, test_tdate,
            eltype, true)
        row_PM25 = test_station(model_path, rmf_df, stn_df2, :PM25, stn_code, stn_name,
            sample_size, output_size, res_dir, "PM25_rm_$(r_fea_str)", test_fdate, test_tdate,
            eltype, true)

        push!(row_PM10, rm_features_str[idx])
        push!(row_PM25, rm_features_str[idx])

        push!(fea_stats_df, row_PM10)
        push!(fea_stats_df, row_PM25)

        ProgressMeter.next!(p);
    end

    row_PM10 = test_station(model_path, stn_df, :PM10, stn_code, stn_name,
            sample_size, output_size, res_dir, "PM10_rm_FULL", test_fdate, test_tdate, false)
    row_PM25 = test_station(model_path, stn_df, :PM25, stn_code, stn_name,
        sample_size, output_size, res_dir, "PM25_rm_FULL", test_fdate, test_tdate, false)

    push!(row_PM10, "FULL")
    push!(row_PM25, "FULL")
    push!(fea_stats_df, row_PM10)
    push!(fea_stats_df, row_PM25)
    CSV.write(res_dir * "feature_stats.csv", fea_stats_df)

end

function post_forecast(input_path::String, model_path::String, res_dir::String,
    sample_size::Integer, output_size::Integer,
    test_fdate::ZonedDateTime, test_tdate::ZonedDateTime,
    eltype::DataType)

    stn_codes = [111121, 111123, 111131, 111141, 111142,
        111151, 111152, 111161, 111171, 111181,
        111191, 111201, 111212, 111221, 111231,
        111241, 111251, 111261, 111262, 111273,
        111274, 111281, 111291, 111301, 111311]
    stn_names = ["중구", "종로구", "용산구", "광진구", "성동구",
        "중랑구", "동대문구", "성북구", "도봉구", "은평구",
        "서대문구", "마포구", "강서구", "구로구", "영등포구",
        "동작구", "관악구", "강남구", "서초구", "송파구",
        "강동구", "금천구", "강북구", "양천구", "노원구"]

    fore_stats_df = DataFrame(
        code = Int64[],
        name = String[],
        lat = Real[],
        lon = Real[],
        colname = String[],
        fore_all = Real[],
        fore_high = Real[])

    @info "    Test with other station"
    p = Progress(length(stn_codes), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)

    stations = (; zip(Symbol.(stn_names), stn_codes)...)
    df = load_data_DNN("/input/jongro_seoul.csv", stations)
    df = filter_raw_data(df, test_fdate, test_tdate)

    for (stn_code, stn_name) in Base.Iterators.zip(stn_codes, stn_names)
        stn_df = filter_station_DNN(df, stn_code)

        row_PM10 = test_classification(model_path, stn_df, :PM10, stn_code, stn_name,
            sample_size, output_size, res_dir,
            "PM10_" * string(stn_code) * "_$(stn_name)", test_fdate, test_tdate,
            eltype, false)
        row_PM25 = test_classification(model_path, stn_df, :PM25, stn_code, stn_name,
            sample_size, output_size, res_dir,
            "PM25_" * string(stn_code) * "_$(stn_name)", test_fdate, test_tdate,
            eltype, false)

        push!(fore_stats_df, row_PM10)
        push!(fore_stats_df, row_PM25)
        ProgressMeter.next!(p);
    end

    CSV.write(res_dir * "forecast_stats.csv", fore_stats_df)
end

function run()
    input_path = "/input/input.csv"
    model_path = "/mnt/"

    sample_size = 48
    output_size = 24

    # For GPU
    eltype::DataType = Float64

    test_fdate = ZonedDateTime(2018, 1, 1, 1, tz"Asia/Seoul")
    test_tdate = ZonedDateTime(2018, 12, 31, 23, tz"Asia/Seoul")

    @info "Postprocessing per station"
    res_dir = "/mnt/post/station/"
    Base.Filesystem.mkpath(res_dir)
    mkpath_resdir(res_dir, output_size)
    post_station(input_path, model_path, res_dir,
        sample_size, output_size, test_fdate, test_tdate, eltype)

    #=
    @info "Postprocessing per feature removing"
    res_dir = "/mnt/post/feature/"
    Base.Filesystem.mkpath(res_dir)
    mkpath_resdir(res_dir, output_size)
    post_feature(input_path, model_path, res_dir,
        sample_size, output_size, test_fdate, test_tdate, eltype)

    @info "Postprocessing for forecasting"
    res_dir = "/mnt/post/forecast/"
    Base.Filesystem.mkpath(res_dir)
    mkpath_resdir(res_dir, output_size)
    post_forecast(input_path, model_path, res_dir,
        sample_size, output_size, test_fdate, test_tdate, eltype)
    =#
end

run()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

const name = "CGAL"
const version = v"5.0"

# Collection of sources required to build CGAL
const sources = [
    "https://github.com/CGAL/cgal.git" => "6e9a233d552262d61d22677fa2ae9d5b69f413e8",
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/benlorenz/boostBuilder/releases/download/v1.71.0-1/build_boost.v1.71.0.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/GMP-v6.1.2-1/build_GMP.v6.1.2.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/MPFR-v4.0.2-1/build_MPFR.v4.0.2.jl",
]

# Bash recipe for building across all platforms
const script = raw"""
## pre-build setup
# exit on error
set -eu

# check c++ standard reported by the compiler
# CGAL uses CMake's try_run to check if it needs to link with Boost.Thread
# depending on the c++ standard supported by the compiler. From c++11 onwards,
# CGAL doesn't require Boost.Thread (by default since CGAL 5.0)
__need_boost_thread=0
__cplusplus=$($CXX -x c++ -dM -E - </dev/null | grep __cplusplus | grep -o '[0-9]*')

## configure build
mkdir -p "$WORKSPACE/srcdir/build" && cd "$WORKSPACE/srcdir/build"

CMAKE_FLAGS=(
  ## cmake specific
  -DCMAKE_TOOLCHAIN_FILE="/opt/$target/$target.toolchain"
  -DCMAKE_BUILD_TYPE=Release
  -DCMAKE_INSTALL_PREFIX="$prefix"
  -DCMAKE_FIND_ROOT_PATH="$prefix"
  ## cgal specific
  -DCGAL_HEADER_ONLY=OFF
  -DWITH_CGAL_Core=ON
  -DWITH_CGAL_ImageIO=OFF
  -DWITH_CGAL_Qt5=OFF
  # try_run doesn't like cross-compilation: this is required
  -DCGAL_test_cpp_version_RUN_RES=$__need_boost_thread
  -DCGAL_test_cpp_version_RUN_RES__TRYRUN_OUTPUT=$__cplusplus
)

cmake ${CMAKE_FLAGS[@]} ../cgal*/

## and away we go..
cmake --build . --config Release --target install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
const platforms = filter(expand_gcc_versions(supported_platforms())) do p
    p.compiler_abi.gcc_version > :gcc6
end

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libCGAL", :libCGAL),
    LibraryProduct(prefix, "libCGAL_Core", :libCGAL_Core),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: rgcv/CGALBuilder

max_issues_repo_name: rgcv/CGALBuilder

max_forks_repo_name: rgcv/CGALBuilder


text:
immutable Float
    s::Uint64
    e::Int32
    de::Int32
end

Float() = Float(0,0,0)
Float(x,y) = Float(x,y,int32(0))
Float(d::FloatingPoint) = Float(_significand(d), _exponent(d))

# Consts
const Float10MSBits = 0xFFC0000000000000 # used normalize(Float)
const FloatSignMask = 0x8000000000000000 # used in normalize(Float)
const FloatSignificandSize = int32(64)

function normalize(v::Float)
    f = v.s
    e = v.e
    while (f & Float10MSBits) == 0
        f <<= 10
        e -= 10
    end
    while (f & FloatSignMask) == 0
        f <<= 1
        e -= 1
    end
    return Float(f,e)
end
function normalize(v::Float64)
    s = _significand(v); e = _exponent(v)
    while (s & HiddenBit(Float64)) == 0
        s <<= uint64(1)
        e -= int32(1)
    end
    s <<= uint64(FloatSignificandSize - SignificandSize(Float64))
    e -=  int32( FloatSignificandSize - SignificandSize(Float64))
    return Float(s, e)
end

# Float64
DenormalExponent(::Type{Float64}) = int32(-ExponentBias(Float64) + 1)
ExponentMask(::Type{Float64}) = 0x7FF0000000000000
PhysicalSignificandSize(::Type{Float64}) = int32(52)
SignificandSize(::Type{Float64}) = int32(53)
ExponentBias(::Type{Float64}) = int32(0x3FF + PhysicalSignificandSize(Float64))
SignificandMask(::Type{Float64}) = 0x000FFFFFFFFFFFFF
HiddenBit(::Type{Float64}) = 0x0010000000000000
uint_t(d::Float64) = reinterpret(Uint64,d)
# Float32
DenormalExponent(::Type{Float32}) = int32(-ExponentBias(Float32) + 1)
ExponentMask(::Type{Float32}) = 0x7F800000
PhysicalSignificandSize(::Type{Float32}) = int32(23)
SignificandSize(::Type{Float32}) = int32(24)
ExponentBias(::Type{Float32}) = int32(0x7F + PhysicalSignificandSize(Float32))
SignificandMask(::Type{Float32}) = 0x007FFFFF
HiddenBit(::Type{Float32}) = 0x00800000
uint_t(d::Float32) = reinterpret(Uint32,d)
# Float16
# TODO

function _exponent{T<:FloatingPoint}(d::T)
  isdenormal(d) && return DenormalExponent(T)
  biased_e::Int32 = int32((uint_t(d) & ExponentMask(T)) >> PhysicalSignificandSize(T))
  return int32(biased_e - ExponentBias(T))
end
function _significand{T<:FloatingPoint}(d::T)
  s = uint_t(d) & SignificandMask(T)
  return !isdenormal(d) ? s + HiddenBit(T) : s
end
isdenormal{T<:FloatingPoint}(d::T) = (uint_t(d) & ExponentMask(T)) == 0

function normalizedbound(f::FloatingPoint)
    v = Float(_significand(f),_exponent(f))
    m_plus = normalize(Float((v.s << 1) + 1, v.e - 1))
    if lowerboundaryiscloser(f)
        m_minus = Float((v.s << 2) - 1, v.e - 2)
    else
        m_minus = Float((v.s << 1) - 1, v.e - 1)
    end
    return Float(m_minus.s << (m_minus.e - m_plus.e), m_plus.e), m_plus
end
function lowerboundaryiscloser{T<:FloatingPoint}(f::T)
    physical_significand_is_zero = (uint_t(f) & SignificandMask(T)) == 0
    return physical_significand_is_zero && (_exponent(f) != DenormalExponent(T))
end

(-)(a::Float,b::Float) = Float(a.s - b.s,a.e,a.de)

const FloatM32 = 0xFFFFFFFF

function (*)(this::Float,other::Float)
    a::Uint64 = this.s >> 32
    b::Uint64 = this.s & FloatM32
    c::Uint64 = other.s >> 32
    d::Uint64 = other.s & FloatM32
    ac::Uint64 = a * c
    bc::Uint64 = b * c
    ad::Uint64 = a * d
    bd::Uint64 = b * d
    tmp::Uint64 = (bd >> 32) + (ad & FloatM32) + (bc & FloatM32)
    # By adding 1U << 31 to tmp we round the final result.
    # Halfway cases will be round up.
    tmp += uint64(1) << 31
    result_f::Uint64 = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32)
    return Float(result_f,this.e + other.e + 64,this.de)
end

const CachedPowers = Float[
  Float(0xfa8fd5a0081c0288, -1220, -348),
  Float(0xbaaee17fa23ebf76, -1193, -340),
  Float(0x8b16fb203055ac76, -1166, -332),
  Float(0xcf42894a5dce35ea, -1140, -324),
  Float(0x9a6bb0aa55653b2d, -1113, -316),
  Float(0xe61acf033d1a45df, -1087, -308),
  Float(0xab70fe17c79ac6ca, -1060, -300),
  Float(0xff77b1fcbebcdc4f, -1034, -292),
  Float(0xbe5691ef416bd60c, -1007, -284),
  Float(0x8dd01fad907ffc3c, -980, -276),
  Float(0xd3515c2831559a83, -954, -268),
  Float(0x9d71ac8fada6c9b5, -927, -260),
  Float(0xea9c227723ee8bcb, -901, -252),
  Float(0xaecc49914078536d, -874, -244),
  Float(0x823c12795db6ce57, -847, -236),
  Float(0xc21094364dfb5637, -821, -228),
  Float(0x9096ea6f3848984f, -794, -220),
  Float(0xd77485cb25823ac7, -768, -212),
  Float(0xa086cfcd97bf97f4, -741, -204),
  Float(0xef340a98172aace5, -715, -196),
  Float(0xb23867fb2a35b28e, -688, -188),
  Float(0x84c8d4dfd2c63f3b, -661, -180),
  Float(0xc5dd44271ad3cdba, -635, -172),
  Float(0x936b9fcebb25c996, -608, -164),
  Float(0xdbac6c247d62a584, -582, -156),
  Float(0xa3ab66580d5fdaf6, -555, -148),
  Float(0xf3e2f893dec3f126, -529, -140),
  Float(0xb5b5ada8aaff80b8, -502, -132),
  Float(0x87625f056c7c4a8b, -475, -124),
  Float(0xc9bcff6034c13053, -449, -116),
  Float(0x964e858c91ba2655, -422, -108),
  Float(0xdff9772470297ebd, -396, -100),
  Float(0xa6dfbd9fb8e5b88f, -369, -92),
  Float(0xf8a95fcf88747d94, -343, -84),
  Float(0xb94470938fa89bcf, -316, -76),
  Float(0x8a08f0f8bf0f156b, -289, -68),
  Float(0xcdb02555653131b6, -263, -60),
  Float(0x993fe2c6d07b7fac, -236, -52),
  Float(0xe45c10c42a2b3b06, -210, -44),
  Float(0xaa242499697392d3, -183, -36),
  Float(0xfd87b5f28300ca0e, -157, -28),
  Float(0xbce5086492111aeb, -130, -20),
  Float(0x8cbccc096f5088cc, -103, -12),
  Float(0xd1b71758e219652c, -77, -4),
  Float(0x9c40000000000000, -50, 4),
  Float(0xe8d4a51000000000, -24, 12),
  Float(0xad78ebc5ac620000, 3, 20),
  Float(0x813f3978f8940984, 30, 28),
  Float(0xc097ce7bc90715b3, 56, 36),
  Float(0x8f7e32ce7bea5c70, 83, 44),
  Float(0xd5d238a4abe98068, 109, 52),
  Float(0x9f4f2726179a2245, 136, 60),
  Float(0xed63a231d4c4fb27, 162, 68),
  Float(0xb0de65388cc8ada8, 189, 76),
  Float(0x83c7088e1aab65db, 216, 84),
  Float(0xc45d1df942711d9a, 242, 92),
  Float(0x924d692ca61be758, 269, 100),
  Float(0xda01ee641a708dea, 295, 108),
  Float(0xa26da3999aef774a, 322, 116),
  Float(0xf209787bb47d6b85, 348, 124),
  Float(0xb454e4a179dd1877, 375, 132),
  Float(0x865b86925b9bc5c2, 402, 140),
  Float(0xc83553c5c8965d3d, 428, 148),
  Float(0x952ab45cfa97a0b3, 455, 156),
  Float(0xde469fbd99a05fe3, 481, 164),
  Float(0xa59bc234db398c25, 508, 172),
  Float(0xf6c69a72a3989f5c, 534, 180),
  Float(0xb7dcbf5354e9bece, 561, 188),
  Float(0x88fcf317f22241e2, 588, 196),
  Float(0xcc20ce9bd35c78a5, 614, 204),
  Float(0x98165af37b2153df, 641, 212),
  Float(0xe2a0b5dc971f303a, 667, 220),
  Float(0xa8d9d1535ce3b396, 694, 228),
  Float(0xfb9b7cd9a4a7443c, 720, 236),
  Float(0xbb764c4ca7a44410, 747, 244),
  Float(0x8bab8eefb6409c1a, 774, 252),
  Float(0xd01fef10a657842c, 800, 260),
  Float(0x9b10a4e5e9913129, 827, 268),
  Float(0xe7109bfba19c0c9d, 853, 276),
  Float(0xac2820d9623bf429, 880, 284),
  Float(0x80444b5e7aa7cf85, 907, 292),
  Float(0xbf21e44003acdd2d, 933, 300),
  Float(0x8e679c2f5e44ff8f, 960, 308),
  Float(0xd433179d9c8cb841, 986, 316),
  Float(0x9e19db92b4e31ba9, 1013, 324),
  Float(0xeb96bf6ebadf77d9, 1039, 332),
  Float(0xaf87023b9bf0ee6b, 1066, 340)]

const CachedPowersLength = length(CachedPowers)
const CachedPowersOffset = 348  # -1 * the first decimal_exponent.
const D_1_LOG2_10 = 0.30102999566398114  #  1 / lg(10)
# Difference between the decimal exponents in the table above.
const DecimalExponentDistance = 8
const MinDecimalExponent = -348
const MaxDecimalExponent = 340

function GetCachedPowerForBinaryExponentRange(
            min_exponent::Int,
            max_exponent::Int)
    kQ = 64
    k = ceil((min_exponent+kQ-1)*D_1_LOG2_10)
    foo = CachedPowersOffset
    index = div(foo+itrunc(k)-1,DecimalExponentDistance) + 1
    cp = CachedPowers[index+1]
    return cp
end

function GetCachedPowerForDecimalExponent(
    requested_exponent::Int)
    index = div(requested_exponent + CachedPowersOffset,DecimalExponentDistance)
    cp = CachedPowers[index+1]
    return cp
end