text:
# Julia calling python for vizualization in MIT colletions viewer.
# We should converge with Director, but need the collections renderer code
# in Director before that will happen
#
# This code uses Sudeep's Python drawing utils and transformation code

using Caesar, RoME, KernelDensityEstimate
using IncrementalInference, CloudGraphs, Neo4j
using TransformUtils

using PyCall

# Uncomment out for command line operation
cloudGraph, addrdict = standardcloudgraphsetup(drawdepth=true)
session = addrdict["session"]
DRAWDEPTH = addrdict["draw depth"]=="y" || addrdict["draw depth"]=="yes"

mongoaddress = addrdict["mongo addr"]
collection = "bindata"

# # interactive operation
# session = "SESSTURT21"
# Nparticles = 100
# include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))
# DRAWDEPTH = true

# include("VisualizationUtilities.jl")  # @pyimport getimages as gi

@pyimport pybot.geometry.rigid_transform as bgrt
@pyimport pybot.geometry.quaternion as bgq

# usage example
#	bgq.Quaternion[:from_rpy](0,0,0)

# this resets the Collections Viewer visualization
@pyimport pybot.externals.lcm.draw_utils as bedu
# @pyimport pybot.vision.camera_utils as cu

@pyimport numpy as np
@pyimport cv2 as opencv

# draw_utils.publish_pose_list('apriltag', self.tagwposes, frame_id='origin', texts=ids)

# CAMK [[ 570.34222412    0.          319.5       ]
#  [   0.          570.34222412  239.5       ]
#  [   0.            0.            1.        ]]

CAMK = [[ 570.34222412; 0.0; 319.5]';
 [   0.0; 570.34222412; 239.5]';
 [   0.0; 0.0; 1.0]'];

# dcam = cu.DepthCamera(K=CAMK)
# -np.pi/2, 0, -np.pi/2
dcamjl = DepthCamera(CAMK)
buildmesh!(dcamjl)

temp = bgrt.RigidTransform[:from_rpyxyz](-pi/2, 0, -pi/2, 0, 0, 0.6, axes="sxyz")
# @show temp[:to_matrix]()

# bTc = SE3([0;0;0.6],AngleAxis(-pi/2,[1;0;0])*AngleAxis(-pi/2,[0;0;1]))
bTc = SE3(temp[:to_matrix]())

poseswithdepth = Dict()
poseswithdepth["x1"] = 0 # skip this pose -- there is no big data before ICRA
# poseswithdepth["x52"] = 0 # skip this pose -- there is no big data before ICRA
# poseswithdepth["x53"] = 0 # skip this pose -- there is no big data before ICRA
# poseswithdepth["x54"] = 0 # skip this pose -- there is no big data before ICRA
# poseswithdepth["x55"] = 0 # skip this pose -- there is no big data before ICRA

while true

  # this is being replaced by cloudGraph, added here for development period
  fg = Caesar.initfg(sessionname=session, cloudgraph=cloudGraph)

  IDs = getPoseExVertexNeoIDs(fg.cg.neo4j.connection, session=session, reqbackendset=false);

  println("get local copy of graph")
  if fullLocalGraphCopy!(fg, reqbackendset=false)
  	xx,ll = ls(fg)
  	LD = Array{Array{Float64,1},1}()
  	C = Vector{AbstractString}()
  	for x in xx
  		val = getVal(fg, x)
  		len = size(val,2)
  		[push!(C,"g") for i in 1:len];
  		[push!(LD, vec([val[1:2,i];0])) for i in 1:len];
  	end
  	bedu.publish_cloud("posePts", LD,c=C, frame_id="origin")
  	# draw pose MAP triads
  	# Xpp,Ypp, Thpp, LBLS = get2DPoseMax(fg)
  	poses = []
    cams = []
  	# for i in 1:length(Xpp)
    j=0
  	for x in xx
      j+=1
  		# posei = bgrt.RigidTransform[:from_angle_axis](Thpp[i],[0;0;1],[Xpp[i];Ypp[i];0])
  		vert = getVert(fg,x, api=localapi)
  		X = Float64[]
  		if haskey(vert.attributes, "MAP_est")
  			X = vert.attributes["MAP_est"]
  		else
  			X = getKDEfit(getKDE(vert))
  		end
  		posei = bgrt.RigidTransform[:from_angle_axis](X[3],[0;0;1],[X[1];X[2];0])
  		push!(poses, posei)
      wTb = SE3([X[1];X[2];0], TransformUtils.AngleAxis(X[3],[0;0;1]))
      wTc = wTb*bTc
      camiRT = bgrt.RigidTransform[:from_matrix](matrix(wTc))
      cami = bgrt.Pose[:from_rigid_transform](j, camiRT)

      push!(cams, cami)
  	end
  	bedu.publish_pose_list("MAPposes",poses, frame_id="origin",texts=xx)
    bedu.publish_pose_list("MAPcams",cams, frame_id="origin",reset=false)

    j=0
    for x in xx
      j+=1
      mongk, cvid = getmongokeys(fg, x, IDs)
  		if DRAWDEPTH && haskey(mongk, "depthframe_image") # !haskey(poseswithdepth, x)
        poseswithdepth[x]=1
        mongo_keydepth = bson.ObjectId(mongk["depthframe_image"])
        img, ims = gi.fastdepthimg(db[collection], mongo_keydepth)
        # cvid = -1
        # for id in IDs
        #   if Symbol(fg.g.vertices[id[1]].label) == x
        #     cvid = id[2]
        #     break
        #   end
        # end
        # # fetch depth cloud from mongo for pose
        # cv = CloudGraphs.get_vertex(fg.cg, cvid)
        # bd = CloudGraphs.read_BigData!(cloudGraph, cv)
        # depthcloudpng = nothing
        # for i in bd.dataElements
        #   if i.description == "depthframe-image"
        #     depthcloudpng = i.data
        #     break
        #   end
        # end
        # #interpret data
        # fid = open("tempdepth.png","w")
        # write(fid, depthcloudpng)
        # close(fid)
        # img = opencv.imread("tempdepth.png",2)
        # imgc = map(Float64,img)/1000.0
        #   opencv.imshow("yes", img)
        # calibrate the image # color conversion of points, so we can get pretty pictures...
        # X = dcam[:reconstruct](img)
        X = reconstruct(dcamjl, Array{Float64,2}(img))
        r,c,h = size(X)
        Xd = X[1:3:r,1:3:c,:]
        mask = Xd[:,:,:] .> 4.5
        Xd[mask] = Inf
        # get color information
        rgb = nothing
        seg = nothing
        if haskey(mongk, "keyframe_rgb")
          mongo_key = bson.ObjectId(mongk["keyframe_rgb"])
          rgb, ims = gi.fastrgbimg(db[collection], mongo_key)
          # @show size(rgb)
        end
        if haskey(mongk, "keyframe_segnet")
          mongo_key = bson.ObjectId(mongk["keyframe_segnet"])
          seg, ims = gi.fastrgbimg(db[collection], mongo_key)
        end
        # for i in bd.dataElements
        #   if i.description == "keyframe-image"
        #     rgbpng = i.data
        #   end
        #   if i.description == "segnet-image"
        #     segpng = i.data
        #   end
        # end
        #interpret data
        if rgb != nothing
          # fid = open("temprgb.png","w")
          # write(fid, rgbpng)
          # close(fid)
          # rgb = opencv.imread("temprgb.png")
          # rgb = bgr[:,:,[3;2;1]]
          # publish to viewer via Sudeeps python code
          # @show size(rgb), typeof(rgb)
          rgbss = rgb[1:3:r,1:3:c,:]
          bedu.publish_cloud("depth", Xd, c=rgbss, frame_id="MAPcams",element_id=j, flip_rb=true, reset=false)
        end
        if seg != nothing
          # fid = open("temprgb.png","w")
          # write(fid, segpng)
          # close(fid)
          # seg = opencv.imread("temprgb.png")
          # @show size(seg), typeof(seg)
          segss = seg[1:3:r,1:3:c,:]
          bedu.publish_cloud("segnet", Xd, c=segss, frame_id="MAPcams",element_id=j, flip_rb=true, reset=false)
        end
      end
    end

  	# now do landmarks
  	LD = Array{Array{Float64,1},1}()
  	C = Vector{AbstractString}()
    landms = []
  	for l in ll
      vert = getVert(fg,l, api=localapi)

  		val = getVal(vert)
  		len = size(val,2)
  		[push!(C,"c") for i in 1:len];
  		[push!(LD, vec([val[1:2,i];0])) for i in 1:len];

  		X = Float64[]
  		if haskey(vert.attributes, "MAP_est")
  			X = vert.attributes["MAP_est"]
  		else
  			X = getKDEMax(getKDE(vert))
  		end

      lmi = bgrt.RigidTransform[:from_angle_axis](0,[0;0;1],[X[1];X[2];0])
  		push!(landms, lmi)

  	end
  	bedu.publish_cloud("landmarkPts", LD,c=C, frame_id="origin")
    bedu.publish_pose_list("MAPlandms",landms, frame_id="origin",texts=ll)
  else
    sleep(0.2)
  end
  sleep(0.05)

end

#


text:
is_ch_type(::Val{:Tuple})  = true
can_be_nullable(::Val{Tuple}) = false

remove_vector_type(::Type{Vector{T}}) where {T} = T
remove_vector_type(::Type{CategoricalVector{T}}) where {T} = CategoricalValue{T}
result_type(::Val{:Tuple}, args::TypeAst...) =
    Vector{Tuple{remove_vector_type.(result_type.(args))...}}

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt,
                         ::Val{:Tuple}, args::TypeAst...)
    tuple_columns = read_col_data.(Ref(sock), num_rows, args)
    return tuple.(tuple_columns...)
end

function write_col_data(sock::ClickHouseSock, data::AbstractVector{<:Tuple},
                        ::Val{:Tuple}, args::TypeAst...)
    isempty(data) && return
    length(data[1]) != length(args) &&
            error("Elements count mistmach between input and clickhouse tuples")

    for i in 1:length(args)
        write_col_data(sock, getindex.(data, i), args[i])
    end
end

ext: jl


text:
# %%
using Pkg
Pkg.activate(".") # Create new environment in this folder

<<<<<<>>>>>> c9842476bad2f8e1949a3d95931c4e62023d4134
# %% [markdown]
# We then list our imports

# %%
using Distributions
using LinearAlgebra
using Infiltrator
using StatsBase
using DataFrames
using Plots
using CategoricalArrays
using ShiftedArrays
using Chain
using DataFramesMeta
using LightGraphs
# using TikzGraphs
using Optim
using Random

include("./estimation-code.jl")
using .lewd

# %% [markdown]
# %%
struct define_parameters
    v_sd::Float64 # Firm amenity SD
    α_sd::Float64 # worker type sd
    ψ_sd::Float64 # Firm SD
    vψ_sd::Float64 # Cov(Amenity, Wage)
    cnetw::Float64 # network
    csig::Float64 # cross sectional sd
end

struct define_hyper_parameters
    nk::Int # firm types
    nl::Int # worker types
    nt::Int # time period
    ni::Int # n indiv
    λ::Float64 # moving prob
    δ::Float64 # death prob
    firms_per_type::Int
end

