text:
module TriMeshes

using StaticArrays
using LinearAlgebra
using Polynomials

const Vertex = SVector{3,Float64}
const Point2D = SVector{2,Float64}
const QuadF = SVector{4,Int64}
const QuadI = SVector{4,Int64}
const TriI = SVector{3,Int64}

include("Splines.jl")
include("Spline_exports.jl")

export scale, Edge, Face, Net, vertex!, face!, stl, wrap
export quadRSlice, quadVSlice, outerSlice, innerSlice
export cylinder!, flatDisk!, tube!, flatRing!
export transformVerts!, apply!
export areaXY, magnitude, normalize
export cube!

struct Edge
	from::Integer
	to::Integer
end

struct Face
	AB::Edge
	BC::Edge
	CA::Edge
end

struct Net
	vertices::Vector{Vertex}
	vhashes::Dict{UInt64, Int64}
	faces::Vector{Face}
	edges::Vector{Edge}
	Net() = new(Vector{Vertex}(), Dict{UInt64, Int64}(), Vector{Face}(), Vector{Edge}())
end

vertex!(n::Net, x, y, z) = vertex!(n, Vertex(x, y, z))
vertex!(n::Net, x, y) = vertex!(n, Vertex(x, y, 0.0))
vertex!(v::Vector{Vertex}, x, y, z) = vertex!(v, Vertex(x, y, z))
function vertex!(n::Net, v::Vertex)
	h = hash(v)
	i = get(n.vhashes, h, 0)
	if i == 0
		i = vertex!(n.vertices, v)
		n.vhashes[h] = i
	end
	i
end

function vertex!(vs::Vector{Vertex}, v::Vertex)
	push!(vs, v)
	length(vs)
end

scale(n::Net, x, y, z) = map!(v->Vertex(v.x*x, v.y*y, v.z*z), n.vertices, n.vertices)
scale(n::Net, s::Vertex)  = map!(v->Vertex(v.x*s.x, v.y*s.y, v.z*s.z), n.vertices, n.vertices)

face!(n::Net, a::Int, b::Int, c::Int) = face!(n, Face(Edge(a, b), Edge(b, c), Edge(c, a)))
function face!(n::Net, f::Face)
	push!(n.faces, f)
	length(n.faces)
end

apply!(net, vs, fn) = foreach(v->net.vertices[v] = fn(net.vertices[v]), vs)

function transformVerts!(n::Net, f, vs::Integer, ve::Integer)
	n.vertices[vs:ve] = map(f, n.vertices[vs:ve])
end

function wrap(n::Net, az2r)

	function trans(v::Vertex)
		a = 2pi * v.x / 360
		r = az2r(a, v.z) + v.y
		Vertex(r * cos(a), r * sin(a), v.z)
	end

	map!(trans, n.vertices, n.vertices)
end

quadRSlice(az2r, a, astep, zmin, zmax) = QuadF[
		az2r(a, zmin),
		az2r(a+astep, zmin),
		az2r(a, zmax),
		az2r(a+astep, zmax)
	]

quadVSlice(net::Net, a::Float64, astep, zmin, zmax, rs) = QuadI[
		vertex!(net, rs[1] * cos(a), rs[1] * sin(a), zmin),
		vertex!(net, rs[2] * cos(a+astep), rs[2] * sin(a+astep), zmin),
		vertex!(net, rs[3] * cos(a), rs[3] * sin(a), zmax),
		vertex!(net, rs[4] * cos(a+astep), rs[4] * sin(a+astep), zmax)
	]

function outerSlice(net::Net, az2r::Function, asteps, zmin, zmax)
	vs = zeros(Int, 4asteps)
	astep = 2pi / asteps
	vi = 1
	for a in astep:astep:2pi
		vs[vi:vi+3] = quadVSlice(net, a, astep, zmin, zmax, quadRSlice(az2r, a, astep, zmin, zmax))
		face!(net, vs[vi], vs[vi+1], vs[vi+2])
		face!(net, vs[vi+3], vs[vi+2], vs[vi+1])
		vi += 4
	end
	vs
end

quadVSlice(net::Net, o::Vertex, a::Float64, astep, zstep, rs) = QuadI[
		vertex!(net, o.x+rs[1] * cos(a), o.y+rs[1] * sin(a), o.z),
		vertex!(net, o.x+rs[2] * cos(a+astep), o.y+rs[2] * sin(a+astep), o.z),
		vertex!(net, o.x+rs[3] * cos(a), o.y+rs[3] * sin(a), o.z+zstep),
		vertex!(net, o.x+rs[4] * cos(a+astep), o.y+rs[4] * sin(a+astep), o.z+zstep)
	]

function outerSlice(net::Net, o::Vertex, az2r::Function, asteps, zstep)
	vs = zeros(Int, 4asteps)
	astep = 2pi / asteps
	vi = 1
	for a in astep:astep:2pi
		vs[vi:vi+3] = quadVSlice(net, o, a, astep, zstep, quadRSlice(az2r, a, astep, o.z, o.z+zstep))
		face!(net, vs[vi], vs[vi+1], vs[vi+2])
		face!(net, vs[vi+3], vs[vi+2], vs[vi+1])
		vi += 4
	end
	vs
end

quadVSlice(net, o, t, a, astep, rs) = QuadI[
		vertex!(net, o.x+rs[1] * cos(a), o.y+rs[1] * sin(a), o.z),
		vertex!(net, o.x+rs[2] * cos(a+astep), o.y+rs[2] * sin(a+astep), o.z),
		vertex!(net, t.x+rs[3] * cos(a), t.y+rs[3] * sin(a), t.z),
		vertex!(net, t.x+rs[4] * cos(a+astep), t.y+rs[4] * sin(a+astep), t.z)
	]

function outerSlice(net, o, t, az2r, astart, asteps)
	vs = zeros(Int, 4asteps)
	astep = 2pi / asteps
	vi = 1
	for a in astep+astart:astep:2pi+astart
		vs[vi:vi+3] = quadVSlice(net, o, t, a, astep, quadRSlice(az2r, a, astep, o.z, t.z))
		face!(net, vs[vi], vs[vi+1], vs[vi+2])
		face!(net, vs[vi+3], vs[vi+2], vs[vi+1])
		vi += 4
	end
	vs
end

function innerSlice(net, az2r, astep, zmin, zmax)
	vns = zeros(Int, 4asteps)
	astep = 2pi / asteps
	vi = 1
	for a in astep:astep:2pi
		vns[vi:vi+3] = quadVSlice(net, a, astep, zmin, zmax, quadRSlice(az2r, a, astep, zmin, zmax))
		face!(net, vns[vi+2], vns[vi+1], vns[vi])
		face!(net, vns[vi+1], vns[vi+2], vns[vi+3])
		vi += 4
	end
	vs
end

function cylinder!(net::Net, centreL, centreH, asteps, az2r)
	flat(net, centreL, az2r, v.asteps, :down)
	outerSlice(net, az2r, v.asteps, centreL.z, centreH.z)
	flat(net, centreH, az2r, v.asteps, :up)
end

function flatDisk!(net::Net, centre, az2r, astart, asteps, upordown)
	astep = 2pi/asteps
	for a in astep+astart:astep:2pi+astart
		vn1 = vertex!(net, centre.x + az2r(a, centre.z) * cos(a), centre.y + az2r(a, centre.z) * sin(a), centre.z)
		vn2 = vertex!(net, centre.x, centre.y, centre.z)
		vn3 = vertex!(net, centre.x + az2r(a+astep, centre.z) * cos(a+astep), centre.y + az2r(a+astep, centre.z) * sin(a+astep), centre.z)
		if upordown ==
			face!(net, vn3, vn2, vn1)
		else
			face!(net, vn1, vn2, vn3)
		end
	end
end

function flatRing!(n::Net, originL, astart, asteps, router, rinner, thickness)
	originH = originL + Vertex(0,0,thickness)
	astep = 2pi / asteps
	for a in astep+astart:astep:2pi+astart
		vno = map(p->vertex!(n, p[1] + Vertex(router*cos(p[2]), router*sin(p[2]), 0)), [(originL, a), (originL, a+astep), (originH, a), (originH, a+astep)])
		vni = map(p->vertex!(n, p[1] + Vertex(rinner*cos(p[2]), rinner*sin(p[2]), 0)), [(originL, a), (originL, a+astep), (originH, a), (originH, a+astep)])

		face!(n, vno[1], vno[2], vno[3])
		face!(n, vno[4], vno[3], vno[2])
		face!(n, vni[3], vni[2], vni[1])
		face!(n, vni[2], vni[3], vni[4])

		face!(n, vno[1], vni[1], vni[2])
		face!(n, vni[2], vno[2], vno[1])
		face!(n, vno[4], vni[3], vno[3])
		face!(n, vno[4], vni[4], vni[3])
	end

end

function tube!(net, radius, sides, astart, srcnet, srcvertices, startz, endz)
	for vtxn in srcvertices
		if srcnet.vertices[vtxn[1]].z == startz
			flatDisk!(net, srcnet.vertices[vtxn[1]], (a,z)->radius, astart, sides, :down)
		end
		outerSlice(net, srcnet.vertices[vtxn[1]], srcnet.vertices[vtxn[2]], (a,z)->radius, astart, sides)
		if srcnet.vertices[vtxn[2]].z == endz
			flatDisk!(net, srcnet.vertices[vtxn[2]], (a,z)->radius, astart, sides, :up)
		end
	end

end

export Vertex, angleXY, angleYX, angleYZ, angleZY, angleXZ, angleZX, rotate, translate, transformVerts!, Edge, Face, Net, vertex!, face!, abc, abci, areaXY, magnitude

import Base.show
function show(io::IO,v::Vertex)
       print(io,"Vertex($(round(v.x, digits=4)), $(round(v.y, digits=4)), $(round(v.z, digits=4)))")
end

function LinearAlgebra.dot(v1::Vertex, v2::Vertex)
	vecdot([v1.x, v1.y, v1.z], [v2.x, v2.y, v2.z])
end

function LinearAlgebra.cross(v1::Vertex, v2::Vertex)
	cross([v1.x, v1.y, v1.z], [v2.x, v2.y, v2.z])
end

function LinearAlgebra.normalize(v::Vertex)
      n = normalize([v.x, v.y, v.z])
      Vertex(n[1], n[2], n[3])
end

function angle(x,y)
	# atan(y,x) - take point (x,y) - angle from x axis to point
	if x == y == 0
		0.0
	else
		atan(y,x)
	end
end

function angleXY(v::Vertex)
	# project onto plane XY, angle from X axis to point
	angle(v.x, v.y)
end

function angleYX(v::Vertex)
	# project onto plane XY, angle from Y axis to point
	angle(v.y, v.x)
end

function angleYZ(v::Vertex)
	# project onto plane YZ, angle from Y axis to point
	angle(v.y, v.z)
end

function angleZY(v::Vertex)
	# project onto plane YZ, angle from Z axis to point
	angle(v.z, v.y)
end

function angleXZ(v::Vertex)
	# project onto plane XZ, angle from X axis to point
	angle(v.x, v.z)
end

function angleZX(v::Vertex)
	# project onto plane XZ, angle from Z axis to point
	angle(v.z, v.x)
end

function rotate(v::Vertex, a::Vertex)
	if abs(a.x) > 0
		v = Vertex(v.x, v.y * cos(a.x) - v.z * sin(a.x), v.y * sin(a.x) + v.z * cos(a.x))
	end
	if abs(a.y) > 0
		v = Vertex(v.z * sin(a.y) + v.x * cos(a.y), v.y, v.z * cos(a.y) - v.x * sin(a.y))
	end
	if abs(a.z) > 0
		v = Vertex(v.x * cos(a.z) - v.y * sin(a.z), v.x * sin(a.z) + v.y * cos(a.z), v.z)
	end
	v
end

rotate(x, y, z, xa, ya, za) = rotate(Vertex(x,y,z), Vertex(xa, ya, za))
rotate(v::Vertex, xa, ya, za) = rotate(v, Vertex(xa, ya, za))
rotate(x, y, z, v::Vertex) = rotate(Vertex(x,y,z), v)

translate(v::Vertex, t::Vertex) = v+t

abci(f::Face) = [f.ab.from, f.ab.to, f.ca.from]
abci(n::Net, i::Integer) = abci(n.faces[i])

abc(n::Net, f::Face) = n.vertices[abci(f)]
abc(n::Net, i::Integer) = abc(n, n.faces[i])

magnitude(v::Vertex) = sqrt(v.x^2 + v.y^2 + v.z^2)

function areaXY(a::Vertex, b::Vertex, c::Vertex)

	ab = b-a
	ca = a-c

	h = abs(magnitude(ab) * sin(angleXY(ab) - angleXY(ca)))
	0.5 * h * magnitude(ca)
end
areaXY(n::Net, f::Face) = areaXY(abc(n, f)...)

function quad!(n::Net, v1, v2, v3, v4)
	a = vertex!(n, v1)
	b = vertex!(n, v2)
	c = vertex!(n, v3)
	d = vertex!(n, v4)

	face!(n, a, b, c)
	face!(n, c, d, a)
end

function cube!(n, w, h, d; origin=Vertex(0,0,0))

	v(x,y,z) = Vertex(origin.x+x, origin.y+y, origin.z+z)

	# front
	quad!(n, v(0,0,0), v(0,h,0), v(w,h,0), v(w,0,0))
	# left
	quad!(n, v(0,0,d), v(0,h,d), v(0,h,0), v(0,0,0))
	# back
	quad!(n, v(w,0,d), v(w,h,d), v(0,h,d), v(0,0,d))
	# right
	quad!(n, v(w,0,0), v(w,h,0), v(w,h,d), v(w,0,d))
	# top
	quad!(n, v(0,h,0), v(0,h,d), v(w,h,d), v(w,h,0))
	# bottom
	quad!(n, v(0,0,d), v(0,0,0), v(w,0,d), v(w,0,0))
	nothing;
end

###
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Private functions and macros.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Dictionary that translates table type to environment.
const _latex_table_env = Dict(
    :array     => "array",
    :longtable => "longtable",
    :tabular   => "tabular"
)