struct compute_transition_matrix
    ψ::Vector
    G::Array
    H::Array
    v::Array
    α::Vector

    function compute_transition_matrix(parameters::define_parameters, hyper_parameters::define_hyper_parameters)

        # setting parameters
        v_sd  = parameters.v_sd
        α_sd = parameters.α_sd
        ψ_sd  = parameters.ψ_sd
        vψ_sd = parameters.vψ_sd
        cnetw = parameters.cnetw
        csig  = parameters.csig

        nk = hyper_parameters.nk
        nl = hyper_parameters.nl
        firms_per_type = hyper_parameters.firms_per_type
        # approximate each distribution with some points of support
        # Firm fixed effect:
        function get_firm_characteristics(vψ_sd)
            """
            Input:
            vψ_sd: Covariance between amenities and firm FE.

            Output:
            ψ: Firm fixed effect
            v_a: Amenity value
            v: Total Firm value
            pr_job_offer: Probability of sending a job offer
            """
            Σ = [v_sd vψ_sd; vψ_sd ψ_sd]

            Φ = Normal(0,0.4)

            ψ = quantile.(Φ, (1:nk) / (nk + 1))

            v_a = quantile(Φ[:,1], (1:nk) / (nk + 1))

            # Value of the firm (This will be fixed over time):
            v = ψ .+ v_a # this should be something like firm FE + some amenity v_a

            pr_job_offer = pdf(Normal(0, 1), shuffle(v))

            pr_job_offer = pr_job_offer./sum(pr_job_offer)

            return ψ, v_a, v, pr_job_offer
        end

        # Individual type fixed effect (just for the wage equation):

        function get_individual_characteristics()
            """
            Output:
            α: individual fixed effects
            pr_death: Probability of death is random.
            """

            α = quantile.(Normal(), (1:nl) / (nl + 1))*α_sd

            pr_death = rand()

            return α, pr_death
        end

        ψ, v_a, v, pr_job_offer = get_firm_characteristics(vψ_sd)

        α, pr_death = get_individual_characteristics()

        # Let's create type-specific transition matrices
        # We are going to use joint normals centered on different values
        G = zeros(nl, nk, nk)

        for l in 1:nl, k in 1:nk
            # Transition probability weighted by pr of job offer (which is random)
            G[l, k, :] = pdf( Normal(0, csig), v .- cnetw * v[k])  # Get prob based on firm value
            G[l, k, :] = G[l, k, :] #.* pr_job_offer # Weight by pr of job offer
            G[l, k, k] = 0 # can't get offer from current place
            G[l, k, :] = G[l, k, :] ./ sum(G[l, k, :]) # Normalize probability
        end

        # We then solve for the stationary distribution over psis for each alpha value
        # We apply a crude fixed point approach
        H = ones(nl, nk) ./ nk

        new(ψ, G, H, v, α)

    end

end