# Apply an alignment to a LaTeX table cell.
function _latex_apply_cell_alignment(
    str::String,
    alignment::Symbol,
    j::Int,
    num_printed_cols::Int,
    show_row_number::Bool,
    vlines::Vector{Int},
    left_vline::String,
    mid_vline::String,
    right_vline::String
)
    a = _latex_alignment(alignment)
    aux_j = show_row_number ? j + 1 : j

    # We only need to add left vertical line if it is the first
    # column.
    lvline = if (0 ∈ vlines) && (aux_j - 1 == 0)
        left_vline
    else
        ""
    end

    # For the right vertical line, we must check if it is a mid line
    # or right line.
    if aux_j ∈ vlines
        rvline = (j == num_printed_cols) ? right_vline : mid_vline
    else
        rvline = ""
    end

    # Wrap the data into the multicolumn environment.
    return _latex_envs(str, "multicolumn{1}{$(lvline)$(a)$(rvline)}")
end

# Convert the alignment symbol into a LaTeX alignment string.
function _latex_alignment(s::Symbol)
    if (s == :l) || (s == :L)
        return "l"
    elseif (s == :c) || (s == :C)
        return "c"
    elseif (s == :r) || (s == :R)
        return "r"
    else
        error("Invalid LaTeX alignment symbol: $s.")
    end
end

# Get the LaTeX table description (alignment and vertical columns).
function _latex_table_desc(
    id_cols::AbstractVector,
    alignment::Vector{Symbol},
    show_row_name::Bool,
    row_name_alignment::Symbol,
    show_row_number::Bool,
    row_number_alignment::Symbol,
    vlines::AbstractVector,
    left_vline::AbstractString,
    mid_vline::AbstractString,
    right_vline::AbstractString
)
    Δc = show_row_number ? 1 : 0

    str = "{"

    0 ∈ vlines && (str *= left_vline)

    # Add the alignment information of the row number column if required.
    Δc = 0
    if show_row_number
        Δc += 1
        str *= _latex_alignment(row_number_alignment)
        Δc ∈ vlines && (str *= mid_vline)
    end

    if show_row_name
        Δc += 1
        str *= _latex_alignment(row_name_alignment)
        Δc ∈ vlines && (str *= mid_vline)
    end

    # Process the alignment of all the other columns.
    for i = 1:length(id_cols)
        str *= _latex_alignment(alignment[id_cols[i]])

        if i+Δc ∈ vlines
            if i != length(id_cols)
                str *= mid_vline
            else
                str *= right_vline
            end
        end
    end

    str *= "}"

    return str
end

# Wrap the `text` into LaTeX environment(s).
function _latex_envs(text::AbstractString, envs::Vector{String})
    return _latex_envs(text, envs, length(envs))
end

function _latex_envs(text::AbstractString, env::String)
    if !isempty(text)
        return "\\" * string(env) * "{" * text * "}"
    else
        return ""
    end
end

function _latex_envs(text::AbstractString, envs::Vector{String}, i::Int)
    @inbounds if i > 0
        str = _latex_envs(text, envs[i])
        return _latex_envs(str, envs, i -= 1)
    end

    return text
end

ext: jl


text:
using JuMP; using NCDatasets; using AxisArrays; using CalliopeJuMP.Util;

function load_network_constraints(model_dict)
    sets = model_dict["sets"]
    constraint_dict = Dict()

    if haskey(sets, "loc_techs_symmetric_transmission_constraint")
        constraint_dict["symmetric_transmission_constraint"] = (
            build_constraint(
                ["loc_techs_symmetric_transmission_constraint"],
                "symmetric_transmission", model_dict
            )
        )
    end

    return constraint_dict
end

function symmetric_transmission_constraint_rule(backend_model, set_indices, model_dict)
    """
    Constrain e_cap symmetrically for transmission nodes. Transmission techs only.
    """
    loc_tech = set_indices[1]
    variables = model_dict["variables"]

    loc_tech_remote = model_dict["parameters"]["lookup_remotes"][loc_tech]

    return @constraint(backend_model,
        variables["energy_cap"][loc_tech] == variables["energy_cap"][loc_tech_remote]
    )
end

max_stars_repo_name: calliope-project/CalliopeJuMP

max_issues_repo_name: calliope-project/CalliopeJuMP

max_forks_repo_name: calliope-project/CalliopeJuMP


text:
function applyperipheryics!(system::CVSystem,old=Dict("a"=>0),restart="no")
    numupper = 0;
    if restart == "yes"
        branches = old["branches"];
    end
    for i = 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            if restart == "no"
                system.branches.term[i].P[1,1] = (
                    system.branches.P[i][1,system.solverparams.JL]*mmHgToPa);
                system.branches.term[i].V[1,1] = (
                system.branches.term[i].P[1,1]*system.branches.term[i].C[1] +
                system.branches.term[i].V0[1]);
                system.branches.term[i].P[1,2] = 80*mmHgToPa;
                system.branches.term[i].V[1,2] = (
                system.branches.term[i].P[1,2]*system.branches.term[i].C[2] +
                system.branches.term[i].V0[2]);
                system.branches.term[i].P[1,3] = 40*mmHgToPa;
                system.branches.term[i].V[1,3] = (
                system.branches.term[i].P[1,3]*system.branches.term[i].C[3] +
                system.branches.term[i].V0[3]);
                system.branches.term[i].P[1,4] = 5*mmHgToPa;
                # system.branches.term[i].P[1,4] = 4.9*mmHgToPa;
                system.branches.term[i].V[1,4] = (
                system.branches.term[i].P[1,4]*system.branches.term[i].C[4] +
                system.branches.term[i].V0[4]);
                # system.branches.term[i].P[1,5] = 4.5*mmHgToPa;
                system.branches.term[i].P[1,5] = 4.5*mmHgToPa;
                system.branches.term[i].V[1,5] = (
                system.branches.term[i].P[1,5]*system.branches.term[i].C[5] +
                system.branches.term[i].V0[5]);
            elseif restart == "yes"
                term = branches["term"][i];
                system.branches.term[i].P[1,1] = term["P"][end,1];
                system.branches.term[i].V[1,1] = term["V"][end,1];
                system.branches.term[i].Q[1,1] = term["Q"][end,1];
                system.branches.term[i].P[1,2] = term["P"][end,2];
                system.branches.term[i].V[1,2] = term["V"][end,2];
                system.branches.term[i].Q[1,2] = term["Q"][end,2];
                system.branches.term[i].P[1,3] = term["P"][end,3];
                system.branches.term[i].V[1,3] = term["V"][end,3];
                system.branches.term[i].Q[1,3] = term["Q"][end,3];
                system.branches.term[i].P[1,4] = term["P"][end,4];
                system.branches.term[i].V[1,4] = term["V"][end,4];
                system.branches.term[i].Q[1,4] = term["Q"][end,4];
                system.branches.term[i].P[1,5] = term["P"][end,5];
                system.branches.term[i].V[1,5] = term["V"][end,5];
                system.branches.term[i].Q[1,5] = term["Q"][end,5];
            end
            if system.branches.group[i] == "upper"
                numupper+=1
            end
        end
    end

    if numupper > 0
        if restart == "no"
            system.svc.P[1,1] = 4.5*mmHgToPa;
            system.svc.V[1,1] = system.svc.P[1,1]*system.svc.C + system.svc.V0;
            system.ivc.P[1,1] = 4.5*mmHgToPa;
            system.ivc.V[1,1] = system.ivc.P[1,1]*system.ivc.C + system.ivc.V0;
        elseif restart == "yes"
            svc = old["svc"];
            ivc = old["ivc"];
            system.svc.P[1,1] = svc["P"][end,1];
            system.svc.V[1,1] = svc["V"][end,1];
            system.svc.Q[1,1] = svc["Q"][end,1];
            system.ivc.P[1,1] = ivc["P"][end,1];
            system.ivc.V[1,1] = ivc["V"][end,1];
            system.ivc.Q[1,1] = ivc["Q"][end,1];
        end
    else
        if restart == "no"
            system.ivc.P[1,1] = 3*mmHgToPa;
            system.ivc.V[1,1] = system.ivc.P[1,1]*system.ivc.C + system.ivc.V0;
        elseif restart == "yes"
            ivc = old["ivc"];
            system.ivc.P[1,1] = ivc["P"][end,1];
            system.ivc.V[1,1] = ivc["V"][end,1];
            system.ivc.Q[1,1] = ivc["Q"][end,1];
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
# The BFGS update
# The inverse can either be updated with
#  H = (I-s*y'/y'*s)*H*(I-y*s'/y'*s)/(s'y)^2 - (H*y*s' + s*y'*H)/s'*y
# or
#  H = H + (s'*y + y'*H*y)*(s*s')/(s'*y)^2 - (B*y*s' + s*y'*H)/(s'*y)
#

struct BFGS{T1} <: QuasiNewton{T1}
   approx::T1
end
# function update!(scheme::BFGS, B::DirectApprox, Δx, y)
#    B.A = B.A + y*y'/dot(Δx, y) - B.A*y*y'*B.A/(y'*B.A*y)
# end
# function update!(scheme::BFGS, B::InverseApprox, Δx, y)
#    B.A = (I - Δx*y'/dot(y, Δx))*B.A*(I - y*Δx'/dot(y, Δx)) + Δx*Δx'/dot(y, Δx)
# end
#
function update(H, s, y, scheme::BFGS{<:InverseApprox})
   sy = dot(s, y)
   ρ = inv(sy)

   if isfinite(ρ)
      C = (I - ρ*s*y')
      H = C*H*C' + ρ*s*s'
   end

   H
end
function update(B, s, y, scheme::BFGS{<:DirectApprox})
   ρ = inv(dot(s, y))
   B + ρ*y*y' - B*s*s'*B/(s'*B*s)
end
function update!(H, s, y, scheme::BFGS{<:InverseApprox})
   sy = dot(s, y)
   ρ = inv(sy)

   if isfinite(ρ)
      Hy = H*y
      H .= H .+ ((sy+y'*Hy).*ρ^2)*s*s'
      Hys = Hy*s'
      H .= H .- (Hys + Hys').*ρ
   end
   H
end
# end
# function update!(H, s, y, scheme::BFGS{<:InverseApprox})
#    sy = dot(s, y)
#    ρ = inv(sy)
#
#    if isfinite(ρ)
#       C = (I - ρ.*s*y')
#       H .= C*H*C' + ρ*s*s'
#    end
#    H
# end
function update!(H, s, y, scheme::BFGS{<:DirectApprox})
    ρ = inv(dot(s, y))
    H .+= ρ*y*y' - H*(s*s')*H/(s'*H*s)
end

function update!(A::UniformScaling, s, y, scheme::BFGS{<:InverseApprox})
   update(A, s, y, scheme)
end
function update!(A::UniformScaling, s, y, scheme::BFGS{<:DirectApprox})
   update(A, s, y, scheme)
end


text:
using ParquetWriter

v = reduce(vcat, [readdir(p, join=true) for p in readdir("c:/data/parquet-test\\", join=true)])

@time adf = read_parquet("c:/scratch\\test.parquet");

files = [v for v in v if splitext(v)[2] == ".parquet" && isfile(v)]
@time Threads.@threads for f in files
    @time read_parquet(f, multithreaded = false);
end

files = readdir("C:/git/parquet-data-collection", join=true)

ff = [f for f in files if isfile(f) && splitext(f)[2]==".parquet"]

@time read_parquet.(ff)

@time read_parquet(ff[3]);

@time a= read_parquet(ff[3], multithreaded = false);

using BenchmarkTools
@benchmark read_parquet(ff[3])

ncols(ParFile(ff[3]))
filemetadata = metadata(ff[3])

ep = [@elapsed read_column(ff[3], filemetadata, i) for i in 1:169]

using StatsBase, Statistics
mean(ep)
std(ep)

ff[3]

using ParquetFiles

using ParquetFiles

@time read_parquet(ff[2], multithreaded =false);

for i in 1:ncols(ParFile(ff[2]))
    println(i)
    read_column(ff[2], i)
end

read_column(ff[2], 33)

a = metadata("c:/scratch/test.parquet")

a

using Parquet:ParFile
using Parquet
par = ParFile("c:/scratch\\test.parquet")
filemetadata = metadata("c:/scratch\\test.parquet")
filemetadata.row_groups[1].columns[5]
row_group = filemetadata.row_groups[1]
pages = Parquet.pages(par, row_group.columns[5])

@time a = read_parquet.(ff; multithreaded = false);

@time a = read_parquet.(ff; multithreaded = true);

a

for ff in ff
    try
        read_parquet(ff)
    catch
        println(ff)
    end
end

read_parquet("C:/git/parquet-data-collection/noshowappointments.parquet")

read_column("C:/git/parquet-data-collection/noshowappointments.parquet", 1)

using ParquetWriter

metadata.(files)

mf = metadata(files[2])

@time a = read_parquet(files[2])

using DataFrames|
@time a = read_parquet("c:/scratch/test.parquet");
@time ad = DataFrame(a, copycols=false);

read_column(files[1], 9)

read_parquet(files[1])

read_column(files[1], 2)

i = 0
for f in files
    global i
    try
        read_parquet(f)
        println("ok")
    catch
        i += 1
        println(f)
    end
end

metadata("c:/data/parquet-test/deep_solar_dataset/dsd50p.parquet")

m=metadata("c:/scratch/test1.parquet")

read_parquet("c:/scratch/test.parquet", [:V1])


text:
nanmean(x) = mean(filter(!isnan, x))
nanquantile(x,p) = quantile(filter(!isnan, x),p)

function quality_plot(file::String)
    qs = Vector{Vector{Float64}}(undef,0)
    reader = get_reader(file)
    for record in reader
        push!(qs,Vector{Float64}(quality(record)))
    end
    close(reader)
    q_array = float.(rstack(qs; fill=NaN))
    Iq_array = Int.(rstack(qs))
    max_len = maximum(length.(qs))
    q25,q50,q75 = Vector{Float64}(undef,max_len),Vector{Float64}(undef,max_len),Vector{Float64}(undef,max_len)
    Threads.@threads for i in 1:max_len
        #q25[i],q50[i],q75[i] = nanquantile(q_array[i,:],[0.25,0.5,0.75])
        q50[i] = nanmean(q_array[i,:])
    end
    rows,cols = size(q_array)
    hm_mat = zeros(max_len,41)
     @turbo for row in 1:rows
        for col in 1:cols
                hm_mat[row,Iq_array[row,col]+1] +=1
        end
    end
    fig = Figure()
    ax = Axis(fig[1,1])
    CairoMakie.heatmap!(ax,hm_mat[:,2:end], colormap = (:viridis,0.7))
    CairoMakie.lines!(ax,q50, color = :red)
    return fig
end

function quality_plot(file::String,ax)
    qs = Vector{Vector{Float64}}(undef,0)
    reader = get_reader(file)
    for record in reader
        push!(qs,Vector{Float64}(quality(record)))
    end
    close(reader)
    q_array = float.(rstack(qs; fill=NaN))
    Iq_array = Int.(rstack(qs))
    max_len = maximum(length.(qs))
    q25,q50,q75 = Vector{Float64}(undef,max_len),Vector{Float64}(undef,max_len),Vector{Float64}(undef,max_len)
    Threads.@threads for i in 1:max_len
        #q25[i],q50[i],q75[i] = nanquantile(q_array[i,:],[0.25,0.5,0.75])
        q50[i] = nanmean(q_array[i,:])
    end
    rows,cols = size(q_array)
    hm_mat = zeros(max_len,41)
     @turbo for row in 1:rows
        for col in 1:cols
                hm_mat[row,Iq_array[row,col]+1] +=1
        end
    end

    CairoMakie.heatmap!(ax,hm_mat[:,2:end], colormap = (:viridis,0.7))
    CairoMakie.lines!(ax,q50, color = :red)
end

function quality_plot(dir::String,identifier::String,out_path = "out.pdf";cols = 2, title_len = 28, title_from =false)

    file_list = []
    titles =[]
    for (root, dirs, files) in walkdir(dir)
        for (i,path) in enumerate(joinpath.(root, files))
            if contains(path,identifier)
            push!(file_list,path)
            push!(titles,files[i])
            end
        end
    end
    #n = length(file_list)
    tmps = []
    n = cols*3
    k = length(file_list) ÷ n
    for j in 1:k
        fig = Figure()
        set = (j-1)*n#+1
        xs = repeat(1:cols,(n+cols)÷cols)
        ys = vcat([fill(i,cols) for i in 1:((n+cols)÷cols)]...)
        for i in 1:max(n,length(file_list)-(k*n))
            try
            if title_from == false
                start = max(1,length(titles[i+set])-title_len)
                ax = Axis(fig[ys[i],xs[i]],title =  titles[i+set][start:end])
            else
                start = title_from
                ax = Axis(fig[ys[i],xs[i]],title =  titles[i+set][start:start+title_len])
            end
            quality_plot(file_list[i+set],ax)
            catch e
                println("Error encountered with file $(file_list[i+set])")
                println(e)
            end
        end
        save("tmpqualityplt$(j).pdf",fig)
        push!(tmps,"tmpqualityplt$(j).pdf")
    end
    read(`$(Poppler_jll.pdfunite()) $(tmps) $(out_path)`, String)# join temp pdfs together
    rm.(tmps); # remove temp pds after uniting them
    println("plots saved at $(out_path)")

end


text:
function parseKillerJSON(json_sums)
    sums = []
    for s in json_sums
        indices = Tuple[]
        for ind in s["indices"]
            push!(indices, tuple(ind...))
        end

        if haskey(s, "color")
            push!(sums, (result=s["result"], indices=indices, color=s["color"]))
        else
            push!(sums, (result=s["result"], indices=indices, color="white"))
        end
    end
    return sums
end

@testset "Killer Sudoku" begin

@testset "Killer Sudoku from wikipedia" begin
    com = CS.init()

    grid = zeros(Int, (9,9))

    com_grid = Array{CS.Variable, 2}(undef, 9, 9)
    for (ind,val) in enumerate(grid)
        com_grid[ind] = add_var!(com, 1, 9)
    end

    sums = parseKillerJSON(JSON.parsefile("data/killer_wikipedia"))

    for s in sums
        add_constraint!(com, sum([com_grid[CartesianIndex(ind)] for ind in s.indices]) == s.result)
    end

    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Solved
    @test fulfills_sudoku_constr(com_grid)
    for s in sums
        @test s.result == sum([value(com_grid[CartesianIndex(i)]) for i in s.indices])
    end
end

@testset "Killer Sudoku niallsudoku_5500 with coefficients" begin
    com = CS.init()

    grid = zeros(Int, (9,9))

    com_grid = Array{CS.Variable, 2}(undef, 9, 9)
    for (ind,val) in enumerate(grid)
        com_grid[ind] = add_var!(com, 1, 9)
    end

    sums = parseKillerJSON(JSON.parsefile("data/killer_niallsudoku_5500"))

    # the upper left sum constraint is x+y+z = 10 and the solution is 2+1+7
    # here I change it to 5*x+7*y+z = 24
    add_constraint!(com, 5*com_grid[CartesianIndex(1,1)]+7*com_grid[CartesianIndex(2,1)]+com_grid[CartesianIndex(2,2)] == 24)

    for s in sums[2:end]
        add_constraint!(com, sum([com_grid[CartesianIndex(ind)] for ind in s.indices]) == s.result)
    end

    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Solved
    @test fulfills_sudoku_constr(com_grid)
    @test 5*value(com_grid[CartesianIndex(1,1)])+7*value(com_grid[CartesianIndex(2,1)])+value(com_grid[CartesianIndex(2,2)]) == 24
    for s in sums[2:end]
        @test s.result == sum([value(com_grid[CartesianIndex(i)]) for i in s.indices])
    end
end

@testset "Killer Sudoku niallsudoku_5500 with negative coefficients and -9 to -1" begin
    com = CS.init()

    grid = zeros(Int, (9,9))

    com_grid = Array{CS.Variable, 2}(undef, 9, 9)
    for (ind,val) in enumerate(grid)
        com_grid[ind] = add_var!(com, -9, -1)
    end

    sums = parseKillerJSON(JSON.parsefile("data/killer_niallsudoku_5500"))

    # the upper left sum constraint is x+y+z = 10 and the solution is -2 + -1 + -7
    # here I change it to 5*x-7*y-z = 4
    add_constraint!(com, 5*com_grid[CartesianIndex(1,1)]-7*com_grid[CartesianIndex(2,1)]-com_grid[CartesianIndex(2,2)] == 4)

    for s in sums[2:end]
        add_constraint!(com, sum([com_grid[CartesianIndex(ind)] for ind in s.indices]) == -s.result)
    end

    add_sudoku_constr!(com, com_grid)

    @test solve!(com) == :Solved
    @test fulfills_sudoku_constr(com_grid)
    @test 5*value(com_grid[CartesianIndex(1,1)])-7*value(com_grid[CartesianIndex(2,1)])-value(com_grid[CartesianIndex(2,2)]) == 4
    for s in sums[2:end]
        @test -s.result == sum([value(com_grid[CartesianIndex(i)]) for i in s.indices])
    end
end

end


text:
function RealMatrixFormat()
    return RealMatrixFormat(())
end

function RealMatrixFormat(arg0::JString, arg1::JString, arg2::JString, arg3::JString, arg4::JString, arg5::JString)
    return RealMatrixFormat((JString, JString, JString, JString, JString, JString), arg0, arg1, arg2, arg3, arg4, arg5)
end

function RealMatrixFormat(arg0::JString, arg1::JString, arg2::JString, arg3::JString, arg4::JString, arg5::JString, arg6::NumberFormat)
    return RealMatrixFormat((JString, JString, JString, JString, JString, JString, NumberFormat), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function RealMatrixFormat(arg0::NumberFormat)
    return RealMatrixFormat((NumberFormat,), arg0)
end

function format(obj::RealMatrixFormat, arg0::RealMatrix)
    return jcall(obj, "format", JString, (RealMatrix,), arg0)
end

function format(obj::RealMatrixFormat, arg0::RealMatrix, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (RealMatrix, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function get_available_locales(::Type{RealMatrixFormat})
    return jcall(RealMatrixFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_column_separator(obj::RealMatrixFormat)
    return jcall(obj, "getColumnSeparator", JString, ())
end

function get_format(obj::RealMatrixFormat)
    return jcall(obj, "getFormat", NumberFormat, ())
end

function get_prefix(obj::RealMatrixFormat)
    return jcall(obj, "getPrefix", JString, ())
end

function get_real_matrix_format(::Type{RealMatrixFormat})
    return jcall(RealMatrixFormat, "getRealMatrixFormat", RealMatrixFormat, ())
end

function get_real_matrix_format(::Type{RealMatrixFormat}, arg0::Locale)
    return jcall(RealMatrixFormat, "getRealMatrixFormat", RealMatrixFormat, (Locale,), arg0)
end

function get_row_prefix(obj::RealMatrixFormat)
    return jcall(obj, "getRowPrefix", JString, ())
end

function get_row_separator(obj::RealMatrixFormat)
    return jcall(obj, "getRowSeparator", JString, ())
end

function get_row_suffix(obj::RealMatrixFormat)
    return jcall(obj, "getRowSuffix", JString, ())
end

function get_suffix(obj::RealMatrixFormat)
    return jcall(obj, "getSuffix", JString, ())
end

function parse(obj::RealMatrixFormat, arg0::JString)
    return jcall(obj, "parse", RealMatrix, (JString,), arg0)
end

function parse(obj::RealMatrixFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", RealMatrix, (JString, ParsePosition), arg0, arg1)
end


text:
#
# Function that computes the forces and potential energy,
# using a simple quadratic energy function
#
# modifies vector f, returns u

function uf_QUAD_LL!(n :: Int64, x :: Array{Float64}, f :: Array{Float64}, input :: MDInput)

  # Clear forces, ifirst, and inext arrays
  resetForces(f)
  @. f.ifirst = 0
  @. f.inext = 0

  # Add atoms to their boxes
  for i in 1:n
    ibox = trunc(Int64,x[i,1]/input.cutoff)+1
    jbox = trunc(Int64,x[i,2]/input.cutoff)+1
    inext[i] = ifirst[ibox,jbox]
    ifirst[ibox,jbox] = i
  end

voltar

  # Cycle over boxes computing interactions
  for iboxes in 1:nboxes

    ibox, jbox = setindex2D(iboxex,nboxes)
    i = ifirst[ibox,jbox]
    while i > 0

      # inside box
      celluf!(i,u,f,x,input,ibox,jbox)

      i = inext[i]
    end

  end

end

# Computes the interaction of atom i with atoms of box ibox, jbox

function celluf!(i,u,f,x,input,ibox,jbox,ifirst,inext)

  upair = 0.
  if ibox == jbox
    j = inext[i]
  else
    j = ifirst[ibox,jbox]
  end
  while j > 0

    xj = x[j,1] - x[i,1]
    yj = x[j,2] - x[i,2]
    xj = image(xj,input)
    yj = image(yj,input)

    r2 = xj^2 + yj^2
    if r2 > input.cutoff^2
      continue
    end
    r = sqrt(r2)

    if r < input.sig
      upartial[Threads.threadid()] = upartial[Threads.threadid()] + input.eps*(input.sig-r)^2
      drdxi = -xj/r
      drdyi = -yj/r
      fx = 2*input.eps*(input.sig-r)*drdxi
      fy = 2*input.eps*(input.sig-r)*drdyi
      f[i,1] = f[i,1] + fx
      f[i,2] = f[i,2] + fy
      f[j,1] = f[j,1] - fx
      f[j,2] = f[j,2] - fy
      j = inext[j]
    end

  end

end

  for ipair in 1:npairs

    j = (npairs%n)+1
    i = (npairs-j)/n+1

  Threads.@threads for i in 1:n-1

    for j in i+1:n

      xj = x[j,1] - x[i,1]
      yj = x[j,2] - x[i,2]
      xj = image(xj,input)
      yj = image(yj,input)

      r2 = xj^2 + yj^2
      if r2 > input.cutoff^2
        continue
      end
      r = sqrt(r2)

      if r < input.sig

        upartial[Threads.threadid()] = upartial[Threads.threadid()] + input.eps*(input.sig-r)^2

        drdxi = -xj/r
        drdyi = -yj/r

        fx = 2*input.eps*(input.sig-r)*drdxi
        fy = 2*input.eps*(input.sig-r)*drdyi

        f[i,1] = f[i,1] + fx
        f[i,2] = f[i,2] + fy
        f[j,1] = f[j,1] - fx
        f[j,2] = f[j,2] - fy

      end

    end

  end

  return sum(upartial)

end

function uf_QUAD!(n :: Int64, atoms :: Atoms, f :: Array{Float64}, input :: MDInput)

  x = atoms.x

  @. f = 0.
  upartial = zeros(Float64,Threads.nthreads())
  nenc_partial = zeros(Int64,Threads.nthreads())

  Threads.@threads for i in 1:n-1

    if atoms.status[i] == 2
      continue
    end

    for j in i+1:n

      if atoms.status[j] == 2
        continue
      end

      xj = x[j,1] - x[i,1]
      yj = x[j,2] - x[i,2]
      xj = image(xj,input)
      yj = image(yj,input)

      r2 = xj^2 + yj^2
      if r2 > input.cutoff^2
        continue
      end
      r = sqrt(r2)

      if r < input.sig

        upartial[Threads.threadid()] = upartial[Threads.threadid()] + input.eps*(input.sig-r)^2

        drdxi = -xj/r
        drdyi = -yj/r

        fx = 2*input.eps*(input.sig-r)*drdxi
        fy = 2*input.eps*(input.sig-r)*drdyi

        f[i,1] = f[i,1] + fx
        f[i,2] = f[i,2] + fy
        f[j,1] = f[j,1] - fx
        f[j,2] = f[j,2] - fy

      end

      atoms.status[i], atoms.status[j], encounter = update_status(atoms.status[i],atoms.status[j],r,input)
      if encounter
        nenc_partial[Threads.threadid()] = nenc_partial[Threads.threadid()] + 1
      end

    end

  end

  return sum(upartial), sum(nenc_partial)

end


text:
module TemporaryKeyFrostHelper

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/TemporaryKey" TemporaryKey(x::Integer, y::Integer, directory::String="collectables/FrostHelper/keytemp")

const placements = Ahorn.PlacementDict(
    "Temporary Key (Frost Helper)" => Ahorn.EntityPlacement(
        TemporaryKey
    ),
    "Temporary Key (With Return, Frost Helper)" => Ahorn.EntityPlacement(
        TemporaryKey,
        "point",
        Dict{String, Any}(),
        function(entity::TemporaryKey)
            entity.data["nodes"] = [
                (Int(entity.data["x"]) + 32, Int(entity.data["y"])),
                (Int(entity.data["x"]) + 64, Int(entity.data["y"]))
            ]
        end
    )
)

sprite = "collectables/FrostHelper/keytemp/idle00.png"

Ahorn.nodeLimits(entity::TemporaryKey) = length(get(entity.data, "nodes", [])) == 2 ? (2, 2) : (0, 0)

function Ahorn.selection(entity::TemporaryKey)
    x, y = Ahorn.position(entity)

    if haskey(entity.data, "nodes") && length(entity["nodes"]) >= 2
        controllX, controllY = Int.(entity.data["nodes"][1])
        endX, endY = Int.(entity.data["nodes"][2])

        return [
            Ahorn.getSpriteRectangle(sprite, x, y),
            Ahorn.getSpriteRectangle(sprite, controllX, controllY),
            Ahorn.getSpriteRectangle(sprite, endX, endY)
        ]

    else
        return Ahorn.getSpriteRectangle(sprite, x, y)
    end
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::TemporaryKey)
    px, py = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", Tuple{Int, Int}[])

    for node in nodes
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, px, py, nx, ny, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)
        px, py = nx, ny
    end
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::TemporaryKey, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
# SPDX-License-Identifier: MIT

"""
    struct TropicalSemiring{T<:AbstractFloat} <: Semiring
        val::T
    end

Log-semifield is defined as :
  * ``x \\oplus y \\triangleq \\text{max}(x, y)``
  * ``x \\otimes y \\triangleq x + y``
``\\forall x, y \\in \\mathbb{R}``.
"""
struct TropicalSemiring{T<:AbstractFloat} <: Semiring
    val::T
end

Base.:+(x::TropicalSemiring, y::TropicalSemiring) =
    TropicalSemiring(max(x.val, y.val))
Base.:*(x::TropicalSemiring, y::TropicalSemiring) =
    TropicalSemiring(x.val + y.val)
Base.zero(::Type{TropicalSemiring{T}}) where T = TropicalSemiring{T}(T(-Inf))
Base.one(::Type{TropicalSemiring{T}}) where T = TropicalSemiring{T}(T(0))
Base.isless(x::TropicalSemiring, y::TropicalSemiring) = isless(x.val, y.val)
Base.typemin(x::Type{TropicalSemiring{T}}) where T =
    TropicalSemiring{T}(typemin(T))
Base.typemax(x::Type{TropicalSemiring{T}}) where T =
    TropicalSemiring{T}(typemax(T))
Base.conj(x::Type{TropicalSemiring{T}}) where T = TropicalSemiring{T}(conj(x.val))


text:
using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using CSV
using FreqTables
using CategoricalArrays
using Tables
using StatsModels
cd("C:/Users/jakek/Documents/econ_phd/structural_modeling_course/forked_class/ProblemSets/PS3-gev/")

## 1.
df = CSV.read("nlsw88w.csv")
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occupation

# make cat 8 the base category.
function mlogit(alpha, X, d,Z)
	Zbar = Z .- Z[:,8]
	u = reduce(hcat,[X*alpha[(3*i-2):(3*i)] .+ Zbar[:,i] .*alpha[22] for i in 1:7])
	help = reduce(hcat,[d.==i for i in 1:7])
	loglike = sum(help .*u) - sum(log.(1 .+sum(exp.(u), dims=2)))
	return -loglike
end

# starting values - set intercepts to fraction present in data.
starting = Vector(prop(freqtable(df, :occupation)))
starting = repeat(starting[Not(8)], inner = 3)
starting[Not((1:7) .*3 .-2)] .= 0
# one extra row for the wage coef
starting = append!(starting, 0)

beta_hat_mlogit = optimize(b -> mlogit(b, X, y,Z),starting, BFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
res = reduce(hcat,[vec(beta_hat_mlogit.minimizer)[(3*i-2):(3*i)] for i in 1:7])
res = [res [vec(beta_hat_mlogit.minimizer)[22] 0 0]']
CSV.write("q1.csv", convert(DataFrame,res))
# gamma hat represents how much people value wages relative to non-wage amenities. it also represents how sensitive job choice is to relative wage

## 3.
# construct nests.
df.nest = (3 .* (df.occupation .>7) .+  Int.(df.occupation .<=3) .+ 2 .* (3 .< df.occupation .<=7) )

# make other (3) the base.
# denote the lambdas as the last coef in the set.

function nlogit(alpha, X, d,Z)
	Zbar = Z .- Z[:,8]
	u = reduce(hcat,[exp.((X*alpha[1:3].*(i<=3)+ X*alpha[5:7].*(3<i<=7).+ Zbar[:,i] .*alpha[9]) ./ (alpha[4]*(i<=3)+alpha[8]*(3<i<=7))) for i in 1:8])
	u[:, 8] .= 1
	wc_denom = sum(u[:,1:3], dims=2)
	bc_denom = sum(u[:,4:7], dims=2)
	u[:, 1:3] = (u[:, 1:3] .* wc_denom .^(alpha[4]-1)) ./ (1 .+ bc_denom.^alpha[8] .+wc_denom.^alpha[4])
	u[:, 4:7] = u[:, 4:7] .* bc_denom .^(alpha[8]-1) ./ (1 .+ bc_denom.^alpha[8] .+wc_denom.^alpha[4])
	u[:, 8] = u[:, 8] ./ (1 .+ bc_denom.^alpha[8] .+wc_denom.^alpha[4])
	u = log.(u)
	help = reduce(hcat,[d.==i for i in 1:8])
	loglike = sum( help .*u )
	return -loglike
end

# starting values - set intercepts to fraction present in data.
starting = Vector(prop(freqtable(df, :nest)))
starting = repeat(starting[Not(3)], inner = 4)
starting[Not((1:2) .*4 .-3)] .= 0
starting[[4 8]] .= 1
# one extra row for the wage coef
starting = append!(starting, 0)
beta_hat_nlogit = optimize(b -> nlogit(b, X, y,Z),starting, BFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
res = reduce(hcat,[vec(beta_hat_nlogit.minimizer)[(4*i-3):(4*i)] for i in 1:2])
res = [res [vec(beta_hat_nlogit.minimizer)[9] 0 0 0]']
CSV.write("q3.csv", convert(DataFrame,res))


text:
"""
    Beta
* Parametrization
* Score
* Fisher Information
* `time_varying_params` map.
* Default link
"""
Beta

function score!(score_til::Matrix{T}, y::T, ::Type{Beta}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = log(y) + digamma(param[t, 1] + param[t, 2]) - digamma(param[t, 1])
    score_til[t, 2] = log(1 - y) + digamma(param[t, 1] + param[t, 2]) - digamma(param[t, 2])
    return
end

function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{Beta}, param::Matrix{T}, t::Int) where T
    minus_trigamma_a_b = -trigamma(param[t, 1] + param[t, 2])
    aux.fisher[1, 1] = trigamma(param[t, 1]) + minus_trigamma_a_b
    aux.fisher[2, 2] = trigamma(param[t, 2]) + minus_trigamma_a_b
    aux.fisher[2, 1] = minus_trigamma_a_b
    aux.fisher[1, 2] = minus_trigamma_a_b
    return
end

function log_likelihood(::Type{Beta}, y::Vector{T}, param::Matrix{T}, n::Int) where T
    loglik = 0.0
    for t in 1:n
        loglik += (param[t, 1] - 1)*log(y[t]) + (param[t, 2] - 1)*log(1 - y[t]) - logbeta(param[t, 1], param[t, 2])
    end
    return -loglik
end

# Links
function link!(param_tilde::Matrix{T}, ::Type{Beta}, param::Matrix{T}, t::Int) where T
    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
    param_tilde[t, 2] = link(LogLink, param[t, 2], zero(T))
    return
end
function unlink!(param::Matrix{T}, ::Type{Beta}, param_tilde::Matrix{T}, t::Int) where T
    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
    param[t, 2] = unlink(LogLink, param_tilde[t, 2], zero(T))
    return
end
function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{Beta}, param::Matrix{T}, t::Int) where T
    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
    aux.jac[2] = jacobian_link(LogLink, param[t, 2], zero(T))
    return
end

# utils
function update_dist(::Type{Beta}, param::Matrix{T}, t::Int) where T
    small_threshold!(param, SMALL_NUM, t)
    return Beta(param[t, 1], param[t, 2])
end

function params_sdm(d::Beta)
    return Distributions.params(d)
end

function num_params(::Type{Beta})
    return 2
end


text:
#=
# ЛЕКЦИЯ 4

## 1. Система типов Julia

- Иерархия типов
Во главе всей иерархии типов - тип Any

Тип Number - стоит во главе всех числовых типов (и является подтипом Any)
У типа Number есть всего 2 подтипа:
    Complex
    Real

Подтипами типа Real, например, являются типы:
    AbstractFloat # - тип всех чисел с плавающей точкой
    AbstractIrrational # - тип всех иррациональных чисел
    Integer # - тип всех целых чисел
    Rational # - тип всех рациональных чисел

Подтипами типа Integer, например, являются типы:
    Bool # - логический тип
    Signed # - знаковые целые
    Unsigned # - беззнаковые целые

Подтипами типа Signed, например, являются типы:
    BigInt # - "длинные" целые со знаком
    Int128 # - 128 битные целые со знаком
    Int64  # - 64 битные целые со знаком - синоним Int
    Int32  # - 32 битные целые со знаком
    Int16  # - 16 битные целые со знаком
    Int8   # - 8 битные целые со знаком

Типы
    BigInt, Int128, Int64, Int32, Int16, Int8

- это примеры КОНКРЕНТНЫХ типов.

Все остальные рассмотренные типы были АБСТРАКТНЫМИ.

Конкретные типы - это типы, стоящие на самой нижней иерархической ступеньке.
Все остальные типы - это абстрактные типы.

Любой объект (переменная или значение) имеет какой-либо конкретный тип.
Создать объект абстрактного типа невозможно, т.к. абстрактный тип не предполагает
наличия каких-либо конкретных данных, абстрактные типы служат исключительно для создания
иерархии типов.

Некоторые из абстрактных типов являются параметрическими.
Например, параметрическим является абстрактный тип Complex:

    Complex{Float64} # - это уже получится конктетный тип (т.к. в данном случае параметр - есть конкретный тип Float64)
    Complex{Real} # - а это по-прежнему абстрактный тип (т.к. в данном случае параметр - есть абстрактный тип)

Конкретные типы не могут иметь подтипов (т.к. они стоят на САМОЙ НИЖНЕЙ иерархической ступеньке).

Для выяснение соподчиненности типов имеются встроенные функции:

    supertype(TypeName) # - возвращает родительский тип для типа TypeName
    subtypes(TypeName) # - возвращает вектор (список) дочерних типов для типа TypeName

Кроме того, имеется бинарная операция <: (>:), например:

    Int64 <: Integer # - вернёт true
    Int64 >: Integer # - вернёт false
    Int64 <: Real # - тоже вернёт true

Система типов Julia является расширяемой, т.е. в иерархию изначально существующих (встроенных) типов.
Новые типы добавляемые в систему типов, называются пользовательскими.

Пользовательские типы могут быть как конкретными (определяются на основе конструкции struct),
так и абстрактными (abstract type), в том числе и параметрическими.

Все стандартные контейнерные типы данных, такие как Vector, или Matrix (оба этих типа являются просто частным
случаем типа Array) являются параметрическими. Например, определять вектор типа Int (т.е. Int64) надо так:

    Vector{Int}

Но, ВНИМАНИЕ, определение:
    Vector{Integer}
- это было бы неправильным определением типа, правильно в этом случае было бы так:
    Vector{<:Integer}
Для этого правила есть объективные причины.

## 2. Обобщенный подход к программированию

Подход к программированию (парадигма программирования), согласно которому описание данных и алгоритмов (функций)
делается так, что бы это описание было применимо ни к какому-то одному конкретному типу данных
(или к набору конкретных типов данных), а было бы применимо к различным типам данных, причем для этого
не требовалось бы вносить никакие изменения в это описание.

Обобщенное описание типов данных сводится к определению абстрактных типов, в том числе к определению параметрических
типов. Но к этому мы ещё вернёмся.

Рассмотрим пример проектирвания обобщённой функции.

Пусть требуется написать функцию, возвращающую определитель матрица 2*2.
Если элементы матрицы имеют тип Int, то определение функция может иметь вид:
=#
function det(A::Matrix{Int})
    @assert size(A) == (2,2) # встроенная функция size возвращает кортеж из размеров массива
    # - это проверка условия, что матрица размера 2*2
    # (в противном случае произойдет прерывание функции и вывод сообщения об ошибке)
    return A[1,1]*A[2,2]-A[1,2]*A[2,1]
end

#=
А если теперь понадобилось бы вычислить определитель матрицы, элементы которой имеют тип Float64,
то понадобилось бы новое определение функции:
=#
function det(A::Matrix{Float})
    @assert size(A) == (2,2)
    # - это проверка условия, что матрица размера 2*2
    # (в противном случае произойдет прерывание функции и вывод сообщения об ошибке)
    return A[1,1]*A[2,2]-A[1,2]*A[2,1]
end

#=
Как видим, тело новой функции осталось прежним, изменилась только аннотация типа в заголовке функции.

При этом, если понадобится рассмотреть матрицы каких-либо еще других типов, например,
Rational{Int} или Comple{Float64}, то понадобятся и новые определения функции. При этом по внешнему виду
тело функции будет оставаться всё тем же.

Для того, что бы не плодить большое число описаний этой функции, её описание можно было бы сделать
обобщенным (generic).
Для этого можно было бы просто отказаться от аннотации типа её аргумента:
=#
function det(A) # det(A::Any) - это то же самое
    @assert size(A) == (2,2)
    A[1,1]*A[2,2]-A[1,2]*A[2,1] # слово return здесь не обязательно
end

#=
Если интерфейс фактического тип параметра A включает опреации `*`, `-`, то данная функция может быть вычислена.
Например:
=#
> det([1 2;3 4]) # -> -2 (тип Int64)
> det([1.0 2.0; 3.0 4.0])  # -> -2.0 (тип Float64)

#=
Но в следующем случае:
=#
> det(["1" "2"; "3" "4"])
#=
произойдет ошибка, вызванная тем, что опрерация `-` не определена для типа String
(операция `*` - в Julia это операция конкатенации строк, так что резултатом A[1,1]*A[2,2] будет "14",
а результатом A[1,2]*A[2,1] будет "23")
=#

#=
Если тип аргумента обобщенной функции аннотирован как Any (в этом случае можно сказать, что тип не аннотирован),
то из заголовка такой функции не ясно, каким интерфейсом этот аргумент должен обладать.
В этом случае говорят, что требуемый интерфейс явно не описан.
Такую типизацию аргументов функции, без явных требований к интерфейсу, принято называть "утиной" типизацией.

Такое название происходит от образного выражения "если нечто ходит как утка и крякает как утка, то это нечто
и есть утка".

Применительно к нашему примеру, если для фактического типа аргумента определены операции `*` и `-`, то этот
аргумент является допустимым (является той самой "уткой").

Однако мы могли бы явно описать требования к интерфейсу аргументу нашей функции, определив его тип так:
=#
function det(A::Vector{<:Number}) # det(A::Any) - это то же самое
    @assert size(A) == (2,2)
    A[1,1]*A[2,2]-A[1,2]*A[2,1] # слово return здесь не обязательно
end
#=
Здесь мы явно указали, что элементами матрицы должны быть числа, а для любых чисел должны быть определены
операции `*` и `-`
=#

#=
С практической точки зрения разница между утиной типизацией и явным указанием типа аргумена Matrix{:<Number}
не столь уж велика.
Так в первом случае при вызове
=#
> det(["1" "2"; "3" "4"])
#=
ошибка произойдет при попытке скомпилировать выражение A[1,1]*A[2,2]-A[1,2]*A[2,1]

А во втором случае ошибка произойдет при попытке скомпилировать выражение det(["1" "2"; "3" "4"]).

Иногда, правда второй случай более предпочтителен в виду того, что в противном случае ошибка при попытке
компиляции может происходить где-то достаточно "глубоко" в цепочке вызовов различных функций, и это может
несколько затруднять отладку.

С другой стороны, иногда, при использовании сложных типов, их полное описаение может оказаться довольно
громоздким, и  в таких случаях его часто упрощают (делают более общим, чем требуется) либо вообще оставляют
соответствующий аргумент функции без аннотации.

### Замечание

Вместо описания
=#
det(A::Matrix{<:Number}) = A[1,1]*A[2,2]-A[1,2]*A[2,1]
#=
В языке Julia возможно эквивалентное опивание:
=#
det(A::Matrix{T}) where {T <:Number} = A[1,1]*A[2,2]-A[1,2]*A[2,1]
#=
Здесь разница лишь в том, что теперь в теле функции есть возможность использовать тип T.
В рассмотренном примере нам это в общем-то ни к чему, поэтому рассморим другой простой пример.
Пусть требуется написать обобщенную функцию, суммирующую элементы вектора
=#
function sum(vector::Vector{T}) where{T<:Number}
    s=T(0)
    # здесь 0 типа Int64 преобразуется к нулю фактического типа T
    #(иначе переменная s будет иметь не стабильный тип, что плохо, с точки зрения генерации эффетивного кода
    # компилятором Julia)
    for a in vector
        s =+ a
    end
    return s
end
#=
В противном случае пришлось бы делать так:
=#
function sum(vector::Vector{<:Number})
    s=eltype(vector)(0)
    # встроенная функция eltype(vector) возвращает тип элементов данного массива vector
    for a in vector
        s =+ a
    end
    return s
end

#=
В действительности функция sum является встроенной и её возможная реализация приведена просто в качестве примера.

## 3. Пользовательские типы данных

Конкретные (не абстрактные) пользовательские типы проектируются на основе конструкции struct.

С помощью этой конструкции определяются составные типы данных, состоящие из некоторого количества полей
(формально число полей в структуре может быть любым, в том числе и равным 0).
Например:
=#

struct NameType
    a::T1
    b::T2
end

#=
- это определение нового типа данных с именем NameType (имя типа может быть практически любым,
но следует придерживаться указанного стиля именования).
Поля этого составного типа данных имеют типы T1, Т2 соответственно.
Типы полей могут быть любыми, но желательно, чтобы они были конкретными, а не абстрактными, что обеспечит
возможность компилятору Julia генерировать более эффективный код.

Что бы создать объет данного типа нужно воспользоваться соответствующим конструктором типа. Конструктор типа
- это специальная функция, которая создает объет соответствующего типа и возвращает ссылку на этот объект
(именно - ссылку, а не сам объект). Имя конструктрора всегда совпадает с именем типа
(так же как и для встроенных типов).

Если конструктор не был специально определен, то всегда имеется конструктор по умолчанию.
В данном случае конструктором по умолчанию можно было бы воспользоваться так:
=#
> x=NameType(0,0)
#=
Тогда, при условии, что значение 0 может быть приведено к типам T1 и T2,
поле `a` объекта `x` получит значение T1(0), а поле `b` - значение T2(0).

Получить значения этих полей можно, используя составное имя, так:
=#
> x.a
> x.b
#=
Но изменить значения этих аолей не получится, например:
=#
> x.a = T1(1)
#=
Изменять значения полей структуры в данном случае нельзя, потому что структура была объявлена (не явно),
как неизменяемая (с точки зрения возможности компиляции в эффективный код неизменяемые структуры
предпочтительнее изменяемых).
Для создания изменяемых структур используется специальное слово `mutable`, например:
=#
mutable struct NameType2
    a::Int
    b::Int
end

> y=NameType2(0,0)

> y.a=1 # такое присваивание возможно

> y=NameType1(0,0)

> y.a=1 # это присваивание не возможно

#=
Перейдем теперь к более содержательным примерам.

### Пример 1

Пусть Робот на ходится где-то на поле между внутренними перегородками или прямолинейной формы или в форме
прямоугольников (внутренние перегородки не касаются друг друга и внешней перегородки).

И пусть требуется сначала переместить Робота в один из углов поля, затем выполнить некоторые действия
и вернуться в исходный угол, а затем от туда вернуться в исходное положение.

С такой задачей мы уже встречались ни один раз. Вот как можно унифицировать эту задачу, спроектировав специальный
тип данных, в котором будет представлен путь возвращения Робота из стартового угла в исходное положение.
=#
struct BackPath
    sides::NTuple{2,HorizonSide}
    path::Vector{Int}
# Определим теперь специальный конструктор для этого типа (конструктор по умочанию нас не устроит):

    function BackPath(robot, sides::NTuple{2,HorizonSide}=(Sud,West))
        local path=Int[]
        while !isborder(robot,sides[1]) || !isborder(robot,sides[2])
            for s in sides
                push!(path, movements!(robots,s))
            end
        end
        return new(reverse(inverse.(sides)), reverse(path))
    end
end
#=
При определении конструктора всегда используется встроенная фунуция `new`, у которой должно быть
столько же аргументов, сколько имеется полей в структуре, и в результате выполнения котрой будет создан
объект соответствующего типа, полям которого будут присвоены значения аргументов функции new: первому полю
присваивается значение первого аргумента, второму полу - значение вторго аргумента, и т.д.  по порядку.

Так же для данного типа определим ещё следующую функцию, с помощью которой Робот может быть перемещен
в исходное положение путем прохождения ранее пройденного пути в обратном направлении:
=#
function back!(robot, backpath::BackPath)
    i=1
    for n in backpath.path
        movements!(robot, backpath.sides[i], n)
        i = i%2 + 1
    end
end

md"""
Общая схема решения задач с использованием объекта типа BackPath:
"""
#УТВ: Робот в некоторой начальной клетке среди прямоугольных перегородок
    back_path = BackPath(r)
    #УТВ: Робот - в Ю-З углу
    #..... требуемые действия
    #УТВ: Робот снова находится в Ю-З углу
    back!(r, back_path)
    #УТВ: Робот - в исходном положении

md"""
### Пример 2

Для решения некоторых задач может понадобится знать Декартовы координаты положения Робота относительно некоторого
"начального" его положения.
Для этого спроектируем сначала следующий тип данных на основе  изменяемой (мутабельной) структуры:
"""
mutable struct Coord
    x::Int
    y::Int
    Coord() = new(0,0) # - такой конструктор определен для удобства использования
    Coord(x::Int,y::Int) = new(x,y)
end
md"""
Конструкторов может быть несколько, и конструкторы можно определять внутри конструкции struct (все остальные функции
для типа) должны определяться как внешние функции, т.е. за пределами конструкции struct.

Так же для данного типа определим ещё слеующие две функции. Первая из них `move!` изменяет значение текущих координат
при ВООБРАЖАЕМОМ перемещении Робота в заданном направлении на 1 шаг.

Вторая `get_coord` - возвращает кортеж с текущими координатами.
"""
function move!(coord::Coord, side::HorizonSide)
    if side==Nord
        coord.y += 1
    elseif side==Sud
        coord.y -= 1
    elseif side==Ost
        coord.x += 1
    else #if side==West
        coord.x -= 1
    end
end

get_coord(coord::Coord) = (coord.x, coord.y)

md"""
Наличие функции `get_coord` избавляет от необходимости помнить (знать) внутренне устройство объекта типа
структуры `Coord`.

Эти две функции `move!` и `get_coord` составляют интерфейс типа `Coord`.

Взаимодействие с объектами данного типа только посредством его интерфейса,
исключая прямое взаимодействие с полями объекта, есть соблюдение принципа ИНКАПСУЛЯЦИИ
данных (принцип сокрытия данных в объете).

Это очень важный принцип и ему слеует следовать.

Например, использовать разработанный тип можно так:
"""
> coord=Coord()
> move!(coords,Nord)
> get_coords(coords) # -> (0,1)

md"""
## 4. Иерархия пользовательских типов

Для того чтобы получить Робота, который при своих перемещениях будет отслеживать свои координаты
потребуется  "скрестить" два типа:

Robot + Coord

T.к. оба типа являются конкретными, то сделать это можно следующими способами: путем композиции типов или путем
проектирования дополнительного абстрактного типа с последующим наследованием от него.

### Композиция типов

Следующая структура представляет собой композицию конкретного типа Robot и конкретного типа Coord:
"""
struct XYRobot
    robot::Robot
    coord::Coord
end

md"""
остается только определить интерфейс для этого типа
(чем мы и займемся в следующий раз)
"""


text:
module IsingModelTest

    export IsingModel, Ensemble, Plot2D, Gif2D, GetUpSpins, CentralMoments, EnsProbUpSpins

    using Colors, Plots, StatsBase, Parameters, Distributions, StatsBase, LinearAlgebra;

    """
    Define the Ising model struct.
    """
    @with_kw struct IsingModel
        N::Int64 = 100
        h::Float64
        J::Float64
        p0::Float64 = 0.5 # initial prob of spin up
        periodic::Int64 = 1
        mf::Int64 = 0
        fully_connected::Int64 = 0
    end

    """
    Define the struct for an ensemble of sims.
    """
    @with_kw struct Ensemble
        IM::IsingModel
        e_size::Int64 = 1 # size of ensemble
        method::String = "Glauber" # or "SSA"
        tₘ::Real = 10^5 # no units of time
        Δt::Real = 10^2 # no units of time
        ts::Real = 1.0 # timescale: real time τ = ts*t
        sims = SimRun(method,e_size,tₘ,Δt,ts,IM)
    end

    """
    Choice of algorithm
    """
    function SimRun(alg::String,e_size::Int64,tₘ::Real,Δt::Real,ts::Real,IM::IsingModel)
        if alg == "Glauber"
            return Glauber(e_size,tₘ,Δt,ts,IM)
        elseif alg == "SSA"
            return SSA(e_size,tₘ,Δt,ts,IM)
        else
            error(alg, " is not an appropriate choice of algorithm.
            Please choose either 'Glauber' or 'SSA'.")
        end
    end

    """
    Make the initial NxN matrix of spins.
    """
    function S_make(p0::Float64, N::Int64)
        if p0 > 1 || p0 < 0
            error("p must be bounded between 0 and 1")
        else
            return ones(N,N)-2*sample([1,0],Weights([1-p0,p0]),(N,N))
        end
    end

    """
    Define the effective field at each spin.
    """
    function hₑ(i::Int64, j::Int64, S::Matrix{Float64}, IM::IsingModel, periodic::Int64)
        @unpack J, h, N, fully_connected = IM;
        if fully_connected == 1
            return 2*h + 2*(J/N)*(sum(S)-S[i,j]) # J scaled by N.
        elseif fully_connected == 0 && periodic == 0
            if i==1 && j==1
                return 2*h + 2*J*(S[i+1,j]+S[i,j+1])
            elseif i==N && j==N
                return 2*h + 2*J*(S[i-1,j]+S[i,j-1])
            elseif i==1 && j==N
                return 2*h + 2*J*(S[i+1,j]+S[i,j-1])
            elseif i==N && j==1
                return 2*h + 2*J*(S[i-1,j]+S[i,j+1])
            elseif i==1
                return 2*h + 2*J*(S[i+1,j]+S[i,j+1]+S[i,j-1])
            elseif j==1
                return 2*h + 2*J*(S[i+1,j]+S[i-1,j]+S[i,j+1])
            elseif i==N
                return 2*h + 2*J*(S[i-1,j]+S[i,j+1]+S[i,j-1])
            elseif j==N
                return 2*h + 2*J*(S[i+1,j]+S[i-1,j]+S[i,j-1])
            else
                return 2*h + 2*J*(S[i+1,j]+S[i-1,j]+S[i,j+1]+S[i,j-1])
            end
        else # is periodic BC
            return 2*h + 2*J*(S[mod1(i+1,N),j]+S[mod1(i-1,N),j]+S[i,mod1(j+1,N)]+S[i,mod1(j-1,N)])
        end
    end

    """
    Effective field in MF setting.
    """
    function hₘ(i::Int64, j::Int64, S::Matrix{Float64}, IM::IsingModel, periodic::Int64)
        @unpack J, h, N, fully_connected = IM;
        m = mean(S);
        Jm = J*m;
        if fully_connected == 1
            return 2*h + 2*J*m # J scaled by N.
        elseif fully_connected == 0 && periodic == 0
            if i==1 && j==1
                return 2*h + 2*2*Jm
            elseif i==N && j==N
                return 2*h + 2*2*Jm
            elseif i==1 && j==N
                return 2*h + 2*2*Jm
            elseif i==N && j==1
                return 2*h + 2*2*Jm
            elseif i==1
                return 2*h + 2*3*Jm
            elseif j==1
                return 2*h + 2*3*Jm
            elseif i==N
                return 2*h + 2*3*Jm
            elseif j==N
                return 2*h + 2*3*Jm
            else
                return 2*h + 2*4*Jm
            end
        else # is periodic BC
            return 2*h + 2*4*Jm
        end
    end

    """
    Energy difference calculation
    """
    function ΔH(i::Int64, j::Int64, S::Matrix{Float64}, IM::IsingModel, periodic::Int64=1, mf::Int64=0)
        if mf == 0
            Δe = hₑ(i, j, S, IM, periodic)*S[i,j]
        elseif mf == 1
            Δe = hₘ(i, j, S, IM, periodic)*S[i,j]
        else
            error("mf must be 1 or 0.")
        end
        return Δe
    end

    """
    Logistic rule.
    """
    function LogisticRule(E₀::Float64)
        return (1+exp(E₀))^-1
    end

    """
    Glauber simulation.
    """
    function Glauber(ens_its::Int64,tₘ::Int64,Δt::Int64,μ₁::Float64,IM::IsingModel)
        @unpack p0, N, mf, periodic = IM;
        ens_store = Array{Vector{Matrix{Float64}}}(undef, ens_its);
        for it in 1:ens_its
            S = S_make(p0, N); # make initial config
            S_temp = copy(S); # use S_temp for updating.
            store = Array{Matrix{Float64}}(undef, tₘ÷Δt + 1);
            store[1] = copy(S_temp);
            for n in 1:tₘ+1
                i,j = Int.(map(round,rand(2).*(N-1))) .+ [1,1]; # agent chosen at random
                E₀ = ΔH(i,j,S_temp,IM,periodic,mf) # periodic BC
                pr = LogisticRule(E₀);
                weights = [pr,1-pr]; # flip, don't flip
                flip = sample([1,0],Weights(weights));
                if flip == 1 # if chosen to flip perform
                    S_temp[i,j] = -S_temp[i,j]
                end
                if n%Δt == 0 # store every Δt steps
                    store[(n÷Δt)] = copy(S_temp)
                end
            end
            ens_store[it] = store
        end
        return ens_store
    end

    """
    Propensity for SSA
    """
    function propensity(S::Matrix{Float64},IM::IsingModel,γ::Float64)
        @unpack N, mf, periodic = IM;
        γ = γ/(N^2); # normalise the timescale to match Glauber.
        props = zeros(size(S))
        for i in 1:size(S)[1]
            for j in 1:size(S)[2]
                E₀ = ΔH(i,j,S,IM,periodic,mf); # periodic BC
                pr = LogisticRule(E₀);
                props[i,j] = γ*pr;
            end
        end
        return props
    end

    """
    SSA simulations
    """
    function SSA(ens_its::Int64,τₘ::Real,dτ::Real,γ::Real,IM::IsingModel)
        @unpack p0, N, mf, periodic = IM;
        ens_store = Array{Vector{Matrix{Float64}}}(undef, ens_its);
        for it in 1:ens_its
            S = S_make(p0, N); # make initial config
            S_temp = copy(S); # use S_temp for updating.
            store = Array{Matrix{Float64}}(undef, Int(τₘ÷dτ)+1);
            τ = 0.0;
            sim = 0;
            while τ < τₘ
                props = propensity(S_temp,IM,γ);
                props_flat = vec(transpose(props)); # concatted rows
                w0 = sum(props);
                r1, r2 = rand(2);
                u = (1/w0)*log(1/r1);
                k = findfirst(x -> x>=r2*w0,cumsum(props_flat));
                n = 0;
                while τ+n*dτ < τ+u
                    if τ+n*dτ > τₘ
                        store[Int(τₘ÷dτ)+1] = copy(S_temp)
                        break
                    else
                        store[Int((τ+n*dτ)÷dτ)+1] = copy(S_temp)
                        n+=1;
                    end
                end
                new_i = ((k-1)÷N)+1;
                new_j = k-(new_i-1)*N;
                S_temp[new_i,new_j] = -S_temp[new_i,new_j];
                τ += u;
                if τ > τₘ
                    store[Int(τₘ÷dτ)+1] = copy(S_temp)
                end
                sim += 1;
                if sim%1000==0
                    println(τ)
                end
            end
            ens_store[it] = store;
        end
        return ens_store
    end

    """
    Plot slice at some time.
    """
    function Plot2D(t::Int64, m::Int64, ens::Ensemble) # m = sim set
        @unpack method, Δt, tₘ, ts, sims = ens;
        if method == "Glauber" # then t is timeslice
            plot(Gray.(sims[m][t]))
        elseif method == "SSA" # then t is the actual time
            n = Int((t/ts)÷Δt)
            if n <= tₘ
                plot(Gray.(sims[m][n]))
            else
                error("Plot time must be within the simulated times.")
            end
        else
            error("Simulation method chosen is invalid.")
        end
    end

    """
    IM gif maker.
    """
    function Gif2D(m::Int64, ens::Ensemble, fname::String, frames::Real=40)
        @unpack sims = ens;
        anim = @animate for i ∈ 1:length(sims[m])
            plot(Gray.(sims[m][i]))
        end
        return gif(anim, join([fname,".gif"]), fps = frames)
    end

    """
    Take the trajectories of up spins from each set of sims in an ensemble.
    """
    function GetUpSpins(ens::Ensemble, norm::Int64=0)
        @unpack IM, sims = ens
        @unpack N = IM
        store = Array{Float64,2}(undef, length(sims), length(sims[1]));#Array{Array{Float64,1}}
        for m in 1:length(sims)
            for t in 1:length(sims[1])
                store[m,t] = count(j->(j==1.0),sims[m][t])
            end
        end
        if norm == 0
            return store
        elseif norm == 1
            return store/N^2
        else
            error("norm must take a value of 0 or 1.")
        end
    end

    """
    Get the central moments from the ensemble
    """
    function CentralMoments(ens::Ensemble,n::Int64,norm::Int64=0)
        upspins = GetUpSpins(ens, norm)
        store = Array{Float64,1}(undef, length(upspins[1,:]))
        if n == 1
            for t in 1:length(upspins[1,:])
                store[t] = mean(upspins[:,t])
            end
        else
            for t in 1:length(upspins[1,:])
                store[t] = moment(upspins[:,t],n)
            end
        end
        # need to sort out times
        @unpack Δt, ts = ens;
        no_unit_t = LinRange(0:Δt:length(upspins[1,:])-1);
        times = no_unit_t * ts;
        return (times, store)
    end

    """
    Get prob distribution at all times for an ensemble
    """
    function EnsProbUpSpins(ens::Ensemble, T::Int64)
        UpSpins = GetUpSpins(ens)
        mod_bins = LinRange(minimum(UpSpins[:,T])-0.5:1:maximum(UpSpins[:,T])+0.5);
        mid_pts = LinRange(minimum(UpSpins[:,T]):1:maximum(UpSpins[:,T]));
        bin_vals = normalize(fit(Histogram, UpSpins[:,T], mod_bins), mode=:probability).weights;
        return (mid_pts, bin_vals)
    end

end # module end

max_stars_repo_name: jamesholehouse/Ising

max_issues_repo_name: jamesholehouse/Ising

max_forks_repo_name: jamesholehouse/Ising


text:
function dc_networkflow(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, PTDF::PTDFArray)

    fbr = m[:fbr]
    name_index = m[:fbr].axes[1]
    time_index = m[:fbr].axes[2]
    nodes = PTDF.axes[2].val
    node_set = 1:length(nodes)

    ts_flow = gemm('N', 'N', PTDF.data, netinjection.timeseries_active)

    remove_undef!(netinjection.var_active)

    branchflow = JuMP.JuMPArray(Array{ConstraintRef}(undef, length(name_index), time_index[end]), name_index, time_index)

    #TODO: Make consistent with the use of AxisArrays. This syntax doesn't exploit it properly
    for t in time_index
        for (ix,branch) in enumerate(name_index)
            branchflow[branch,t] = @constraint(m, fbr[branch,t] == sum(netinjection.var_active[i,t]*PTDF.data[ix,i] for i in node_set) - ts_flow[ix,t])
        end
    end

    JuMP.register_object(m, :branchflow, branchflow)

    return m
end


text:
# linear function with [-0.01, 0.01] uniform distribution
function myLinear(T::Type, indim::Int, outdim::Int)
    r = T(0.01)
    w = rand(-r, r, outdim, indim)
    b = fill(T(0), outdim, 1)
    Linear(zerograd!(Var(w)), (zerograd!(Var(b))))
end

function scale!(embed::Embedding)
    mean = sum(w -> sum(w.data), embed.ws)
    std = sum(w -> dot(w.data, w.data), embed.ws)
    count = length(embed.ws) * length(embed.ws[1].data)
    mean = mean / count
    std = sqrt(std / count - mean^2)
    info("Scaling word embeddings:")
    info("(mean = $mean, std = $std) -> (mean = 0.0, std = 1.0)")
    for w in embed.ws
        w.data -= mean
        w.data /= std
    end
end

function update!(opt::Union{SGD,MyAdaGrad}, gold::Vector{Int}, pred::Var)
    l = crossentropy(gold, pred) * (1 / length(gold))
    loss = sum(l.data)
    vars = gradient!(l)
    for v in vars
        if isa(v.f, Merlin.Functor)
            # l2 normalization
            BLAS.axpy!(Float32(10e-8), v.data, v.grad)
            Merlin.update!(v.f, opt)
        end
    end
    loss
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
"Gets a tuple of the types in a Union."
@inline union_types(T) = (T, )
@inline function union_types(u::Union)
    # Supposedly, 'u.a' is the first type and 'u.b' is a Union of the other types.
    # However, sometimes a is the Union, and b is the type.
    # So I pass both sides through 'union_types()' to be sure.
    (union_types(u.a)..., union_types(u.b)...)
end

"
When deserializing a union of types, this determines the order that the types are tried in.
Lower-priority values are tried first.
"
union_ordering(T::Type)::Float64 =
    # If this is an AbstractType, give it a specific priority.
    # Otherwise, unknown/user-made types are given the lowest priority.
    (StructTypes.StructType(T) isa StructTypes.AbstractType) ?
        union_ordering(Val(:abstract)) :
        +Inf
#
union_ordering(::Type{<:Enum}) = 0.0
union_ordering(::Type{<:Integer}) = 1.0
union_ordering(::Type{<:AbstractFloat}) = 2.0
union_ordering(::Type{<:Dates.TimeType}) = 3.0
union_ordering(::Type{Symbol}) = 4.0
union_ordering(::Type{<:AbstractString}) = 5.0
union_ordering(::Val{:abstract}) = 6.0
#TODO: Make the dictionary and vector priorities dynamic based on the priority of their input types
union_ordering(::Type{<:AbstractDict}) = 7.0
union_ordering(::Type{<:AbstractVector}) = 8.0
const MAX_BUILTIN_UNION_PRIORITY = union_ordering(Vector)

"Organizes the types in a union by their `union_ordering()`."
@inline union_parse_order(U::Union) = sort(collect(union_types(U)), by=union_ordering)
union_parse_order(T::Type) = (T, )
#NOTE: originally I used TupleTools.sort(), but it breaks
#    for unions containing DataTypes and UnionAll's, because
#    those are two different types of things, so the tuple isn't homogenous.

"
Tries parsing a type from a union, returning `Nothing` if it failed
    or `Some{T}` if it succeeded.
By default, simply performs the conversion and catches any exception that happens.
"
function union_try_parse(converter, toml_data, union_type)::Union{Nothing, Some}
    try
        return Some(converter(toml_data, union_type))
    catch e
        return nothing
    end
end
# Optimize the try-parse for simple types.
union_try_parse(_, x::T, ::Type{T2}) where {T2, T<:T2} = Some(x)
union_try_parse(_, f::Real, N::Type{<:Real}) =
    (hasmethod(N, tuple(typeof(f))) && (f >= typemin(N)) && (f <= typemax(N))) ?
        Some(N(f)) :
        nothing
union_try_parse(_, f::Real, I::Type{<:Integer}) =
    (isinteger(f) && (f >= typemin(I)) && (f <= typemax(I)) && hasmethod(I, tuple(typeof(f)))) ?
        Some(I(f)) :
        nothing
union_try_parse(_, s::AbstractString, N::Type{<:Real}) = let parsed = tryparse(N, s)
    return isnothing(parsed) ? nothing : Some(parsed)
end
union_try_parse(_, s::AbstractString, I::Type{<:Integer}) = let parsed = tryparse(Float64, s)
    return (!isnothing(parsed) && isinteger(parsed)) ? Some(I(parsed)) : nothing
end
union_try_parse(_, x, ::Type{Symbol}) =
    hasmethod(Symbol, tuple(typeof(x))) ?
        Some(Symbol(x)) :
        nothing
union_try_parse(c, x::AbstractString, ::Type{Nothing}) =
    (x == c.null_string) ?
        Some(nothing) :
        nothing


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

const moon_d_lng = SVector(0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0, 4, 0, 4, 2,
                           2, 1, 1, 2, 2, 4, 2, 0, 2, 2, 1, 2, 0, 0, 2, 2, 2, 4, 0, 3,
                           2, 4, 0, 2, 2, 2, 4, 0, 4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2, 2)
const moon_M_lng = SVector(0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0, 0, 0, 0,
                           1, 1, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, -2, 1, 2, -2, 0, 0,
                           -1, 0, 0, 1, -1, 2, 2, 1, -1, 0, 0, -1, 0, 1, 0, 1, 0, 0,
                           -1, 2, 1, 0, 0)
const moon_Mprime_lng = SVector(1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1, 1, -1,
                                3, -2, -1, 0, -1, 0, 1, 2, 0, -3, -2, -1, -2, 1, 0, 2,
                                0, -1, 1, 0, -1, 2, -1, 1, -2, -1, -1, -2, 0, 1, 4, 0,
                         -2, 0, 2, 1, -2, -3, 2, 1, -1, 3, -1)
const moon_F_lng = SVector(0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2, 0, 0, 0,
                           0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, -2, 2, 0,
                           2, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, 0, 0,
                           0, 0, -2)
const moon_sin_lng = [6288774, 1274027, 658314, 213618, -185116, -114332, 58793,
                      57066, 53322, 45758, -40923, -34720, -30383, 15327,
                      -12528, 10980, 10675, 10034, 8548, -7888, -6766,
                      -5163, 4987, 4036, 3994, 3861, 3665, -2689, -2602, 2390,
                      -2348, 2236, -2120, -2069, 2048, -1773, -1595, 1215,
                      -1110, -892, -810, 759, -713, -700, 691, 596, 549, 537,
                      520, -487, -399, -381, 351, -340, 330, 327, -323, 299, 294,
                      0]
const moon_cos_lng = [-20905355, -3699111, -2955968, -569925, 48888, -3149,
                      246158, -152138, -170733, -204586, -129620, 108743,
                      104755, 10321, 0, 79661, -34782, -23210, -21636, 24208,
                      30824, -8379, -16675, -12831, -10445, -11650, 14403,
                      -7003, 0, 10056, 6322, -9884, 5751, 0, -4950, 4130, 0,
                      -3958, 0, 3258, 2616, -1897, -2117, 2354, 0, 0, -1423,
                      -1117, -1571, -1739, 0, -4421, 0, 0, 0, 0, 1165, 0, 0,
                      8752]
const moon_d_lat = SVector(0, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 0, 4, 0,
                           0, 0, 1, 0, 0, 0, 1, 0, 4, 4, 0, 4, 2, 2, 2, 2, 0, 2, 2, 2,
                           2, 4, 2, 2, 0, 2, 1, 1, 0, 2, 1, 2, 0, 4, 4, 1, 4, 1, 4, 2)
const moon_M_lat = SVector(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, -1, -1, 1, 0,
                           1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0,
                           1, 1, 0, -1, -2, 0, 1, 1, 1, 1, 1, 0, -1, 1, 0, -1, 0, 0, 0,
                           -1, -2)
const moon_Mprime_lat = SVector(0, 1, 1, 0, -1, -1, 0, 2, 1, 2, 0, -2, 1, 0, -1, 0, -1,
                                -1, -1, 0, 0, -1, 0, 1, 1, 0, 0, 3, 0, -1, 1,  -2, 0,
                                2, 1, -2, 3, 2, -3, -1, 0, 0, 1, 0, 1, 1, 0, 0, -2, -1,
                                1, -2, 2, -2, -1, 1, 1, -1, 0, 0)
const moon_F_lat = SVector(1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1,
                           -1, -1, 1, 3, 1, 1, 1, -1, -1, -1, 1, -1, 1, -3, 1, -3, -1,
                           -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 3, -1, -1, 1, -1, -1, 1,
                           -1, 1, -1, -1, -1, -1, -1, -1, 1)
const moon_sin_lat = [5128122, 280602, 277693, 173237, 55413, 46271, 32573,
                      17198, 9266, 8822, 8216, 4324, 4200, -3359, 2463, 2211,
                      2065, -1870, 1828, -1794, -1749, -1565, -1491, -1475,
                      -1410, -1344, -1335, 1107, 1021, 833, 777, 671, 607, 596,
                      491, -451, 439, 422, 421, -366, -351, 331, 315, 302, -283,
                      -229, 223, 223, -220, -220, -185, 181, -177, 176, 166,
                      -164, 132, -119, 115, 107]

function _moonpos(jd::T, radians::Bool) where {T<:AbstractFloat}
    # Number of Julian centuries since 2000-01-01T12:00:00
    t = (jd - J2000) / JULIANCENTURY
    # Mean longitude of the moon referred to mean equinox of the date
    Lprimed = mod(@evalpoly(t, 218.3164477, 481267.88123421,
                            -0.0015786, inv(538841), -inv(6.5194e7)), 360)
    Lprime = deg2rad(Lprimed)
    # Mean elongation of the Moon
    d = deg2rad(mod(@evalpoly(t, 297.8501921, 445267.1114034, -0.0018819,
                              inv(545868), -inv(1.13065e8)), 360))
    # Sun's mean anomaly
    M = deg2rad(mod(@evalpoly(t, 357.5291092, 35999.0502909, -0.0001536,
                              inv(2.449e7)), 360))
    # Moon's mean anomaly
    Mprime = deg2rad(mod(@evalpoly(t, 134.9633964, 477198.8675055,
                                   0.0087414, inv(6.9699e4),
                                   -inv(1.4712e7)), 360))
    # Moon's argument of latitude
    F = deg2rad(mod(@evalpoly(t, 93.2720950, 483202.0175233, -0.0036539,
                              -inv(3.526e7), inv(8.6331e8)), 360))
    # Eccentricity of Earth's orbit around the Sun
    E = @evalpoly t 1 -0.002516 -7.4e-6
    E2 = E*E
    # Additional arguments
    A1 = deg2rad(119.75 + 131.849*t)
    A2 = deg2rad(53.09 + 479264.290*t)
    A3 = deg2rad(313.45 + 481266.484*t)
    suml_add = 3958*sin(A1) + 1962*sin(Lprime - F) + 318*sin(A2)
    sumb_add = -2235*sin(Lprime) + 382*sin(A3) + 175*sin(A1 - F) +
        175*sin(A1 + F) + 127*sin(Lprime - Mprime) -
        115*sin(Lprime + Mprime)
    # Sum the period terms
    sinlng = Vector{T}(moon_sin_lng)
    coslng = Vector{T}(moon_cos_lng)
    sinlat = Vector{T}(moon_sin_lat)
    for i in eachindex(moon_M_lng)
        if abs(moon_M_lng[i]) == 1
            sinlng[i] = E*sinlng[i]
            coslng[i] = E*coslng[i]
        elseif abs(moon_M_lng[i]) == 2
            sinlng[i] = E2*sinlng[i]
            coslng[i] = E2*coslng[i]
        end
        if abs(moon_M_lat[i]) == 1
            sinlat[i] = E*sinlat[i]
        elseif abs(moon_M_lat[i]) == 2
            sinlat[i] = E2*sinlat[i]
        end
    end
    arg = moon_d_lng * d + moon_M_lng * M + moon_Mprime_lng * Mprime + moon_F_lng * F
    geolong = Lprimed + (dot(sinlng, sin.(arg)) + suml_add)/1000000
    dis = 385000.56 + dot(coslng, cos.(arg))/1000
    arg = moon_d_lat * d + moon_M_lat * M + moon_Mprime_lat * Mprime + moon_F_lat * F
    geolat = (dot(sinlat, sin.(arg)) + sumb_add) / 1000000
    nlong, elong = nutate(jd)
    geolong = mod(geolong + nlong/3600, 360)
    λ = deg2rad(geolong)
    β = deg2rad(geolat)
    # Find mean obliquity and convert λ, β to right ascension and declination.
    ɛ = ten(23, 26) + @evalpoly(t/100, 21.448, -4680.93, -1.55, 1999.25, -51.38,
                                -249.67, -39.05, 7.12, 27.87, 5.79, 2.45)/3600
    ɛ = deg2rad(ɛ + elong/3600)
    sin_β, cos_β = sincos(β)
    sin_ɛ, cos_ɛ = sincos(ɛ)
    sin_λ, cos_λ = sincos(λ)
    ra = mod2pi(atan(sin_λ * cos_ɛ - tan(β) * sin_ɛ, cos_λ))
    dec = asin(sin_β * cos_ɛ + cos_β * sin_ɛ * sin_λ)
    if radians
        return ra, dec, dis, λ, β
    else
        return rad2deg(ra), rad2deg(dec), dis, geolong, geolat
    end
end

"""
    moonpos(jd[, radians=true]) -> ra, dec, dis, geolong, geolat

### Purpose ###

Compute the right ascension and declination of the Moon at specified Julian
date.

### Arguments ###

* `jd`: the Julian ephemeris date.  It can be either a scalar or an array
* `radians` (optional boolean keyword): if set to `true`, then all output
  angular quantities are given in radians rather than degrees.  The default is
  `false`

### Output ###

The 5-tuple `(ra, dec, dis, geolong, geolat)`:

* `ra`: apparent right ascension of the Moon in degrees, referred to the true
  equator of the specified date(s)
* `dec`: the declination of the Moon in degrees
* `dis`: the distance between the centre of the Earth and the centre of the Moon
  in kilometers
* `geolong`: apparent longitude of the moon in degrees, referred to the ecliptic
  of the specified date(s)
* `geolat`: apparent longitude of the moon in degrees, referred to the ecliptic
  of the specified date(s)

If `jd` is an array, then all output quantities are arrays of the same length as
`jd`.

### Method ###

Derived from the Chapront ELP2000/82 Lunar Theory (Chapront-Touze' and Chapront,
1983, 124, 50), as described by Jean Meeus in Chapter 47 of ``Astronomical
Algorithms'' (Willmann-Bell, Richmond), 2nd edition, 1998.  Meeus quotes an
approximate accuracy of 10" in longitude and 4" in latitude, but he does not
give the time range for this accuracy.

Comparison of the IDL procedure with the example in ``Astronomical Algorithms''
reveals a very small discrepancy (~1 km) in the distance computation, but no
difference in the position calculation.

### Example ###

(1) Find the position of the moon on April 12, 1992

```jldoctest
julia> using AstroLib

julia> jd = jdcnv(1992, 4, 12);

julia> adstring(moonpos(jd)[1:2],precision=1)
" 08 58 45.23  +13 46 06.1"
```

This is within 1" from the position given in the Astronomical Almanac.

(2) Plot the Earth-moon distance during 2016 with sampling of 6 hours.  Use
[PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting

```julia
using PyPlot
points = DateTime(2016):Dates.Hour(6):DateTime(2017);
plot(points, moonpos(jdcnv.(points))[3])
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
moonpos(jd::Real; radians::Bool=false) = _moonpos(float(jd), radians)

function moonpos(jd::AbstractArray{J}; radians::Bool=false) where {J<:Real}
    typej = float(J)
    ra = similar(jd, typej)
    dec = similar(jd, typej)
    dis = similar(jd, typej)
    long = similar(jd, typej)
    lat = similar(jd, typej)
    for i in eachindex(jd)
        ra[i], dec[i], dis[i], long[i], lat[i] = moonpos(jd[i], radians=radians)
    end
    return ra, dec, dis, long, lat
end


text:
using Knet,ProgressMeter
include(Knet.dir("data/imagenet.jl"))
include(Knet.dir("examples/resnet/resnet.jl"))

function loadparams(atype)
    mname       = "imagenet-resnet-101-dag"
    model       = matconvnet(mname)
    avgimg      = Array{Float32}(model["meta"]["normalization"]["averageImage"])
    w, ms       = ResNet.get_params(model["params"], atype)
    deleteat!(w,283:length(w))
    deleteat!(ms,189:length(ms))
    return w,ms,avgimg
end

function resnet101(w,x,ms; mode=1)
    # layer 1
    conv1 = ResNet.reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = ResNet.pool(conv1; window=3, stride=2)
    # layer 2,3
    r2 = ResNet.reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = ResNet.reslayerx5(w[34:72], r2, ms; mode=mode)
    return ResNet.reslayerx5(w[73:282], r3, ms; mode=mode)
end

function imgdata(img,averageImage)
    a0 = RGB.(load(img))
    new_size = ntuple(i->div(size(a0,i)*224,minimum(size(a0))),2)
    a1 = Images.imresize(a0, new_size)
    i1 = div(size(a1,1)-224,2)
    j1 = div(size(a1,2)-224,2)
    b1 = a1[i1+1:i1+224,j1+1:j1+224]
    c1 = permutedims(channelview(b1),(3,2,1))
    d1 = convert(Array{Float32},c1)
    e1 = reshape(d1, (224,224,3,1))
    f1 = (255 * e1 .- averageImage)
    g1 = permutedims(f1, [2,1,3,4])
end

function extract(dhome,set,params;atype=Array{Float32})
    w,ms,avgimg = params
    sethome = joinpath(dhome,"images",set)*"/"
    binfile = "data/"*set*".bin"
    files   = readdir(sethome)
    p = Progress(length(files))
    for file in files
        if file[end-2:end] == "png"
            x = convert(atype,imgdata(sethome*file,avgimg))
            y = convert(Array{Float32},resnet101(w,x,ms; mode=1))
            open(binfile,"a+") do f
                write(f,y)
            end
        end
        next!(p)
    end
end

function extract_features(dhome)
    atype = gpu()<0 ? Array{Float32} : KnetArray{Float32}
    params = loadparams(atype)
    println("Extracting training features")
    extract(dhome,"train",params;atype=atype)
    println("Extracting validation features")
    extract(dhome,"val",params;atype=atype)
end

max_stars_repo_name: feyzaakyurek/MAC-Network

max_issues_repo_name: feyzaakyurek/MAC-Network

max_forks_repo_name: feyzaakyurek/MAC-Network


text:
######################################################################
# Unit Tests for Pia.jl library
# -----
# Licensed under Apache License 2.0

include("./Test.Futures.jl")
include("./Test.Lockables.jl")
include("./Test.Logging.jl")
include("./Test.Signals.jl")
include("./Test.Voidables.jl")
include("./Test.Pia.jl")


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

doc"""An `IntervalBox` is an $N$-dimensional rectangular box, given
by a Cartesian product of $N$ `Interval`s.
"""
struct IntervalBox{N,T} <: StaticVector{N, Interval{T}}
    data::NTuple{N,Interval{T}}
end

# StaticArrays.Size{N,T}(::Type{IntervalBox{N,T}}) = StaticArrays.Size(N) # @pure not needed, I think...
Base.getindex(a::IntervalBox, i::Int) = a.data[i]

IntervalBox(x::Interval) = IntervalBox( (x,) )  # single interval treated as tuple with one element

## arithmetic operations
# Note that standard arithmetic operations are implemented automatically by FixedSizeArrays.jl

mid(X::IntervalBox) = mid.(X)

## set operations

# TODO: Update to use generator
⊆(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} =
    all(X .⊆ Y)

∩(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} =
    IntervalBox(X .∩ Y)
∪(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} =
    IntervalBox(X .∪ Y)

#=
On Julia 0.6 can now write
∩{N,T}(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) = IntervalBox(NTuple{N, Interval{Float64}}( (X[i] ∩ Y[i]) for i in 1:N))
=#

isempty(X::IntervalBox) = any(isempty, X)

diam(X::IntervalBox) = maximum(diam.(X))

emptyinterval(X::IntervalBox{N,T}) where {N,T} = IntervalBox(emptyinterval.(X))

import Base.×
×(a::Interval...) = IntervalBox(a...)
×(a::Interval, b::IntervalBox) = IntervalBox(a, b...)
×(a::IntervalBox, b::Interval) = IntervalBox(a..., b)
×(a::IntervalBox, b::IntervalBox) = IntervalBox(a..., b...)

IntervalBox(x::Interval, ::Type{Val{n}}) where {n} = IntervalBox(SVector(ntuple(i->x, Val{n})))

IntervalBox(x::Interval, n::Int) = IntervalBox(x, Val{n})


text:
using Knet, CUDA

pd_knet(y::KnetArray, x::KnetArray, perm) = permutedims!(y,x,perm)
pd_cpux(y::KnetArray, x::KnetArray, perm) = copyto!(y, permutedims(Array(x),perm))
pd_cuxx(y::KnetArray, x::KnetArray, perm) = (permutedims!(cu(y),cu(x),perm); y)
pd_kern(y::KnetArray, x::KnetArray, perm) = nothing # call specific kernel
pd_tran(y::KnetArray, x::KnetArray, perm) = Knet._transpose!(y,x)

function CUDA.cu(x::KnetArray{T}) where {T}
    p = CUDA.CuPtr{T}(UInt(x.ptr.ptr))
    Base.unsafe_wrap(CuArray{T}, p, size(x); own=false)
end

# ## best not to use CUDA.jl memory manager
# function Knet.ka(x::CuArray{T,N}) where {T,N}
#     p = Base.bitcast(Knet.Cptr, x.buf.ptr)
#     k = Knet.KnetPtr(p, sizeof(x), gpu(), x)
#     # finalizer(identity, k) # hacky way to avoid gc? gives error in running finalizer
#     KnetArray{T,N}(k, size(x))
# end

function pd_time(f, x, p, a=nothing, n=max(10,div(1000000,length(x))))
    println(f,size(x),p,n)
    y = similar(x, size(x)[[p...]])
    f(y,x,p)
    if a != nothing; @assert y == a; end
    @time begin
        for i in 1:n
            f(y,x,p)
        end
        Knet.cudaDeviceSynchronize()
    end
    @time begin
        for i in 1:n
            f(y,x,p)
        end
        Knet.cudaDeviceSynchronize()
    end
    return y
end

@info "Simple transpose"
p = (2,1)
x = KnetArray(rand(Float32,100,100))
a = pd_time(pd_cpux, x, p)
pd_time(pd_tran, x, p, a)
pd_time(pd_knet, x, p, a) # confirm this uses transpose
pd_time(pd_cuxx, x, p, a)
#pd_time(pd_kern, x, p)
x = y = a = nothing; Knet.gc(); #Knet.gpuinfo()

@info "ILKER-DRAW"
p = (2,1,3)
x = KnetArray(rand(Float32,28,5,100))
a = pd_time(pd_cpux, x, p)
pd_time(pd_knet, x, p, a)
pd_time(pd_cuxx, x, p, a)
#pd_time(pd_kern, x, p, a)
x = y = a = nothing; Knet.gc(); #Knet.gpuinfo()

@info "ILKER-ATTN"
p = (2,1,3)
x = KnetArray(rand(Float32,196,512,32))
a = pd_time(pd_cpux, x, p)
pd_time(pd_knet, x, p, a)
pd_time(pd_cuxx, x, p, a)
#pd_time(pd_kern, x, p, a)
x = y = a = nothing; Knet.gc(); #Knet.gpuinfo()

@info "OSMAN-BERT-1"
p = (2,1,3)
x = KnetArray(rand(Float32,768,512,16))
a = pd_time(pd_cpux, x, p)
pd_time(pd_knet, x, p, a)
pd_time(pd_cuxx, x, p, a)
#pd_time(pd_kern, x, p, a)
x = a = nothing; Knet.gc(); #Knet.gpuinfo()

@info "OSMAN-BERT-2"
x = KnetArray(rand(Float32,64,512,12,16))
p = (1,3,2,4)
a = pd_time(pd_cpux, x, p)
pd_time(pd_knet, x, p, a)
pd_time(pd_cuxx, x, p, a)
#pd_time(pd_kern, x, p, a)
x = a = nothing; Knet.gc(); #Knet.gpuinfo()

nothing


text:
using UUIDs

include("rewrite.jl")

function docs_url(deployment_url, name, uuid, version)
    isempty(deployment_url) && return nothing
    return join([deployment_url, DocumentationGenerator.get_docs_dir(name, uuid), version], '/')
end

get_docs_dir(name, uuid) = get_docs_dir(name, UUID(uuid))
get_docs_dir(name, uuid::UUID) = joinpath(name, Base.package_slug(uuid, 5))

function find_free_x_servernum()
    i = 99
    while isfile("/tmp/.X$(i)-lock")
        i += 1
    end

    return i
end


text:
export compute_BEM_entry

# In this file we implement the `compute_BEM_entry` function.
# This function determines whether the corresponding integral is regular,
# singular or nearly singular. It differentiates between collocation and
# Galerkin discretization to do so.
# In the end, the function invokes `projectionintegral` (for collocation)
# or `doubleprojectionintegral` (for Galerkin).

"Compute the BEM matrix entry A[i,j]."
compute_BEM_entry(op::DenseBEMOperator, i::Int, j::Int, quad = op.quad) =
    compute_BEM_entry(i, (src(op), j), samplingoperator(op), integraloperator(op), quad)

# Catch collocation/galerkin based on the sampling operator
compute_BEM_entry(i, ϕ_j, sampling::GridSampling, intop, quad) =
    collocation_BEM_entry(grid(sampling), i, ϕ_j, intop, quad)

compute_BEM_entry(i, ϕ_j, sampling::ProjectionSampling, intop, quad) =
    galerkin_BEM_entry((dictionary(sampling), i), ϕ_j, measure(sampling), intop, quad)

##############
# Collocation
##############

collocation_BEM_entry(grid::AbstractGrid, i, ϕ_j, intop, quad) =
    collocation_BEM_entry(grid[i], ϕ_j, intop, quad)

"Does the collocation entry correspond to a (nearly) singular integral, up to the given distance threshold `δ`?"
issingular_coll(sing, x, ϕ_j) =
    issingular_coll(sing, x, ϕ_j, DomainSets.default_tolerance(support(ϕ_j...)))

issingular_coll(::NoSingularity, x, ϕ_j, δ) = false

issingular_coll(sing, x, ϕ_j, δ) =
    _issingular_coll(x, support(ϕ_j...), δ)
_issingular_coll(x, d::AbstractInterval, δ) = approx_in(x, d, δ)
_issingular_coll(x, d::MappedDomain, δ) = _issingular_coll(inverse_map(d, x), superdomain(d), δ)
# Line below catches a corner case
_issingular_coll(x, d::PeriodicInterval, δ) =
    approx_in(x, d, δ) || approx_in(x+period(d), d, δ) || approx_in(x-period(d), d, δ)

"Which threshold decides whether an integral is nearly singular?"
nearsingular_threshold(ϕ_j) = nearsingular_threshold(support(ϕ_j...))
# By default we choose twice the length of a single basis function
nearsingular_threshold(d::AbstractInterval) = 2DomainSets.width(d)
nearsingular_threshold(d::PeriodicInterval) = 2sum(map(DomainSets.width, components(d)))
# we don't know the radius of a general domain...
nearsingular_threshold(d::Domain) = 0.01

function collocation_BEM_entry(x, ϕ_j, intop, quad, sing = singularity(intop))
    if issingular_coll(sing, x, ϕ_j)
        singular_collocation_entry(singular(quad), x, ϕ_j, intop, sing)
    elseif issingular_coll(sing, x, ϕ_j, nearsingular_threshold(ϕ_j))
        nearly_singular_collocation_entry(nearlysingular(quad), x, ϕ_j, intop, sing)
    else
        regular_collocation_entry(regular(quad), x, ϕ_j, intop)
    end
end

collocation_singularity(x, sing) = sing
collocation_singularity(x, sing::LogSingularDiagonal) = LogSingPoint(x)

singular_collocation_entry(quad, x, ϕ_j, intop, sing) =
    projectionintegral(singular(quad), y -> eval_kernel(intop, x, y),
            ϕ_j..., measure(intop), collocation_singularity(x, sing))

nearly_singular_collocation_entry(quad, x, ϕ_j, intop, sing) =
    projectionintegral(nearlysingular(quad), y -> eval_kernel(intop, x, y),
            ϕ_j..., measure(intop), collocation_singularity(x, sing))

regular_collocation_entry(quad, x, ϕ_j, intop) =
    projectionintegral(regular(quad), y -> eval_kernel(intop, x, y), ϕ_j..., measure(intop))

##############
# Galerkin
##############

"Does the Galerkin entry correspond to a (nearly) singular integral, up to the given distance threshold `δ`?"
issingular_galerkin(sing, ϕ_i, ϕ_j) =
    issingular_galerkin(sing, ϕ_i, ϕ_j, DomainSets.default_tolerance(support(ϕ_j...)))

issingular_galerkin(::NoSingularity, ϕ_i, ϕ_j, δ) = false
issingular_galerkin(sing, ϕ_i, ϕ_j, δ) =
    _issingular_galerkin(sing, ϕ_i, ϕ_j, δ, support(ϕ_i...), support(ϕ_j...))
_issingular_galerkin(sing, ϕ_i, ϕ_j, δ, domain1, domain2) =
    distance(domain1, domain2) < δ

function distance(domain1::AbstractInterval, domain2::AbstractInterval)
    a1 = leftendpoint(domain1)
    b1 = rightendpoint(domain1)
    a2 = leftendpoint(domain2)
    b2 = rightendpoint(domain2)
    if a2 > b1
        a2-b1
    elseif a1 > b2
        a1-b2
    else
        zero(a1)
    end
end

function distance(domain1::PeriodicInterval, domain2::PeriodicInterval)
    @assert domain1.periodicdomain == domain2.periodicdomain
    period = width(domain1.periodicdomain)
    sub1 = domain1.subdomain
    sub2 = domain2.subdomain
    min(distance(sub1, sub2), distance(sub1 .+ period,sub2), distance(sub1 .- period,sub2))
end

distance(domain1::UnitInterval, domain2::PeriodicInterval) = zero(eltype(domain1))

function galerkin_BEM_entry(ϕ_i, ϕ_j, μ_projection, intop, quad, sing = singularity(intop))
    if issingular_galerkin(sing, ϕ_i, ϕ_j)
        singular_galerkin_BEM_entry(singular(quad), sing, ϕ_i, μ_projection, ϕ_j, intop)
    elseif issingular_galerkin(sing, ϕ_i, ϕ_j, nearsingular_threshold(ϕ_j))
        nearly_singular_galerkin_BEM_entry(nearlysingular(quad), sing, ϕ_i, μ_projection, ϕ_j, intop)
    else
        regular_galerkin_BEM_entry(regular(quad), ϕ_i, μ_projection, ϕ_j, intop)
    end
end

galerkin_singularity(sing) = sing
galerkin_singularity(sing::LogSingularDiagonal) = SingularDiagonal()

singular_galerkin_BEM_entry(quad, sing, ϕ_i, μ_projection, ϕ_j, intop) =
    doubleprojection(singular(quad), (x,y) -> eval_kernel(intop, x, y),
            ϕ_j..., measure(intop), ϕ_i..., μ_projection, galerkin_singularity(sing))

nearly_singular_galerkin_BEM_entry(quad, sing, ϕ_i, μ_projection, ϕ_j, intop) =
    doubleprojection(nearlysingular(quad), (x,y) -> eval_kernel(intop, x, y),
            ϕ_j..., measure(intop), ϕ_i..., μ_projection, galerkin_singularity(sing))

regular_galerkin_BEM_entry(quad, ϕ_i, μ_projection, ϕ_j, intop) =
    doubleprojection(regular(quad), (x,y) -> eval_kernel(intop, x, y), ϕ_j...,
            measure(intop), ϕ_i..., μ_projection)

#######################
# Ranges and subblocks
#######################

# There is overhead in the computation of a single entry that can be
# saved if we are computing a range of entries.

compute_BEM_range!(result, A::DenseBEMOperator, I, J, quad = A.quad) =
    compute_BEM_range!(result, I, J, quad, A.src, A.sampling, A.intop)

compute_BEM_range!(result, I, J, quad, dict, sampling::GridSampling, intop) =
    collocation_BEM_range!(result, grid(sampling), I, quad, dict, J, intop, measure(intop))

collocation_BEM_range!(result, grid::AbstractGrid, i::Int, quad, dict, J, intop, measure) =
    collocation_BEM_range!(result, grid[i], i, quad, dict, J, intop, measure)

function collocation_BEM_range!(result, x, i, quad, dict, J, intop, measure)
    sing = singularity(intop)
    for j in J
        result[i,j] = collocation_BEM_entry(x, quad, dict, j, intop, measure, sing)
    end
    result
end

function collocation_BEM_range!(result, grid::AbstractGrid, I, quad, dict, J, intop, measure)
    for i in I
        collocation_BEM_range!(result, grid, i, quad, dict, J, intop, measure)
    end
    result
end