struct simulation_draw

    df::DataFrame
    hyper_parameters::define_hyper_parameters
    transition_matrix::compute_transition_matrix
    kk::Matrix

    function simulation_draw(hyper_parameters::define_hyper_parameters, transition_matrix::compute_transition_matrix)

        α = transition_matrix.α
        ψ = transition_matrix.ψ
        G = transition_matrix.G
        H = transition_matrix.H
        v = transition_matrix.v

        # Setting parameters
        nk = hyper_parameters.nk
        nl = hyper_parameters.nl
        ni = hyper_parameters.ni
        nt = hyper_parameters.nt
        firms_per_type = hyper_parameters.fpt

        λ  = hyper_parameters.λ
        δ  = hyper_parameters.δ
        w_sigma = 0.2
        firms_per_type = hyper_parameters.firms_per_type

        # We simulate a balanced panel
        ll = zeros(Int64, ni, nt) # Worker type
        kk = zeros(Int64, ni, nt) # Firm type
        id = zeros(Int64, ni, nt) # Individual ID
        spellcount = zeros(Int64, ni, nt) # Employment spell

        id_count = 0

        for i in 1:ni   # Iterate over individuals

            id_count += 1

            # At time 1, we draw random type.
            l = rand(1:nl)
            ll[i,1] = l   # N x T

           # At time 1, we draw random firm, no stationary dist.
            kk[i,1] = sample(1:nk)

            # Insert id of the person:
            id[i,1] = id_count

            # Now iterate over time to figure where they are going
            for t in 2:nt

                kk_past = kk[i,t-1]

                if rand() > δ  # If doesn't die.

                    id[i,t] = id[i,t-1]   #keep the same id
                    ll[i,t] = l          #keep same worker type

                    if rand() < λ  #Consider λ which determines probability of receiving an offer

                        # Here person leaves...

                        # For now get job offer from firms at random
                        ws = fill(1.0/(nk-1), nk)
                        ws[kk_past] = 0
                        # kk_new = sample(1:nk, Weights(G[l, kk_past, :]))     # Get random firm type that makes the offer pr. G | kk[i,t-1] and weighted by vacancy rate.
                        kk_new = sample(1:nk, Weights(ws))     # Get random firm type that makes the offer pr. G | kk[i,t-1] and weighted by vacancy rate.
                        while kk_new == kk_past
                            kk_new = sample(1:nk, Weights(G[l, kk_past, :]))     # Get random firm type that makes the offer pr. G | kk[i,t-1] and weighted by vacancy rate.
                            kk_new = sample(1:nk, Weights(ws))     # Get random firm type that makes the offer pr. G | kk[i,t-1] and weighted by vacancy rate.
                        end
                        if v[kk_past] < v[kk_new] + rand(Gumbel(0, 1))
                            kk[i,t] = kk_new                                    # Otherwise, remain at same job
                            spellcount[i,t] = spellcount[i,t-1] + 1             # Add 1 to spell count
                        else
                            kk[i,t] = kk_past                                   # Otherwise, remain at same job
                            spellcount[i,t] = spellcount[i,t-1]                 # person didn't like the firm too much so he stays
                        end
                    else
                        kk[i,t] = kk_past                                       # Otherwise, remain at same job
                        spellcount[i,t] = spellcount[i,t-1]                     # Do not add to spell count

                    end

                else    # if dies

                    # Fill spot with new guy and guy id
                    id_count += 1
                    id[i,t] = id_count
                    #Draw new type for the guy
                    l = rand(1:nl)
                    ll[i,t] = l
                    # The worker starts at the same firm as the previous guy:
                    kk[i,t] = kk_past
                    spellcount[i,t] = 0
                end

            end

        end

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

        # This part likely to be the same:
        # we remove industries so jj == kk
        jj = kk # Firm identifiers

        contiguous_ids = Dict( unique(jj) .=> 1:length(unique(jj))  )
        jj .= getindex.(Ref(contiguous_ids),jj);

        #----------------------------------------------------------------------------------------
        ii = repeat(1:ni,1,nt)
        tt = repeat((1:nt)',ni,1)
        df = DataFrame(i=ii[:], j=jj[:], l=ll[:], k=kk[:], α=α[ll[:]], ψ=ψ[kk[:]], v = v[kk[:]], t=tt[:], w_id=id[:] , spell=spellcount[:]);

        df[!, :lw] = df.α + df.ψ + w_sigma * rand(Normal(), size(df)[1]);
        df[!, :v_a] = df.v - df.ψ ;

        new(df, hyper_parameters, transition_matrix, kk)

    end
end

# Get the next firm the mover is moving to
# by leading the firm column and extracting
# just the firm and next firm as one observation
function find_firm_links(mover_df::DataFrame)
    firm_link_df = @chain mover_df begin
        sort([:w_id, :t])
        groupby(:w_id)
        transform(:j => lead => :j_next)
        transform([:j, :j_next] => ByRow(==) => :j_j_next_BroadcastFunction)
        subset(:j_j_next_BroadcastFunction => x -> x .== false, skipmissing = true)
        select(:j, :j_next)
    end
    return firm_link_df
end

function create_M_flows(link_df)
    M = @chain link_df begin        # Use link_df, recall this is moves, so j =/= j_next
        groupby([:j, :j_next])      # Group by :j and next job
        combine(nrow => :count, [:j, :j_next] => ((x, y) -> ("M_" .* string.(y) .*"_" .* string.(x))) => :M)
    end
    return M
end

# %%
Random.seed!(1234);
sim_parameters = define_parameters(2.0, 1.0, 1.0, 0.1, 0.2, 0.2)
sim_hyper_parameters = define_hyper_parameters(4, 10, 10, 100_000, 0.8, 0.2, 1)
sim_transition_matrix = compute_transition_matrix(sim_parameters, sim_hyper_parameters)
# %%
sim_draw = simulation_draw(sim_hyper_parameters, sim_transition_matrix)
sim_df = sim_draw.df
df_connected = create_connected_df(sim_df)

# %%
link_df     = find_firm_links(find_movers(sim_df))
M_df        = create_M_flows(link_df)
S_kk, M_0   = create_fixed_point_matrices(sim_df)

# %%
connected_df = lewd.create_connected_df(sim_df)
link_df = lewd.find_firm_links(lewd.find_movers(connected_df))
M_df = lewd.create_M_flows(link_df)
S_kk, M_0 = lewd.create_fixed_point_matrices(connected_df)
S_kk

M_df

# %%
rank = lewd.estimate_rank(S_kk, M_0)
# %%
v_rank = @chain connected_df begin
    groupby(:j)
    combine(:v => unique)
    sort(:v_unique)
end

rank_df = DataFrame([1:length(rank), log.(rank)], [:j, :rank])

comp_df = innerjoin(v_rank, rank_df, on = :j)
comp_df[!, "s_kk"] = diag(S_kk)
sort(comp_df, :s_kk)
# %%
cor(comp_df.v_unique, comp_df.rank)
corspearman(comp_df.v_unique, comp_df.rank)
# %%
cor(comp_df.v_unique, comp_df.s_kk)
corspearman(comp_df.v_unique, comp_df.s_kk)
comp_df

# %%

function estimate_v(df)
    connected_df = lewd.create_connected_df(df)
    S_kk, M_0 = lewd.create_fixed_point_matrices(connected_df)
    v_hat = lewd.estimate_rank(S_kk, M_0)
    return v_hat
end

function sim_many_times(N)
    comp_vec = Vector{DataFrame}(undef, N)
    for i in 1:N
        sim_i = simulation_draw(sim_hyper_parameters, sim_transition_matrix)
        sim_i_df = sim_i.df

        v_rank = @chain connected_df begin
            groupby(:j)
            combine(:v => unique)
            sort(:v_unique)
        end
        v_hat = estimate_v(sim_i_df)
        rank_df = DataFrame([1:length(v_hat), v_hat], [:j, :rank])

        comp_df = innerjoin(v_rank, rank_df, on = :j)
        comp_df[!, "draw_n"] .= i
        comp_vec[i] = comp_df
    end
    return comp_vec
end

many_sim_df = reduce(vcat, sim_many_times(100))

cor_df = @chain many_sim_df begin
    groupby(:draw_n)
    combine([:v_unique, :rank] => ((x, y) -> cor(x, y)) => :cor,
            [:v_unique, :rank] => ((x, y) -> corspearman(x, y)) => :rank_cor)
end
histogram(cor_df.rank_cor, bins = 10)

# %%
function anon_function(param_list)
    # parameters(v_sd, \alpha_sd, ψ_sd, vψ_sd, cnetw, csig)
    parameters = parameters(param_list[1], param_list[2], param_list[3], 0.2, 0.2)
    values = variance_calibration(parameters, sim_hyper_parameters)
    target_values = [
        0.51, # variance of person effect
        0.14, # variance of employer effect
        0.10, # 2 cov(person, employer)
        0.19, # corr(person, employer)
        0.67, # variance log earnings

        0.400, # V^{EE} and \psi (pearson)
        0.530, # V^e and \psi (pearson)
        0.045, # V^{EE} log(size) (pearson)
        0.151, # V^e log(size)(pearson)
        0.093  # ψ log size (pearson)
    ]
    target = [0.084, 0.025, 0.003, 0.138]
    mse = mean((target .- values).^2)
    return mse
end

# %%|
using Optim
results = Optim.optimize(
    anon_function,
    [0.25, 0.25, 0.25, 0.25],
    # LBFGS(),
    NelderMead(),
    Optim.Options(iterations = 10_000)
)

# %%
# %%
using Optim
approx_results = Optim.optimize(
    anon_function,
    [0.25, 0.25, 0.25, 0.25],
    # LBFGS(),
    NelderMead(),
    Optim.Options(g_abstol = 1e-6)
)

approx_calibrated_parameters = parameters(
    Optim.minimizer(approx_results)[1],
    Optim.minimizer(approx_results)[2],
    Optim.minimizer(approx_results)[3],
    0.2,
    Optim.minimizer(approx_results)[4],
    0.2
    # Optim.minimizer(results)[5]
)

# %% [markdown]
# The optimisation function literally returns whatever the initial parameters are but still
# returns a success flag... nice. We <3 identification.
#

# %%
calibrated_transition_matrix  = compute_transition_matrix(calibrated_parameters, sim_hyper_parameters)
calibrated_sim = simulation_draw(sim_hyper_parameters, calibrated_transition_matrix)

values = variance_calibration(calibrated_parameters, sim_hyper_parameters)
target = [0.084, 0.025, 0.003, 0.137]
mse = mean((target .- values).^2)
# %% [markdown]
# Generating calibrated plots - these don't look totally convincing
# %%

p1 = plot(calibrated_transition_matrix.G[1, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[1, :, :]", st=:wireframe)
p2 = plot(calibrated_transition_matrix.G[end, :, :], xlabel="Previous Firm", ylabel="Next Firm", zlabel="G[nl, :, :]", st=:wireframe, right_margin = 10Plots.mm) # right_margin makes sure the figure isn't cut off on the right
plot(p1, p2, layout = (1, 2), size=[600,300])
# %%
plot(calibrated_transition_matrix.H, xlabel="Worker", ylabel="Firm", zlabel="H", st=:wireframe)

# %% [markdown]
# ## Getting connected set:

# %%
##### Ed attempt

# Iterate through our firm link df and create a matrix
# of links. In R I would model.matrix(a ~ b) but idk how to
# do that here

# %%
df_connected = create_connected_df(sim_df)

# %% AKM Estimation:

function akm_estimation(df::DataFrame)
    df[!, :α_hat] .= 0
    df[!, :ψ_hat] .= 0

    delta = Inf
    tol = 1e-5
    transform!(x -> mean(x.lw - x.ψ_hat), groupby(df, :i))
    df.α_hat = df.x1
    transform!(x -> mean(x.lw - x.α_hat), groupby(df, :j))
    df.ψ_hat = df.x1
    mse = mean((df.lw .- df.α_hat .- df.ψ_hat).^2)

    while delta > tol

        transform!(x -> mean(x.lw - x.ψ_hat), groupby(df, :i))
        df.α_hat = df.x1
        transform!(x -> mean(x.lw - x.α_hat), groupby(df, :j))
        df.ψ_hat = df.x1
        new_mse = mean((df.lw .- df.α_hat .- df.ψ_hat).^2)
        delta =  new_mse - mse
        mse = new_mse
    end
    return df
end

# %%
df_connected_results = akm_estimation(df_connected);
# %%

# Q8: Limited mobility bias
#---------------------------------------------------------------------
# %%

function run_mob_bias_simulation(λ::Float64, nt::Int, parameters::parameters)
        λ_hyper_parameters = hyper_parameters(30, 10, λ, nt, 10_000) # gen data changing λ and nt

        λ_transition_matrix = compute_transition_matrix(parameters, λ_hyper_parameters)

        λ_sim = simulation_draw(λ_hyper_parameters, λ_transition_matrix)
        λ_df = λ_sim.df
        λ_df_connected = create_connected_df(λ_df)

        λ_df_connected_results = akm_estimation(λ_df_connected)

        cov_αψ = cov(λ_df_connected_results.α, λ_df_connected_results.ψ)
        cov_αψ_hat = cov(λ_df_connected_results.α_hat,λ_df_connected_results.ψ_hat )

        cov_bias = cov_αψ - cov_αψ_hat
        var_ψ = var(λ_df_connected_results.ψ)
        var_ψ_hat = var(λ_df_connected_results.ψ_hat)
        var_bias = var_ψ - var_ψ_hat
        result_df = DataFrame(
            Dict(
                :cov_αψ => cov_αψ,
                :cov_αψ_hat => cov_αψ_hat,
                :cov_bias => cov_bias,
                :var_ψ => var_ψ,
                :var_ψ_hat => var_ψ_hat,
                :var_bias => var_bias,
                :T => nt,
                :λ => λ
            )
        )

        return  result_df
end

# %%

λ_list = [0.1, 0.3, 0.5, 0.7, 0.9]
nt_list = [5, 6, 8, 10, 15];

λ_bias_df = broadcast(
    x -> run_mob_bias_simulation(x, nt_list[2], sim_parameters),
    λ_list
)

# %%

λ_bias_plot_df = reduce(vcat, λ_bias_df)
plot(
    λ_bias_plot_df[!, :λ],
     λ_bias_plot_df[!, :cov_αψ],
     markershape= :circle,
     linetype = :scatter)
plot!(
    λ_bias_plot_df[!, :λ],
    λ_bias_plot_df[!, :cov_αψ_hat])

# %%
plot(
   λ_bias_plot_df.λ,
   λ_bias_plot_df.var_bias
)

# %%

nt_bias_df = broadcast(
    x -> run_mob_bias_simulation(λ_list[1], x, sim_parameters),
    nt_list
)

nt_bias_df = reduce(vcat, nt_bias_df)

# %%

plot(
    nt_bias_df[!, :T],
     nt_bias_df[!, :cov_αψ],
     markershape= :circle,
     linetype = :scatter)
plot!(
    nt_bias_df[!, :T],
    nt_bias_df[!, :cov_αψ_hat])

# %%
plot(
   nt_bias_df.T,
   nt_bias_df.var_bias
)

#
# %% [markdown]
# ## Limited mobility bias
# %% [markdown]
# We now have everything we need to look at the impact of limited mobility bias. Compute the following:
#
# 1. Compute the estimated variance of firm FE
# 2. Do it for varying levels of mobility λ. For each the number of movers, collect the actual variance and the estimated variance. Run it for different panel lengths: 5, 6, 8, 10, 15.
#
# Question 8
#
#  - Report this in a plot. This should look like the [Andrews et al.](https://www.sciencedirect.com/science/article/pii/S0165176512004272) plot.
# %% [markdown]
# ## Correcting the bias
# %% [markdown]
# Question 9
#
#  - Implement both the exact as well as the approximated (Hutchkinson trace approximation) Andrews et al. type fixed-effect bias correction based on homoskedacity. Redo the plot from Question 6 so that it includes these 2 estimates.
# %% [markdown]
# ## Evidence of learning
#
# Question 10
#
#  - Regress the wage of a worker at time $t$ on his wage at time $t-l$ and the average wage of his co-workers at time $t-l$ for some lags $l$ focusing on workers that did change firms between $t-l$ and $t$.
#  - Comment on the result.
#


text:
module Iterating
using ..CSTParser: EXPR, headof, hastrivia, isoperator, valof, isstringliteral, is_exor, is_lparen, is_rparen

function Base.getindex(x::EXPR, i)
    try
        a = _getindex(x, i)
        if a === nothing
            error("indexing error for $(x.head) expression at $i")
        end
        return a
    catch
        rethrow()
        error("indexing error for $(x.head) expression at $i. args: $(x.args !== nothing ? headof.(x.args) : []) trivia: $(x.trivia !== nothing ? headof.(x.trivia) : [])")
    end
end

function _getindex(x::EXPR, i)
    if headof(x) === :abstract
        _abstract(x, i)
    elseif headof(x) === :as
        odda_event(x, i)
    elseif headof(x) === :block
        _block(x, i)
    elseif headof(x) === :braces
        _braces(x, i)
    elseif headof(x) === :brackets
        oddt_evena(x, i)
    elseif headof(x) === :call
        _call(x, i)
    elseif headof(x) === :comparison || headof(x) === :file
        x.args[i]
    elseif headof(x) === :comprehension
        tat(x, i)
    elseif headof(x) === :const
        _const(x, i)
    elseif headof(x) === :curly
        _curly(x, i)
    elseif headof(x) === :do
        odda_event(x, i)
    elseif headof(x) === :elseif
        _elseif(x, i)
    elseif headof(x) === :errortoken
        x.args[i]
    elseif headof(x) === :export
        oddt_evena(x, i)
    elseif headof(x) === :filter
        _filter(x, i)
    elseif headof(x) === :flatten
        _flatten(x, i)
    elseif headof(x) === :for
        taat(x, i)
    elseif headof(x) === :function || headof(x) === :macro
        _function(x, i)
    elseif headof(x) === :generator
        odda_event(x, i)
    elseif headof(x) === :global || headof(x) === :local
        _global(x, i)
    elseif headof(x) === :if
        if isoperator(first(x.trivia)) # ternary op
            odda_event(x, i)
        else
            _if(x, i)
        end
    elseif headof(x) === :kw
        _kw(x, i)
    elseif headof(x) === :let
        taat(x, i)
    elseif headof(x) === :return
        oddt_evena(x, i)
    elseif headof(x) === :macrocall
        _macrocall(x, i)
    elseif headof(x) === :module
        _module(x, i)
    elseif headof(x) === :outer
        ta(x, i)
    elseif headof(x) === :parameters
        if length(x.args) > 1 && headof(x.args[1]) === :parameters
            if i == length(x)
                x.args[1]
            elseif iseven(i)
                x.trivia[div(i, 2)]
            else
                x.args[div(i + 1, 2) + 1]
            end
        elseif length(x.args) > 1 && headof(x.args[2]) === :parameters
            if i == length(x)
                x.args[2]
            elseif i == 1
                x.args[1]
            elseif isodd(i)
                x.args[div(i + 1, 2) + 1]
            else
                x.trivia[div(i, 2)]
            end
        else
            odda_event(x, i)
        end
        # if hastrivia(x)
        #     odda_event(x, i)
        # else
        #     x.args[i]
        # end
    elseif headof(x) === :primitive
        _primitive(x, i)
    elseif headof(x) === :quote
        _quote(x, i)
    elseif headof(x) === :quotenode
        _quotenode(x, i)
    elseif headof(x) === :ref
        _call(x, i)
    elseif headof(x) === :row || headof(x) === :nrow
        x.args[i]
    elseif headof(x) === :string
        _string(x, i)
    elseif headof(x) === :struct
        _struct(x, i)
    elseif headof(x) === :toplevel
        x.args[i]
    elseif headof(x) === :try
        _try(x, i)
    elseif headof(x) === :tuple
        _tuple(x, i)
    elseif headof(x) === :typed_comprehension
        odda_event(x, i)
    elseif headof(x) === :typed_vcat || headof(x) === :typed_hcat || headof(x) === :typed_ncat
        _typed_vcat(x, i)
    elseif headof(x) === :using || headof(x) === :import
        _using(x, i)
    elseif headof(x) === :vcat || headof(x) === :hcat || headof(x) === :ncat || headof(x) === :bracescat
        _vcat(x, i)
    elseif headof(x) === :vect
        _vect(x, i)
    elseif headof(x) === :where
        _where(x, i)
    elseif headof(x) === :while
        taat(x, i)
    elseif isoperator(headof(x))
        if valof(headof(x)) == ":"
            _colon_in_using(x, i)
        elseif valof(headof(x)) == "."
            _dot(x, i)
        elseif valof(headof(x)) == "=" && hastrivia(x)
            # a loop that has been lowered from in/∈ to = at parse time
            if i == 1
                x.args[1]
            elseif i == 2
                x.trivia[1]
            elseif i == 3
                x.args[2]
            end
        elseif valof(headof(x)) == "->" && headof(x).fullspan == 0
            # anon function as used in a 'do' block
            if i == 1
                x.args[1]
            else
                x.args[2]
            end
        elseif length(x) == 2 && (valof(headof(x)) == "\$" || valof(headof(x)) == "::" || valof(headof(x)) == "<:" || valof(headof(x)) == ">:"  || valof(headof(x)) == "&")
            if i == 1
                x.head
            else
                x.args[1]
            end
        elseif !hastrivia(x)
            if i == 1
                x.args[1]
            elseif i == 2
                x.head
            elseif i == 3
                x.args[2]
            end
        else
            if i == 1
                x.head
            elseif i == 2
                x.trivia[1]
            elseif i == length(x)
                last(x.trivia)
            elseif isodd(i)
                x.args[div(i+1, 2) - 1]
            else
                x.trivia[div(i, 2)]
            end
        end
    else
        error("Indexing $x at $i")
    end
end
Base.iterate(x::EXPR) = length(x) == 0 ? nothing : (x[1], 1)
Base.iterate(x::EXPR, s) = s < length(x) ? (x[s + 1], s + 1) : nothing
Base.firstindex(x::EXPR) = 1
Base.lastindex(x::EXPR) = x.args === nothing ? 0 : length(x)

# Base.setindex!(x::EXPR, val, i) = Base.setindex!(x.args, val, i)
# Base.first(x::EXPR) = x.args === nothing ? nothing : first(x.args)
# Base.last(x::EXPR) = x.args === nothing ? nothing : last(x.args)

function ta(x, i)
    if i == 1
        x.trivia[1]
    elseif i == 2
        x.args[1]
    end
end
function tat(x, i)
    if i == 1
        x.trivia[1]
    elseif i == 2
        x.args[1]
    elseif i == 3
        x.trivia[2]
    end
end

function taat(x, i)
    if i == 1
        x.trivia[1]
    elseif i == 2
        x.args[1]
    elseif i == 3
        x.args[2]
    elseif i == 4
        x.trivia[2]
    end
end

function oddt_evena(x, i)
    if isodd(i)
        x.trivia[div(i + 1, 2)]
    else
        x.args[div(i, 2)]
    end
end

function odda_event(x, i)
    if isodd(i)
        x.args[div(i + 1, 2)]
    else
        x.trivia[div(i, 2)]
    end
end

function _const(x, i)
    if length(x.trivia) === 1
        ta(x, i)
    elseif length(x.trivia) === 2
        #global const
        if i < 3
            x.trivia[i]
        elseif i == 3
            x.args[1]
        end
    end
end

function _global(x, i)
    if hastrivia(x) && (headof(first(x.trivia)) === :GLOBAL || headof(first(x.trivia)) === :LOCAL)
        oddt_evena(x, i)
    else
        odda_event(x, i)
    end
end

function _abstract(x, i)
    if i < 3
        x.trivia[i]
    elseif i == 3
        x.args[1]
    elseif i == 4
        x.trivia[3]
    end
end

function _primitive(x, i)
    if i < 3
        x.trivia[i]
    elseif i == 3
        x.args[1]
    elseif i == 4
        x.args[2]
    elseif i == 5
        x.trivia[3]
    end
end

function _struct(x, i)
    if length(x.trivia) == 2
        if i == 1
            x.trivia[1]
        elseif 1 < i < 5
            x.args[i - 1]
        elseif i == 5
            x.trivia[2]
        end
    elseif length(x.trivia) == 3 # mutable
        if i < 3
            x.trivia[i]
        elseif 2 < i < 6
            x.args[i - 2]
        elseif i == 6
            x.trivia[3]
        end
    end
end

function _block(x, i)
    if hastrivia(x) # We have a begin block
        if headof(x.trivia[1]) === :BEGIN
            if i == 1
                x.trivia[1]
            elseif 1 < i < length(x)
                x.args[i - 1]
            elseif i == length(x)
                x.trivia[2]
            end
        elseif headof(x.trivia[1]) === :COMMA # comma sep list as in :let block or iterators
            odda_event(x, i)
        end
    else
        x.args[i]
    end
end

function _quote(x, i)
    if x.trivia !== nothing && length(x.trivia) == 1
        if i == 1
            x.trivia[1]
        elseif i == 2
            x.args[1]
        end
    elseif x.trivia !== nothing && length(x.trivia) == 2
        tat(x, i)
    elseif x.trivia === nothing
        x.args[i]
    else
        error()
    end
end

function _quotenode(x, i)
    if hastrivia(x)
        if i == 1
            x.trivia[1]
        elseif i == 2
            x.args[1]
        end
    elseif i == 1
        x.args[1]
    end
end

function _function(x, i)
    if length(x.args) == 1
        tat(x, i)
    else length(x.args) == 2
        taat(x, i)
    end
end

function _braces(x, i)
    if length(x.args) > 0 && headof(x.args[1]) === :parameters
        if i == 1
            x.trivia[1]
        elseif i == length(x)
            last(x.trivia)
        elseif i == length(x) - 1
            x.args[1]
        elseif iseven(i)
            x.args[div(i, 2) + 1]
        else
            x.trivia[div(i + 1, 2)]
        end
    else
        if i == length(x)
            last(x.trivia)
        else
            oddt_evena(x, i)
        end
    end
end

function _curly(x, i)
    if i == 1
        x.args[1]
    elseif i == length(x)
        last(x.trivia)
    elseif length(x.args) > 1 && headof(x.args[2]) === :parameters
        if i == length(x) - 1
            x.args[2]
        elseif isodd(i)
            x.args[div(i + 1, 2) + 1]
        else
            x.trivia[div(i, 2)]
        end
    else
        if isodd(i)
            x.args[div(i + 1, 2)]
        else
            x.trivia[div(i, 2)]
        end
    end
end

function _using(x, i)
    oddt_evena(x, i)
end

function _colon_in_using(x, i)
    if i == 1
        x.args[1]
    elseif i == 2
        x.head
    elseif isodd(i)
        x.args[div(i + 1, 2)]
    else
        x.trivia[div(i, 2) - 1]
    end
end

function count_continuous(f, itr)
    cnt = 0
    for x in itr
        !f(x) && break
        cnt += 1
    end
    return cnt
end

function _dot(x, i)
    if x.head.span == 0 # Empty dot op, in using statement
        ndots = count_continuous(a -> valof(a) == ".", x.args)
        if i <= ndots
            x.args[i]
        elseif iseven(ndots)
            if isodd(i)
                x.args[ndots + div(i - ndots + 1, 2)]
            else
                x.trivia[div(i - ndots, 2)]
            end
        else
            if iseven(i)
                x.args[ndots + div(i - ndots + 1, 2)]
            else
                x.trivia[div(i - ndots, 2)]
            end
        end
    elseif length(x) == 2
        if i == 1
            x.head
        elseif i == 2
            x.args[1]
        end
    else
        if i == 1
            x.args[1]
        elseif i == 2
            x.head
        elseif i == 3
            x.args[2]
        end
    end
end

function _call(x, i)
    if isoperator(x.args[1])
        if length(x) == 2 # unary op call
            x.args[i]
        elseif length(x) == 3 && !hastrivia(x)# binary
            if i == 1
                x.args[2]
            elseif i == 2
                x.args[1]
            elseif i == 3
                x.args[3]
            end
        elseif hastrivia(x) && headof(x.trivia[1]) === :LPAREN # op call w/ brackets
            _curly(x, i)
        else # chained call: a + b + c + d + ...
            if i == 1
                x.args[2]
            elseif i == 2
                x.args[1]
            elseif i == 3
                x.args[3]
            elseif hastrivia(x)
                if iseven(i)
                    x.trivia[div(i - 2, 2)]
                else
                    x.args[div(i - 3, 2) + 3]
                end
            end
        end
    elseif hastrivia(x)
        _curly(x, i)
    else
        x.args[i]
    end
end

function _kw(x, i)
    if i == 1
        x.args[1]
    elseif i == 2
        x.trivia[1]
    elseif i == 3
        x.args[2]
    end
end

function _tuple(x, i)
    hasparams = x.args !== nothing && length(x.args) > 0 && headof(x.args[1]) === :parameters
    if hasparams
        if i == 1
            first(x.trivia)
        elseif i == length(x)
            last(x.trivia)
        elseif i == length(x) - 1
            first(x.args)
        elseif isodd(i)
            x.trivia[div(i + 1, 2)]
        else
            x.args[div(i, 2) + 1]
        end
    else
        if isempty(x.args)
            x.trivia[i]
        elseif hastrivia(x)
            if first(x.trivia).head === :LPAREN
                if i == length(x)
                    last(x.trivia)
                else
                    oddt_evena(x, i)
                end
            else
                odda_event(x, i)
            end
        else
            x.args[i]
        end
        # if isempty(x.args)
        #     x.trivia[i]
        # elseif length(x.trivia) == length(x.args) - 1  # No brackets, no trailing comma
        #     odda_event(x, i)
        # elseif length(x.trivia) - 1 == length(x.args) # Brackets, no trailing comma
        #     oddt_evena(x, i)
        # elseif length(x.trivia) - 2 == length(x.args) # Brackets, trailing comma
        #     if i == length(x)
        #         last(x.trivia)
        #     elseif i == length(x) - 1
        #         x.trivia[end-1]
        #     else
        #         oddt_evena(x, i)
        #     end
        # else
        #     odda_event(x, i)
        # end
    end
end

function _if(x, i)
    if length(x) == 4 # if c expr end
        taat(x, i)
    elseif length(x) == 5 # if c expr elseif... end
    if i == 1
        x.trivia[1]
    elseif i == 2
        x.args[1]
    elseif i == 3
        x.args[2]
    elseif i == 4
        x.args[3]
    elseif i == 5
        x.trivia[2]
    end
    elseif length(x) == 6 # if c expr else expr end
        if i == 1
            x.trivia[1]
        elseif i == 2
            x.args[1]
        elseif i == 3
            x.args[2]
        elseif i == 4
            x.trivia[2]
        elseif i == 5
            x.args[3]
        elseif i == 6
            x.trivia[3]
        end
    end
end

function _elseif(x, i)
    if length(x) == 3 || length(x) == 5
        if i == 1
            x.trivia[1]
        elseif i == 2
            x.args[1]
        elseif i == 3
            x.args[2]
        elseif i == 4
            x.trivia[2]
        elseif i == 5
            x.args[3]
        end
    elseif length(x) == 4
        if i == 1
            x.trivia[1]
        else
            x.args[i - 1]
        end
    end
end

function _string(x, i)
    # TODO: this is mega slow
    ai, ti = 1,1
    arg = isstringliteral(x.args[1])# && sizeof(valof(x.args[1])) !== x.args[1].fullspan
    isinterpolant = !arg
    bracket = false
    if hastrivia(x) && (x.trivia[1].head === :STRING || x.trivia[1].head === :TRIPLESTRING) && isempty(x.trivia[1].val)
        if i == 1
            return x.trivia[1]
        end
        arg = false
    end
    if i == length(x) && hastrivia(x) && (last(x.trivia).head === :STRING || last(x.trivia).head === :TRIPLESTRING || last(x.trivia).head === :errortoken) && isempty(last(x.trivia).val)
        return last(x.trivia)
    end
    for j = 1:i
        if j == i
            return arg ? x.args[ai] : x.trivia[ti]
        end

        if arg
            ai += 1
            if isinterpolant
                arg = !bracket
                if ai <= length(x.args) && !(isstringliteral(x.args[ai]) || x.args[ai].head === :errortoken)
                    # interpolated value immediately followed by xor
                    arg = false
                end
                bracket = false
                isinterpolant = false

            else
                arg = false
            end
        else
            if is_exor(x.trivia[ti]) || (x.trivia[ti].head === :errortoken && x.trivia[ti].args[1].head === :OPERATOR && valof(x.trivia[ti].args[1]) == "\$")
                if ti < length(x.trivia) && is_lparen(x.trivia[ti + 1])
                    #
                else
                    isinterpolant = true
                    arg = true
                end
            elseif is_lparen(x.trivia[ti])
                isinterpolant = true
                arg = true
                bracket = true
            elseif is_rparen(x.trivia[ti])
                isinterpolant = false
                arg = !(ai <= length(x.args) && !(isstringliteral(x.args[ai])) || x.args[ai].head == :errortoken)
                bracket = false
            end
            ti += 1
        end
    end
end

function _macrocall(x, i)
    if !hastrivia(x)
        x.args[i]
    else
        if i < 3
            x.args[i]
        elseif i == length(x)
            last(x.trivia)
        elseif length(x.args) > 2 && headof(x.args[3]) === :parameters
            if i == length(x) - 1
                x.args[3]
            elseif i == length(x)
                last(x.trivia)
            elseif isodd(i)
                x.trivia[div(i, 2)]
            else
                x.args[div(i + 1, 2) + 2]
            end
        else
            if isodd(i)
                x.trivia[div(i, 2)]
            else
                x.args[div(i + 1, 2) + 1]
            end
        end
    end
end

function _vect(x, i)
    _braces(x, i)
end

function _vcat(x, i)
    if i == 1
        x.trivia[1]
    elseif i == length(x)
        x.trivia[2]
    else
        x.args[i - 1]
    end
end

function _typed_vcat(x, i)
    if i == 1
        x.args[1]
    elseif i == 2
        x.trivia[1]
    elseif i == length(x)
        x.trivia[2]
    else
        x.args[i - 1]
    end

end

function _module(x, i)
    if i == 1
        x.trivia[1]
    elseif i == 2
        x.args[1]
    elseif i == 3
        x.args[2]
    elseif i == 4
        x.args[3]
    elseif i == 5
        x.trivia[2]
    end
end

function _filter(x, i)
    if i == length(x)
        x.args[1]
    elseif iseven(i)
        x.trivia[div(i, 2)]
    else
        x.args[div(i + 1, 2) + 1]
    end
end

function _try(x, i)
    if i == 1
        x.trivia[1]
    elseif i == 2
        x.args[1]
    elseif length(x) == 6 # try expr catch e expr end
        if i == 3
            x.trivia[2]
        elseif i == 4
            x.args[2]
        elseif i == 5
            x.args[3]
        elseif i == 6
            x.trivia[3]
        end
    elseif length(x) == 7 # try expr catch finally end
        if i == 3
            x.args[2]
        elseif i == 4
            x.args[3]
        elseif i == 5
            x.trivia[2]
        elseif i == 6
            x.args[4]
        elseif i == 7
            x.trivia[3]
        end
    elseif length(x) == 8 # 'try expr catch e expr finally expr end' or 'else end'
        if i == 3
            x.trivia[2]
        elseif i == 4
            x.args[2]
        elseif i == 5
            x.args[3]
        elseif i == 6
            x.trivia[3]
        elseif i == 7
            x.args[4]
        elseif i == 8
            x.trivia[4]
        end
    elseif length(x) == 10 # try expr catch e expr else expr finally expr end
        if i == 3
            x.trivia[2]
        elseif i == 4
            x.args[2]
        elseif i == 5
            x.args[3]
        elseif i == 6
            x.trivia[4]
        elseif i == 7
            x.args[5]
        elseif i == 8
            x.trivia[3]
        elseif i == 9
            x.args[4]
        elseif i == 10
            x.trivia[5]
        end
    end
end

function _where(x, i)
    if i == 1
        x.args[1]
    elseif i == 2
        x.trivia[1]
    elseif i == 3 && length(x) == 3 && length(x.trivia) == 1
        return x.args[2]
    elseif i == length(x)
        last(x.trivia)
    elseif length(x.args) > 1 && headof(x.args[2]) === :parameters
        if i == length(x) - 1
            x.args[2]
        elseif isodd(i)
            x.trivia[div(i + 1, 2)]
        else
            x.args[div(i, 2) + 1]
        end
    else
        oddt_evena(x, i)
    end
end

function _flatten(x, i)
    lhs = _flatten_lhs(x)
    lhs[i]
end

function _flatten_lhs(x, ret = [])
    if x.args[1].head === :generator || x.args[1].head === :flatten
        if headof(x) !== :flatten
            for i = 2:length(x)
                push!(ret, x[i])
            end
        end
        _flatten_lhs(x.args[1], ret)
    else
        for i = 2:length(x)
            push!(ret, x[i])
        end
        pushfirst!(ret, x.args[1])
    end
end

end


text:
"""

    mpnmodel(S::Int64, Co::Float64, forbidden::Float64)

Return `UnipartiteNetwork` with links assigned according to minimum
potential niche model for given number of `S`, connectivity `Co` and
probability of `forbidden` link occurence.

> Allesina, S., Alonso, D. and Pascual, M. (2008) ‘A General Model for Food Web
> Structure’, Science, 320(5876), pp. 658–661. doi: 10.1126/science.1156269.

See also: `nichemodel`, `cascademodel`, `nestedhierarchymodel`
"""
function mpnmodel(S::Int64, Co::Float64, forbidden::Float64)

    Co >= 0.5 && throw(ArgumentError("The connectance cannot be larger than 0.5"))

    # Beta distribution parameter
    β = 1.0/(2.0*Co)-1.0

    # Pre-allocate the network
    A = UnipartiteNetwork(zeros(Bool, (S, S)))

    # Generate body size
    n = sort(rand(Uniform(0.0, 1.0), S))

    # Pre-allocate centroids
    c = zeros(Float64, S)

    # Generate random ranges
    r = n .* rand(Beta(1.0, β), S)

    # Generate random centroids
    for s in 1:S
        c[s] = rand(Uniform(r[s]/2, n[s]))
    end

    # The smallest species has a body size and range of 0
    for small_species_index in findall(x -> x == minimum(n), n)
        n[small_species_index] = 0.0
        r[small_species_index] = 0.0
    end

    # Testing
    # -------------------------------------------------------------------------
    # ranges = Dict{Int64, Tuple}()
    # -------------------------------------------------------------------------

    for consumer in 1:S

        # For testing
        # ---------------------------------------------------------------------
        # u = c[consumer] + r[consumer]
        # l = c[consumer] - r[consumer]
        # diet = findall(x -> (x > l) & (x < u), n)
        # #println("Species ", diet, " fall into diet of consumer ", consumer)
        # length(diet) != 0 && (ranges[consumer] = (diet[1], diet[end]))
        # ---------------------------------------------------------------------

        for resource in 2:(S-1)

            lower = (n[resource] > (c[consumer] - r[consumer]))
            lowerminus = (n[resource-1] > (c[consumer] - r[consumer]))
            upper = n[resource] < (c[consumer] + r[consumer])
            upperplus = n[resource+1] < (c[consumer] + r[consumer])

            if (upper & lower)

                (rand() < (1 - forbidden)) && (A[consumer, resource] = true)

                # Take care of first and last resource if they belong to the range
                ((resource-1 == 1) & lowerminus) && (A[consumer, resource-1] = true)
                ((resource+1 == S) & upperplus) && (A[consumer, resource+1] = true)

                # Edges of range
                lowerminus || (A[consumer, resource] = true)
                upperplus || (A[consumer, resource] = true)

            end

        end

    end

    return A

end

"""

    mpnmodel(N::T) where {T<: UnipartiteNetwork}

Applied to `UnipartiteNetwork` return its randomized version.
"""
function mpnmodel(N::T) where {T<: UnipartiteNetwork}

    # Estimate proportion of forbidden links

    ## Find ordering
    ## sort(N)

    ## Get the value
    forbidden = 0

    # Simulate the network
    return mpn(species(N), connectance(N), forbidden)

end

"""

    mpnmodel(parameters::Tuple)

Parameters tuple can also be provided in the form (S::Int64, Co::Float64,
forbidden::Float64).
"""
function mpnmodel(parameters::Tuple)
    return mpnmodel(parameters[1], parameters[2],parameters[3])
end

"""

    mpnmodel(mpnmodel(S::Int64, L::Int64, forbidden::Float64))

Average number of links can be specified instead of connectance.

"""
function mpnmodel(S::Int64, L::Int64, forbidden::Float64)

    Co = (L/(S*S))
    return mpnmodel(S, Co, forbidden)

end


text:
# Evolve all single output 2-input or 3-input goals, or evolve a chosen number of random 4-input goals.

# Example run:
# >  julia -p 4  -L EvolveGoals.jl -e "run_evolve_goals(\"prun\",3,16,20,3,4)"
# Runs on 4 cores with numinputs=3, random number seed=16, genome size=20, 3 runs/goal, 4 random goals.
# Output to the file  out/prun3_16_20_3_4.csv

# Evolves all single-output 2-input and single-output 3-input goals if numinputs is set to 2 or 3
#      and number_random_goals has the default value of 0.
# Evolves random single-ouput 4-input goals if numimputs is set to 4.
#      In this case, number_random_goals should be set in the call to run_evolve_goals.
# The chromosome length (numlevels) and the number of runs per goal can be specified.
# The random number seed can also be specified.
# To evolve all 2-input goals, run run_evolve_goals() with the second argument equal to 2
# To evolve all 3-input goals, run run_evolve_goals() with the second argument equal to 3
# To run using multiple processors, start Julia using the "-p" option, such as "julia -p 4"
#   to start with 4 processors
include("../src/CartesianGP.jl")
using CartesianGP

# macro similar to the @time macro except that it writes to the ost stream rather than the stdout stream.
# This definition must precede its use below.
macro otime(ex)
           quote
               local t0 = time()
               local val = $(esc(ex))
               local t1 = time()
               println(ost,"elapsed time: ", t1-t0, " seconds")
               val
           end
       end

# A high-level function to run function evolve_goals.
# it sets the summary to the filename, and max_gens to default values
# First example:  run_evolve_goals("prun",2,13,10,20)
#   calls "evolve_goals" with numinputs=2, rseed=13,numlevels=10,runs_per_goal=20
#   and the output going to the file: prun2_13_10_20.csv
# Second example:  run_evolve_goals("prun",4,13,100,5,10)
#   calls "evolve_goals" with numinputs=4, rseed=13,numlevels=100,runs_per_goal=5,number_random_goals=10
#   and the output going to the file: prun2_13_10_20.csv
function run_evolve_goals(filename_prefix,numinputs=2,rseed=2,num_levels=2,runs_per_goal=20,number_random_goals=0)
    directory = "out/"  # subdirectory for output CSV file
    filetype = ".csv"   # file extension
    underscore = "_"
    if number_random_goals == 0
        filename = "$directory$filename_prefix$numinputs$underscore$rseed$underscore$num_levels$underscore$runs_per_goal$filetype"
    else
        filename = "$directory$filename_prefix$numinputs$underscore$rseed$underscore$num_levels$underscore$runs_per_goal$underscore$number_random_goals$filetype"
    end
    ost = open(filename,"w")
    if numinputs == 2
        max_gens = 10000
    elseif numinputs == 3
        max_gens = 20000
    else   # numinputs == 4 or higher
        max_gens = 70000
    end
    params = Params(numinputs,1,1,num_levels,num_levels,raman_funcs)
    for id in procs()
        remotecall(id,set_globals,params,max_gens)
    end
    @otime ost = evolve_goals(ost,filename,params,runs_per_goal,max_gens,number_random_goals,rseed)
    close(ost)
end

# Sets the global variables params and max_gens
function set_globals(pp,mg)
    global params = pp
    global max_gens = mg
end

# This is the function used by map or pmap to evolve each goal in goal_list
function p_mu_lambda(g)
    global params
    global max_gens
    (ch,gens) = mu_lambda(params,g,max_gens)
    n = ch.number_active_nodes
    (g.truth_table,gens,n)
end

# Primitive functions used in this simulation.  Same as those used by Raman and Wagner.
raman_funcs = [AND, OR, XOR, NAND, NOR]

# Computes the parameters used for this example
function Params(numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs )
    targetfitness = 1.0
    fitfunc = hamming_max
    mu = 1
    lambda = 4
    mutrate = 0.05
    return Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
end

# Computes the goal list which includes repeated elements for multiple trials per goal
function goal_list_setup(p,runs_per_goal,maxgens,number_random_goals=0)
    global max_gens = maxgens
    num_goals = 0
    if p.numinputs <= 3
        if number_random_goals == 0
            num_goals = 2^2^p.numinputs   # for numinputs == 2 or 3, evolve all goals
            goal_list = [Goal(p.numinputs,(convert(BitString,div(i,runs_per_goal)),)) for i in 0:num_goals*runs_per_goal-1]
        else # if number_random_goals != 0 then evolve random goals
            num_goals = number_random_goals
            BitString_list = [convert(BitString,rand(0:2^2^p.numinputs-1)) for i in 1:num_goals]
            goal_list = [Goal(p.numinputs,(BitString_list[div(i,runs_per_goal)+1],)) for i in 0:num_goals*runs_per_goal-1]
        end
    else
        if number_random_goals == 0
            num_goals = 8  # a default value that is computationally feasible for numinputs==4.
        else
            num_goals = number_random_goals
        end
        BitString_list = [convert(BitString,rand(0:2^2^p.numinputs-1)) for i in 1:num_goals]
        goal_list = [Goal(p.numinputs,(BitString_list[div(i,runs_per_goal)+1],)) for i in 0:num_goals*runs_per_goal-1]
    end
    return num_goals,goal_list
end

# Computes averages array which has rows corresponding to goals and 3 columns:
#    1. average number of generations per goal (not counting runs where the goal was not evolved)
#    2. average number of active nodes used per goal
#    3. number of successful (generations less than max_gens and goal was evolved) runs per goal
function averages(ng,runs_per_goal,max_gens,r)
    # Convert r which is a 1-D array with multiple entries for the runs for each goal
    #   to a 2-D array with runs_per_goal rows and num_goals columns.
    s = reshape(r,runs_per_goal,ng)
    average_array = zeros(ng,3)
    runs_per_goal_array = zeros(ng)
    for j in 1:ng
        for i in 1:runs_per_goal
            # k==1 is generations, k==2 is number of active nodes
            if s[i,j][2] == max_gens
                # do not add to average_array[j,1]
                average_array[j,2] += s[i,j][3]
                # do not add to runs_per_goal_array[j]
            else
                runs_per_goal_array[j] += 1
                for k = 1:2
                    average_array[j,k] += s[i,j][k+1]
                end
            end
        end
    end
    for j in 1:ng
        if runs_per_goal_array[j] > 0
            average_array[j,1] = average_array[j,1]/runs_per_goal_array[j]
        else
            average_array[j,1] = 0
        end
        average_array[j,2] = average_array[j,2]/runs_per_goal
        average_array[j,3] = runs_per_goal_array[j]  # store number of successful runs in column 3
    end
    average_array
end

# Runs "mu_lambda" or "evolve" to evolve the goals on goal_list.
# Does runs_per_goal evolutions of each goal.
# For each goal, the output is the goal, the average number of generations, and the average number of active nodes,
#    and the number of succesful runs.  A run is successful if the goal is found in less than max_gens.
# Creates a CSV file with the parameter settings followed by the above described output for each run.
# Also writes this information to stdout so that the user can see the progress made.
# TO DO:  Figure out how to display the information to stdout as the process runs.
function evolve_goals( outstream::IOStream, summary::AbstractString, p, runs_per_goal, max_gens, number_random_goals, rseed)
    num_goals = 0
    goal_list = Void
    num_goals,goal_list = goal_list_setup(p, runs_per_goal, max_gens, number_random_goals )
    println(outstream,summary)
    println(outstream,Dates.now())
    print(outstream,  "host:          ",readall(`hostname`))
    println(outstream,"num processes: ",length(procs()))
    print(outstream,readall(`julia -v`))
    println(outstream,"num inputs:    ",p.numinputs)
    println(outstream,"num outputs:   ",p.numoutputs)
    println(outstream,"num per level: ",p.numperlevel)
    println(outstream,"num levels:    ",p.numlevels)
    println(outstream,"levels back:   ",p.numlevelsback)
    println(outstream,"num goals:     ",num_goals)
    println(outstream,"mu:            ",p.mu)
    println(outstream,"lambda:        ",p.lambda)
    println(outstream,"mutrate:       ",p.mutrate)
    println(outstream,"rand seed:     ",rseed)
    println(outstream,"fitfunc:       ",p.fitfunc)
    # reshape in the following line converts to a row vector to get rid of commas for CSV output
    println(outstream,"funcs: ",reshape([f.name for f in p.funcs],(1,length(p.funcs))))
    println(outstream,"max gens:      ",max_gens),
    println(outstream,"runs_per_goal: ",runs_per_goal),
    srand(rseed)
    sum_gens = 0
    if length(procs()) > 1
        r = pmap(p_mu_lambda,goal_list)
    else
        r = map(p_mu_lambda,goal_list)
    end
    average_array = averages(num_goals,runs_per_goal,max_gens,r)
    @printf("goal,  gens, active,  successful_runs\n")
    println(outstream,"goal, ave_gens, ave_active, succesful_runs") # header line for CSV
    for i in 1:num_goals
        sum_gens += average_array[i,1]
        @printf(outstream,"%#x, %6.1f, %6.1f, %4.0f\n",goal_list[i*runs_per_goal].truth_table[1],average_array[i,1],average_array[i,2],average_array[i,3])
        @printf("%#x, %6.1f, %6.1f, %4.0f\n",goal_list[i*runs_per_goal].truth_table[1],average_array[i,1],average_array[i,2],average_array[i,3])
    end
    average_gens = convert(AbstractFloat,sum_gens)/num_goals
    println(outstream,"average gens: ",average_gens)
    println("average gens: ",average_gens)
    outstream
end


text:
using OrdinaryDiffEq
using Trixi

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

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

initial_condition   = initial_condition_shock_bubble

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

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

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

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

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

summary_callback    = SummaryCallback()

analysis_interval   = 100
analysis_callback   = AnalysisCallback(semi, interval=analysis_interval)

alive_callback      = AliveCallback(analysis_interval=analysis_interval)

save_solution       = SaveSolutionCallback(interval=100,
                                           save_initial_solution=true,
                                           save_final_solution=true,
                                           solution_variables=:primitive)

amr_indicator = IndicatorHennemannGassner(semi,
                                          alpha_max=1.0,
                                          alpha_min=0.0001,
                                          alpha_smooth=true,
                                          variable=density_pressure)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                          base_level=5,
                                          max_level =9, max_threshold=0.01)

amr_callback = AMRCallback(semi, amr_controller,
                               interval=5,
                               adapt_initial_condition=true,
                               adapt_initial_condition_only_refine=true)

stepsize_callback   = StepsizeCallback(cfl=0.3)

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

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


text:
using HDF5
using SARA: TemperatureProfile
savepath = "SARA/NatCom2020/outer/data/"
xrd_file = "Bi2O3_19F44_01_outer_xrd_gradients.h5"
# xrd_file = "Bi2O3_19F44_01_outer_xrd_gradients_no_input_noise.h5"

# loading xrd data
xrd_f = h5open(savepath * xrd_file, "r")
xrd_temperatures = read(xrd_f, "temperatures")
xrd_dwelltimes = read(xrd_f, "dwelltimes")
xrd_gradients = read(xrd_f, "gradients")
xrd_var_gradients = read(xrd_f, "var(gradients)")
xrd_coefficients = read(xrd_f, "coefficients")
close(xrd_f)

xrd_f = h5open(savepath * "Bi2O3_19F44_01_outer_xrd_data.h5", "r")
xrd_positions = read(xrd_f, "positions")
close(xrd_f)

####### loading optical data
optical_file = "Bi2O3_19F44_01_outer_optical_gradients.h5"
optical_f = h5open(savepath * optical_file, "r")
optical_temperatures = read(optical_f, "temperatures")
optical_dwelltimes = read(optical_f, "dwelltimes")
optical_gradients = read(optical_f, "gradients")
optical_var_gradients = read(optical_f, "var(gradients)")
close(optical_f)

include("../inner/inner_load_data.jl")
path = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/Bi2O3_19F44_01/"
# path = ".."
file = "Bi2O3_19F44_01_inner_loop_data.json"
f = JSON.parsefile(path * file)
optical_positions, optical_coefficients, optical_rescaling_parameters,
        optical_inner_temperatures, optical_inner_dwelltimes = load_data(f)

using Plots
plotly()
dT = 50 # uncertainty in peak temperature in C
dx = 100/1000 # uncertainty in position in mm
P = TemperatureProfile(dT, dx)
function get_match(T_max, log10_τ, temperatures, dwelltimes)
   j = 1
   while true
      if T_max == temperatures[j] && log10_τ == dwelltimes[j]
          break
      end
      j += 1
   end
   return j
end

using SARA: iSARA, get_temperature_process
using GaussianDistributions
for i in 3:3
    T_max = xrd_temperatures[end, i] + 10 # adding 10 because of offset from T peak for outer temperature
    log10_τ = xrd_dwelltimes[end, i]
    j = get_match(T_max, log10_τ, optical_temperatures[end, :].+10, optical_dwelltimes[end, :])
    println("T_max = $T_max, $(optical_temperatures[end, j]+10), $(optical_inner_temperatures[j])")
    println("log10_τ = $log10_τ, $(optical_dwelltimes[end, j]), $(optical_inner_dwelltimes[j])")

    σ = .01
    l = .1
    θ = optical_rescaling_parameters[j]
    k = iSARA(l, θ)

    # coefficients as a function of position
    ncoeff = 3
    coeff_ind = 1:3 #[1, 2, 3, 4] #8, 16] # coefficient indices
    p2 = plot(layout = grid(1, ncoeff)) #@layout [a b c d]) # xlabel = "position (mm)", title = "coefficients"

    # GP for xrd data
    xrd_coeff = xrd_coefficients[coeff_ind, :, i]
    xrd_pos = xrd_positions[:, i]
    c = [conditional(Gaussian(k), xrd_pos, xrd_coeff[ci, :], σ^2) for ci in coeff_ind]
    cx = [mean(ci)(xi) for ci in c, xi in xrd_pos]
    plot!(xrd_pos, xrd_coeff', label = "xrd")
    plot!(xrd_pos, cx', label = "xrd gp")

    # GP for optical data
    opt_coeff = [optical_coefficients[j][ci][cj] for ci in coeff_ind, cj in 1:151]
    opt_pos = optical_positions[j]
    # cent, fom = center_of_symmetry(opt_coeff)
    # @. opt_pos -= opt_pos[cent]
    c = [conditional(Gaussian(k), opt_pos, opt_coeff[ci, :], σ^2) for ci in coeff_ind]
    cx = [mean(ci)(xi) for ci in c, xi in opt_pos]
    plot!(opt_pos, opt_coeff', label = "optical")
    plot!(opt_pos, cx', label = "opt gp")

    # coefficients as a function of temperature
    T_profile = P(T_max, log10_τ)

    xrd_temp = T_profile.(xrd_positions[:, i])
    opt_temp = T_profile.(optical_positions[j])
    p3 = plot(layout = grid(1, ncoeff)) # xlabel = "position (mm)", title = "coefficients"

    # xrd
    input_noise = Val(true)
    c = [get_temperature_process(Gaussian(k), xrd_pos, xrd_coeff[ci, :], σ, P, T_max, log10_τ, input_noise) for ci in coeff_ind]
    cx = [mean(ci)(xi) for ci in c, xi in xrd_temp]
    plot!(xrd_temp, xrd_coeff', label = "xrd")
    plot!(xrd_temp, cx', label = "xrd gp")

    # optical
    c = [get_temperature_process(Gaussian(k), opt_pos, opt_coeff[ci, :], σ, P, T_max, log10_τ, input_noise) for ci in coeff_ind]
    cx = [mean(ci)(xi) for ci in c, xi in opt_temp]
    plot!(opt_temp, opt_coeff', label = "optical")
    plot!(opt_temp, cx', label = "opt gp")

    p1 = plot(xlabel = "temperature (C)", title = "GP gradients", yaxis = :log10)
    plot!(xrd_temperatures[:, i], xrd_gradients[:, i],
            #ribbon = 2sqrt.(xrd_var_gradients[:, i]),
            label = "xrd")
    plot!(optical_temperatures[:, j], optical_gradients[:, j],
            #ribbon = 2sqrt.(optical_var_gradients[:, j]),
            label = "optical")

    local l = @layout [a; b; c]
    plot(p2, p3, p1, layout = l, legend = false) #, title = "T_max = $T_max, log10_τ = $log10_τ")
    gui()
end


text:
# Command: Is Element Enabled
"""
    isenabled(element::Element)::Bool

Is Element Enabled.
"""
function isenabled(element::Element)::Bool
    @unpack addr, id = element.session
    element_id = element.id
    response = HTTP.get(
        "$addr/session/$id/element/$element_id/enabled",
        [("Content-Type" => "application/json; charset=utf-8")],
    )
    @assert response.status == 200
    JSON3.read(response.body).value
end


text:
using JuMP

using LinearAlgebra

function parse_box_qp(filename::String)
    open(filename) do fp
        n = parse(Int, readline(fp))
        c = parse.(Float64, split(readline(fp)))
        @assert n == length(c)
        Q = zeros(Float64, n, n)
        for i = 1:n
            Q[:, i] = parse.(Float64, split(readline(fp)))
        end
        return Q, c
    end
end

function formulate_box_qp(Q::Matrix, c::Vector; epigraph::Bool = true)
    n = length(c)
    @assert n == size(Q, 1) == size(Q, 2)
    model = JuMP.Model()
    JuMP.@variable(model, 0 ≤ x[1:n] ≤ 1)
    if epigraph
        JuMP.@variable(model, y)
        JuMP.@constraint(model, y ≥ -0.5dot(x, Q * x) - dot(c, x))
        JuMP.@objective(model, Min, y)
    else
        JuMP.@objective(model, Min, -0.5dot(x, Q * x) - dot(c, x))
    end
    return JuMP.backend(model)
end


text:
"""
    columns_widths

Get the width of each column in a `Table`.

If user passed a `columns_widths` argument use that, otherwise
compute the width of each column based on the size of its contents
including the header and footer.
"""
function calc_columns_widths(
    N_cols::Int,
    N_rows::Int,
    columns_widths::Union{Nothing,Vector,Int},
    show_header::Bool,
    header,
    tb,
    sch,
    footer,
    hpad,
)
    if !isnothing(columns_widths)
        columns_widths = expand(columns_widths, N_cols)
        return columns_widths
    end

    headers_widths = show_header ? collect(width.(header)) : zeros(N_cols)
    data_widths = collect(map(c -> max(width.(tb[c])...), sch.names))
    footers_widths = isnothing(footer) ? zeros(N_cols) : collect(width.(footer))
    widths = hcat(headers_widths, data_widths, footers_widths)
    hpad = isa(hpad, Int) ? repeat([hpad], N_rows) : hpad
    widths = Int.([mapslices(x -> max(x...), widths; dims = 2)...] .+ hpad * 2)
    return widths
end

"""
    rows_heights

Get the height of each row in a `Table`
"""
function rows_heights(N_rows::Int, show_header::Bool, header, rows_values, footer, vpad)
    headers_height = show_header ? max(height.(header)...) : 0
    data_heights = collect(map(r -> max(height.(r)...), rows_values))
    footers_height = isnothing(footer) ? 0 : max(height.(footer)...)
    vpad = isa(vpad, Int) ? repeat([vpad], N_rows) : vpad
    heights = [headers_height, data_heights..., footers_height] .+ vpad .* 2
    return heights
end

"""
    expand

Expand single `Table` arguments into a vector if necessary.
"""
expand(v::Vector, N::Int) = v
expand(v::Union{Symbol,String,Int}, N::Int) = repeat([v], N)

"""
    assert_table_arguments

Check that arguments passed to `Table` match the Table's shape.
Single arguments are fine, but when a vector is passed it should have the appropriate size.
"""
function assert_table_arguments(
    N_cols,
    N_rows,
    show_header,
    header,
    header_style,
    header_justify,
    columns_style,
    columns_justify,
    columns_widths,
    footer,
    footer_style,
    footer_justify,
    hpad,
    vpad,
)
    # check header
    problems = []
    if !isnothing(header) && show_header
        length(header) != N_cols &&
            push!(problems, "Got header with length $(length(header)), expected $N_cols")
        (!isa(header_style, String) && length(header_style) != N_cols) && push!(
            problems,
            "Got header_style with length $(length(header_style)), expected $N_cols",
        )
        (!isa(header_justify, Symbol) && length(header_justify) != N_cols) && push!(
            problems,
            "Got header_justify with length $(length(header_justify)), expected $N_cols",
        )
    end

    # check columns
    (!isa(columns_style, String) && length(columns_style) != N_cols) && push!(
        problems,
        "Got columns_style with length $(length(columns_style)), expected $N_cols",
    )
    (!isa(columns_justify, Symbol) && length(columns_justify) != N_cols) && push!(
        problems,
        "Got columns_justify with length $(length(columns_justify)), expected $N_cols",
    )
    columns_widths isa Vector &&
        length(columns_widths) != N_cols &&
        push!(
            problems,
            "Got columns_widths with length $(length(columns_widths)), expected $N_cols",
        )

    # check footer
    if !isnothing(footer)
        isa(footer, Function) ||
            (length(footer)) != N_cols && push!(
                problems,
                "Got footer with length $(length(footer)), expected $N_cols",
            )
        (!isa(footer_style, String) && length(footer_style) != N_cols) && push!(
            problems,
            "Got footer_style with length $(length(footer_style)), expected $N_cols",
        )
        (!isa(footer_justify, Symbol) && length(footer_justify) != N_cols) && push!(
            problems,
            "Got footer_justify with length $(length(footer_justify)), expected $N_cols",
        )
    end

    # check hpad
    (!isa(hpad, Int) && length(hpad) != N_cols) &&
        push!(problems, "Got hpad with length $(length(hpad)), expected $N_cols")
    (!isa(vpad, Int) && length(vpad) != N_rows) &&
        push!(problems, "Got vpad with length $(length(vpad)), expected $N_rows")

    # if there were problems, alert user and fail gracefully
    if length(problems) > 0
        @warn "Failed to create Term.Table"
        warn_color = orange
        tprintln.("  {$warn_color}" .* problems .* "{/$warn_color}"; highlight = true)
    end
    return length(problems) == 0
end


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

mutable struct RollingUpgradeProgressInfo <: SwaggerModel
    successfulInstanceCount::Any # spec type: Union{ Nothing, Int32 } # spec name: successfulInstanceCount
    failedInstanceCount::Any # spec type: Union{ Nothing, Int32 } # spec name: failedInstanceCount
    inProgressInstanceCount::Any # spec type: Union{ Nothing, Int32 } # spec name: inProgressInstanceCount
    pendingInstanceCount::Any # spec type: Union{ Nothing, Int32 } # spec name: pendingInstanceCount

    function RollingUpgradeProgressInfo(;successfulInstanceCount=nothing, failedInstanceCount=nothing, inProgressInstanceCount=nothing, pendingInstanceCount=nothing)
        o = new()
        validate_property(RollingUpgradeProgressInfo, Symbol("successfulInstanceCount"), successfulInstanceCount)
        setfield!(o, Symbol("successfulInstanceCount"), successfulInstanceCount)
        validate_property(RollingUpgradeProgressInfo, Symbol("failedInstanceCount"), failedInstanceCount)
        setfield!(o, Symbol("failedInstanceCount"), failedInstanceCount)
        validate_property(RollingUpgradeProgressInfo, Symbol("inProgressInstanceCount"), inProgressInstanceCount)
        setfield!(o, Symbol("inProgressInstanceCount"), inProgressInstanceCount)
        validate_property(RollingUpgradeProgressInfo, Symbol("pendingInstanceCount"), pendingInstanceCount)
        setfield!(o, Symbol("pendingInstanceCount"), pendingInstanceCount)
        o
    end
end # type RollingUpgradeProgressInfo

const _property_map_RollingUpgradeProgressInfo = Dict{Symbol,Symbol}(Symbol("successfulInstanceCount")=>Symbol("successfulInstanceCount"), Symbol("failedInstanceCount")=>Symbol("failedInstanceCount"), Symbol("inProgressInstanceCount")=>Symbol("inProgressInstanceCount"), Symbol("pendingInstanceCount")=>Symbol("pendingInstanceCount"))
const _property_types_RollingUpgradeProgressInfo = Dict{Symbol,String}(Symbol("successfulInstanceCount")=>"Int32", Symbol("failedInstanceCount")=>"Int32", Symbol("inProgressInstanceCount")=>"Int32", Symbol("pendingInstanceCount")=>"Int32")
Base.propertynames(::Type{ RollingUpgradeProgressInfo }) = collect(keys(_property_map_RollingUpgradeProgressInfo))
Swagger.property_type(::Type{ RollingUpgradeProgressInfo }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_RollingUpgradeProgressInfo[name]))}
Swagger.field_name(::Type{ RollingUpgradeProgressInfo }, property_name::Symbol) =  _property_map_RollingUpgradeProgressInfo[property_name]

function check_required(o::RollingUpgradeProgressInfo)
    true
end

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


text:
export BroadCast

"""
`BroadCast(A::AbstractOperator, dim_out...)`

BroadCast the codomain dimensions of an `AbstractOperator`.

```julia
julia> A = Eye(2)
I  ℝ^2 -> ℝ^2

julia> B = BroadCast(A,(2,3))
.I  ℝ^2 -> ℝ^(2, 3)

julia> B*[1.;2.]
2×3 Array{Float64,2}:
 1.0  1.0  1.0
 2.0  2.0  2.0

```

"""
struct BroadCast{N,
		 L <: AbstractOperator,
		 T <: AbstractArray,
		 D <: AbstractArray,
		 M,
		 C <: NTuple{M,Colon},
		 I <: CartesianIndices
		 } <: AbstractOperator
	A::L
	dim_out::NTuple{N,Int}
	bufC::T
	bufD::D
	cols::C
	idxs::I

	function BroadCast(A::L,dim_out::NTuple{N,Int},bufC::T, bufD::D) where {N,
								      L<:AbstractOperator,
								      T<:AbstractArray,
								      D<:AbstractArray
								      }
		Base.Broadcast.check_broadcast_shape(dim_out,size(A,1))
		if size(A,1) != (1,)
			M = length(size(A,1))
			cols = ([Colon() for i = 1:M]...,)
			idxs = CartesianIndices((dim_out[M+1:end]...,))
			new{N,L,T,D,M,typeof(cols),typeof(idxs)}(A,dim_out,bufC,bufD,cols,idxs)
		else #singleton case
			M = 0
			idxs = CartesianIndices((1,))
			new{N,L,T,D,M,NTuple{0,Colon},typeof(idxs)}(A,dim_out,bufC,bufD,(),idxs)
		end

	end
end

# Constructors

BroadCast(A::L, dim_out::NTuple{N,Int}) where {N,L<:AbstractOperator} =
BroadCast(A, dim_out, zeros(codomainType(A),size(A,1)), zeros(domainType(A),size(A,2)) )

# Mappings

function mul!(y::CC, R::BroadCast{N,L,T,D,M}, b::DD) where {N,L,T,D,M,CC,DD}
	mul!(R.bufC, R.A, b)
	y .= R.bufC
end

function mul!(y::CC, A::AdjointOperator{BroadCast{N,L,T,D,M,C,I}}, b::DD) where {N,L,T,D,M,C,I,CC,DD}
    R = A.A
	fill!(y, 0.)
	for i in R.idxs
		@views mul!(R.bufD, R.A', b[R.cols...,i.I...])
		y .+= R.bufD
	end
	return y
end

#singleton
function mul!(y::CC, A::AdjointOperator{BroadCast{N,L,T,D,0,C,I}}, b::DD) where {N,L,T,D,C,I,CC,DD}
    R = A.A
	fill!(y, 0.)
	bii = zeros(eltype(b),1)
	for bi in b
		bii[1] = bi
		mul!(R.bufD, R.A', bii)
		y .+= R.bufD
	end
	return y
end

#TODO make this more general
#length(dim_out) == size(A,1) e.g. a .= b; size(a) = (m,n) size(b) = (1,n) matrix out, column in
function mul!(y::CC, A::AdjointOperator{BroadCast{2,L,T,D,2,C,I}}, b::DD) where {L,T,D,C,I,CC,DD}
    R = A.A
	fill!(y, 0.)
	for i = 1:size(b,1)
		@views mul!(R.bufD, R.A, b[i,:]')
		y .+= R.bufD
	end
	return y
end

#singleton Eye
function mul!(y::CC, R::AdjointOperator{BroadCast{N,L,T,D,0,C,I}}, b::DD) where {N,L<:Eye,T,D,C,I,CC,DD}
	sum!(y,b)
end

# Properties

size(R::BroadCast) = (R.dim_out, size(R.A,2))

  domainType(  R::BroadCast) =   domainType(R.A)
codomainType(  R::BroadCast) = codomainType(R.A)

is_linear(      R::BroadCast) = is_linear(R.A)
is_null(        R::BroadCast) = is_null(R.A)

fun_name(R::BroadCast) = "."*fun_name(R.A)
remove_displacement(B::BroadCast) = BroadCast(remove_displacement(B.A), B.dim_out, B.bufC, B.bufD)


text:
"""
    codas(decisionMat, weights, fs)
Apply CODAS (COmbinative Distance-based ASsessment) method for a given matrix, weights and, type of criteria.
# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alternatives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fs::Array{Function,1}`: m-vector of type of criteria. The benefit criteria shown with "maximum", and the cost criteria shown with "minimum".
 - `tau::Float64`: tau parameter for the algorithm. The default is 0.02.
# Description
codas() applies the CODAS method to rank n alternatives subject to m criteria and criteria type vector.
# Output
- `::CODASResult`: CODASResult object that holds multiple outputs including scores and best index.
# Examples
```julia-repl
julia> decmat
7×5 Array{Float64,2}:
 60.0   0.4   2540.0   500.0   990.0
  6.35  0.15  1016.0  3000.0  1041.0
  6.8   0.1   1727.2  1500.0  1676.0
 10.0   0.2   1000.0  2000.0   965.0
  2.5   0.1    560.0   500.0   915.0
  4.5   0.08  1016.0   350.0   508.0
  3.0   0.1   1778.0  1000.0   920.0

julia> df = DataFrame(decmat)
7×5 DataFrame
 Row │ x1       x2       x3       x4       x5
     │ Float64  Float64  Float64  Float64  Float64
─────┼─────────────────────────────────────────────
   1 │   60.0      0.4    2540.0    500.0    990.0
   2 │    6.35     0.15   1016.0   3000.0   1041.0
   3 │    6.8      0.1    1727.2   1500.0   1676.0
   4 │   10.0      0.2    1000.0   2000.0    965.0
   5 │    2.5      0.1     560.0    500.0    915.0
   6 │    4.5      0.08   1016.0    350.0    508.0
   7 │    3.0      0.1    1778.0   1000.0    920.0
julia> result = codas(df, w, fs);
julia> result.bestIndex
2
julia> result.scores
7-element Array{Float64,1}:
  0.5121764914884954
  1.463300034504913
  1.0715325899642418
 -0.21246799780012637
 -1.8515205523193041
 -1.1716767695713806
  0.18865620373316055

```
# References
Keshavarz Ghorabaee, M., Zavadskas, E. K., Turskis, Z., & Antucheviciene, J. (2016). A new combinative distance-based assessment (CODAS) method for multi-criteria decision-making. Economic Computation & Economic Cybernetics Studies & Research, 50(3), 25-44.
"""
function codas(decisionMat::DataFrame, weight::Array{Float64,1}, fns::Array{Function,1}, tau::Float64=0.02)::CODASResult

    #mat = convert(Matrix, decisionMat)
    mat = Matrix(decisionMat)

    nrows, ncols = size(decisionMat)
    w = unitize(weight)

    colMax = colmaxs(decisionMat)
    colMin = colmins(decisionMat)

    A = similar(decisionMat)

    for i in 1:ncols
        if fns[i] == maximum
            @inbounds A[:, i] = decisionMat[:, i] ./ colMax[i]
        elseif fns[i] == minimum
            @inbounds A[:, i] = colMin[i] ./ decisionMat[:, i]
        end
    end

    wA = similar(A)

    for i in 1:ncols
        @inbounds wA[:, i] = A[:, i] .* w[i]
    end

    wAmin = colmins(wA)

    E = similar(A)
    for i in 1:nrows
        for j in 1:ncols
            E[i,j] = (wA[i,j] .- wAmin[j])^2
        end
    end

    Euc = zeros(Float64, nrows)
    for i in 1:nrows
        Euc[i] = sqrt(sum(E[i,:]))
    end

    T = similar(A)
    for i in 1:nrows
        for j in 1:ncols
            T[i,j] = abs(wA[i,j] .- wAmin[j])
        end
    end

    Tax = zeros(Float64, nrows)
    for i in 1:nrows
        Tax[i] = sum(T[i,:])
    end

    EA = zeros(Float64, nrows, nrows)
    TA = zeros(Float64, nrows, nrows)
    RA = zeros(Float64, nrows, nrows)

    for i in 1:nrows
        for j in 1:nrows
            EA[i,j] = Euc[i] - Euc[j]
            TA[i,j] = Tax[i] - Tax[j]
            if abs(Euc[i] - Euc[j]) < tau
                RA[i,j] = 0
            else
                RA[i,j] = 1
            end
        end
    end

    scores = zeros(Float64, nrows)
    for i in 1:nrows
        scores[i] = sum(EA[i,:] .+ (RA[i,:] .* TA[i,:]))
    end

    rankings = sortperm(scores)
    bestIndex = rankings |> last

    result = CODASResult(
        decisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
# Symmetric triangular distribution
immutable TriangularDist <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function TriangularDist(l::Real, s::Real)
        s > zero(s) || error("scale must be positive")
        new(float64(l), float64(s))
    end
end

TriangularDist(location::Real) = TriangularDist(location, 1.0)
TriangularDist() = TriangularDist(0.0, 1.0)

## Support
isupperbounded(::Union(TriangularDist, Type{TriangularDist})) = true
islowerbounded(::Union(TriangularDist, Type{TriangularDist})) = true
isbounded(::Union(TriangularDist, Type{TriangularDist})) = true

minimum(d::TriangularDist) = d.location - d.scale
maximum(d::TriangularDist) = d.location + d.scale
insupport(d::TriangularDist, x::Real) = minimum(d) <= x <= maximum(d)

## Properties
mean(d::TriangularDist) = d.location
median(d::TriangularDist) = d.location
mode(d::TriangularDist) = d.location

var(d::TriangularDist) = d.scale^2 / 6.0
skewness(d::TriangularDist) = 0.0
kurtosis(d::TriangularDist) = -0.6

entropy(d::TriangularDist) = 0.5 + log(d.scale)

## Functions
function pdf(d::TriangularDist, x::Real)
    if insupport(d, x)
        (1.0 - abs(x-d.location)/d.scale) / d.scale
    else
        0.0
    end
end

function logpdf(d::TriangularDist, x::Real)
    if insupport(d, x)
        log1p(-abs(x-d.location)/d.scale) -log(d.scale)
    else
        -Inf
    end
end

function cdf(d::TriangularDist, x::Real)
    u = (x - d.location) / d.scale
    if u <= -1.0
        0.0
    elseif u <= 0.0
        0.5*(1.0 + u)^2
    elseif u <= 1.0
        1.0 - 0.5*(1.0 - u)^2
    else
        1.0
    end
end
function ccdf(d::TriangularDist, x::Real)
    u = (x - d.location) / d.scale
    if u <= -1.0
        1.0
    elseif u <= 0.0
        1.0 - 0.5*(1.0 + u)^2
    elseif u <= 1.0
        0.5*(1.0 - u)^2
    else
        0.0
    end
end

function logcdf(d::TriangularDist, x::Real)
    u = (x - d.location) / d.scale
    if u <= -1.0
        -Inf
    elseif u <= 0.0
        loghalf + 2.0*log1p(u)
    elseif u <= 1.0
        log1p(-0.5*(1.0 - u)^2)
    else
        0.0
    end
end
function logccdf(d::TriangularDist, x::Real)
    u = (x - d.location) / d.scale
    if u <= -1.0
        0.0
    elseif u <= 0.0
        log1p(- 0.5*(1.0 + u)^2)
    elseif u <= 1.0
        loghalf + 2.0*log1p(-u)
    else
        -Inf
    end
end

function quantile(d::TriangularDist, p::Real)
    @checkquantile p begin
        if p < 0.5
            d.location - d.scale*(1.0 - sqrt(2.0*p))
        else
            d.location + d.scale*(1.0 - sqrt(2.0*(1.0-p)))
        end
    end
end
function cquantile(d::TriangularDist, p::Real)
    @checkquantile p begin
        if p > 0.5
            d.location - d.scale*(1.0 - sqrt(2.0*(1.0-p)))
        else
            d.location + d.scale*(1.0 - sqrt(2.0*p))
        end
    end
end
function invlogcdf(d::TriangularDist, lp::Real)
    @checkinvlogcdf lp begin
        if lp < loghalf
            d.location + d.scale*expm1(0.5*(lp-loghalf))
        else
            d.location + d.scale*(1.0 - sqrt(-2.0*expm1(lp)))
        end
    end
end
function invlogccdf(d::TriangularDist, lp::Real)
    @checkinvlogcdf lp begin
        if lp > loghalf
            d.location - d.scale*(1.0 - sqrt(-2.0*expm1(lp)))
        else
            d.location - d.scale*expm1(0.5*(lp-loghalf))
        end
    end
end

function mgf(d::TriangularDist, t::Real)
    a = d.scale*t
    a == zero(a) && return one(a)
    4.0*exp(d.location*t)*(sinh(0.5*a)/a)^2
end
function cf(d::TriangularDist, t::Real)
    a = d.scale*t
    a == zero(a) && return complex(one(a))
    4.0*exp(im*d.location*t)*(sin(0.5*a)/a)^2
end

## Sampling
function rand(d::TriangularDist)
    ξ1, ξ2 = rand(), rand()
    return d.location + (ξ1 - ξ2) * d.scale
end


text:
include("parameter.jl")

using Random
using StaticArrays

PID = length(ARGS) > 0 ? PID = parse(Int, ARGS[1]) : rand(100_000:999_999)
Random.seed!(PID) #initialize rand(), which is already global and by default using MersenneTwister
@assert PID >= 0 "RNG seed must be >=0"

@flush println("Seed/PID: $PID")

include("./parquetMC/oneBody/group.jl")
using .OneBody
include("./markov.jl")

state, groups, updates=OneBody.init(para, grids)
markov(state, groups, updates, 1000_000, Random.GLOBAL_RNG)

# include("markov.jl")
# Markov.init()

################# Test and Benchmark ############################
# Markov.test()

# import BenchmarkTools: @btime
# for _order = 1:Order
#     println("Benchmark Order $_order")
#     @btime Markov.eval(o) samples = 1 evals = 100 setup = (o = $_order)
# println(sum(Markov.ver4[_order].weight))
# end
# exit()

# println("Start Simulation ...")
# block = 0

# printTimer = StopWatch(PrintTime, Markov.printStatus)
# saveTimer = StopWatch(SaveTime, Markov.save)
# reweightTimer = StopWatch(ReWeightTime, Markov.reweight)
# messageTimer = StopWatch(MessageTime, Markov.save)

# for block = 1:TotalBlock
#     for i = 1:1000_000
#         Curr.step += 1
#         x = rand(Curr.rng)
#         if x < 1.0 / 6.0
#             Markov.increaseOrder()
#         elseif x < 2.0 / 6.0
#             Markov.decreaseOrder()
#         elseif x < 3.0 / 6.0
#             # Markov.changeK()
#         elseif x < 4.0 / 6.0
#             Markov.changeTau()
#         elseif x < 5.0 / 6.0
#             Markov.changeExtK()
#         else
#             Markov.changeExtTau()
#         end
#         # println(Curr.T[1], "-->", Curr.T[LastTidx])

#         i % 8 == 0 && Markov.measure()

#         if i % 1000 == 0
#             check(printTimer)
#             check(saveTimer)
#             check(reweightTimer)
#         end
#     end
# end

# Markov.printStatus()
# Markov.save()
# println("End Simulation. ")

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
#       Following 2 lines of code is an MVV/LVA scheme:
#	val = 128 * PIECEVALUES[capt.getCapt()] + PIECEVALUES[capt.getPiec()];
#	if (capt.isPromotion()) val += 512 * PIECEVALUES[capt.getProm()];

function MVVLVA(c_koma::Int, flag::Int, piece::Int)
    v::Int = 0
    # c_koma::Int = seeMoveCapt(m)
    # piece::Int = seeMovePiece(m) & 0x0f
    capt::Int = c_koma & 0x0f
    piece2::Int = piece & 0x0f

    matCapt::Int = (capt == MJNONE)?0:MATERIAL[capt]
    matPiece::Int = (piece2 == MJNONE)?0:MATERIAL[piece2]

    v = 128 * matCapt + matPiece
    if (flag & FLAG_NARI) == FLAG_NARI
        v += 512 * matPiece
    end

    return v
end

# 早速immutable構文を使ってみるなど(^_^;)

immutable type Move
    move::Uint
    value::Int
    Move(Piece::Int,From::Int,To::Int,Flags::Int,Capt::Int,value::Int) = new(((From&0x000000ff)<<24)|((To&0x000000ff)<<16)|((Capt&0x0000000f)<<12)|((Flags&0x00000007)<<8)|(Piece&0x0000001f),MVVLVA(Capt,Flags,Piece))
end

seeMoveFrom(m::Move)  = int((m.move & 0xff000000)>>>24)
seeMoveTo(m::Move)    = int((m.move & 0x00ff0000)>>>16)
seeMoveCapt(m::Move)  = int((m.move & 0x0000f000)>>>12)
seeMoveFlag(m::Move)  = int((m.move & 0x00000f00)>>>8)
seeMovePiece(m::Move) = int(m.move & 0x0000001f)

const FLAG_NARI = int(0x1) # 成ったときに可能手にフラグを立てる
const FLAG_TORI = int(0x2) # 駒を取ったときに可能手にフラグを立てる
const FLAG_UCHI = int(0x4) # 駒を打ったときに可能手にフラグを立てる

# あとは必要に応じて関数を追加する。先手か後手かなど

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/JapaneseChess

max_forks_repo_name: kimrin/JapaneseChess