text:
# Translation of David Hernandez's nbody.c for integrating hiercharical
# system with BH15 integrator.  Please cite Hernandez & Bertschinger (2015)
# if using this in a paper.

const YEAR  = 365.242
const GNEWT = 39.4845/YEAR^2
const NDIM  = 3
const third = 1./3.
include("kepler_step.jl")
include("init_nbody.jl")

function nbody!(t0::Float64,h::Float64,tmax::Float64,elements::Array{Float64,2},t_out::Array{Float64,1},state_out::Array{Float64,2},nstep::Int64)
fcons = open("fcons.txt","w");
n=8
m=zeros(n)
x=zeros(NDIM,n)
v=zeros(NDIM,n)
L0=zeros(NDIM)
p0=zeros(NDIM)
xcm0=zeros(NDIM)
ssave = zeros(n,n,3)
# Read in the initial orbital elements from a file (this
# should be moved to an input parameter):
elements = readdlm("elements.txt",',')
for i=1:n
  m[i] = elements[i,1]
end
# Initialize the N-body problem using nested hierarchy of Keplerians:
x,v = init_nbody(elements,t0,n)
#    infig1!(m,x,v,n)
# Compute the conserved quantities:
E0=consq!(m,x,v,n,L0,p0,xcm0)
#println("E0: ",E0," L0: ",L0," p0: ",p0," xcm0: ",xcm0)
#read(STDIN,Char)
# Set the time to the initial time:
t = t0
# Set step counter to zero:
i=0
#nstep = 1
# Loop over time steps:
while t < t0+tmax
  # Increment time by the time step:
  t += h
  # Increment counter by one:
  i +=1
  # Carry out a phi^2 mapping step:
  phi2!(x,v,h,m,n)
  # Every nstep steps, output the state of the system:
  if mod(i,nstep) == 0
    state_string = @sprintf("%18.12f",t)
    for j=1:n
      state_string = string(state_string,@sprintf("%18.12f",x[1,j]),@sprintf("%18.12f",x[3,j]),@sprintf("%18.12f",v[1,j]),@sprintf("%18.12f",v[3,j]))
    end
    state_string = state_string*"\n"
    print(fcons,state_string)
  end
  #
end
# Compute the final values of the conserved quantities:
L=zeros(NDIM)
p=zeros(NDIM)
xcm=zeros(NDIM)
E=consq!(m,x,v,n,L,p,xcm)
println("dE/E=",(E0-E)/E0)
# Close the output file:
close(fcons)
return
end

# Advances the center of mass of a binary
function centerm!(m::Array{Float64,1},x::Array{Float64,2},v::Array{Float64,2},delx::Array{Float64,1},delv::Array{Float64,1},i::Int64,j::Int64,h::Float64)
vcm=zeros(NDIM)
mij =m[i] + m[j]
if mij == 0
  for k=1:NDIM
    vcm[k] = (v[k,i]+v[k,j])/2
    x[k,i] += h*vcm[k]
    x[k,j] += h*vcm[k]
  end
else
  for k=1:NDIM
    vcm[k] = (m[i]*v[k,i] + m[j]*v[k,j])/mij
    x[k,i] +=  m[j]/mij*delx[k] + h*vcm[k]
    x[k,j] += -m[i]/mij*delx[k] + h*vcm[k]
    v[k,i] +=  m[j]/mij*delv[k]
    v[k,j] += -m[i]/mij*delv[k]
  end
end
return
end

# Drifts bodies i & j
function driftij!(x::Array{Float64,2},v::Array{Float64,2},i::Int64,j::Int64,h::Float64)
for k=1:NDIM
  x[k,i] += h*v[k,i]
  x[k,j] += h*v[k,j]
end
return
end

# Carries out a Kepler step for bodies i & j
function keplerij!(m::Array{Float64,1},x::Array{Float64,2},v::Array{Float64,2},i::Int64,j::Int64,h::Float64)
# The state vector has: 1 time; 2-4 position; 5-7 velocity; 8 r0; 9 dr0dt; 10 beta; 11 s; 12 ds
# Initial state:
s0 = zeros(Float64,12)
# Final state (after a step):
s = zeros(Float64,12)
delx = zeros(NDIM)
delv = zeros(NDIM)
for k=1:NDIM
  s0[1+k     ] = x[k,i] - x[k,j]
  s0[1+k+NDIM] = v[k,i] - v[k,j]
end
gm = GNEWT*(m[i]+m[j])
if gm == 0
  for k=1:NDIM
    x[k,i] += h*v[k,i]
    x[k,j] += h*v[k,j]
  end
else
  kepler_step!(gm, h, s0, s)
  for k=1:NDIM
    delx[k] = s[1+k] - s0[1+k]
    delv[k] = s[1+NDIM+k] - s0[1+NDIM+k]
  end
# Advance center of mass:
  centerm!(m,x,v,delx,delv,i,j,h)
end
return
end

# Drifts all particles:
function drift!(x::Array{Float64,2},v::Array{Float64,2},h::Float64,n::Int64)
for i=1:n
  for j=1:NDIM
    x[j,i] += h*v[j,i]
  end
end
return
end

# Carries out the phi^2 mapping
function phi2!(x::Array{Float64,2},v::Array{Float64,2},h::Float64,m::Array{Float64,1},n::Int64)
drift!(x,v,h/2,n)
for i=1:n-1
  for j=i+1:n
    driftij!(x,v,i,j,-h/2)
    keplerij!(m,x,v,i,j,h/2)
  end
end
for i=n-1:-1:1
  for j=n:-1:i+1
    keplerij!(m,x,v,i,j,h/2)
    driftij!(x,v,i,j,-h/2)
  end
end
drift!(x,v,h/2,n)
return
end

# Initializes the system for a specific dynamical problem: eccentric binary orbiting
# a larger body.  Units in ~AU & yr.
function infig1!(m::Array{Float64,1},x::Array{Float64,2},v::Array{Float64,2},n::Int64)
a0 = 0.0125
a1 = 1.0
e0 = 0.6
e1 = 0
m[1] = 1e-3
m[2] = 1e-3
m[3] = 1.0
mu0 = (m[1]*m[2])/(m[1]+m[2])
mt0 = m[1]+m[2]
mu1 = (m[3]*mt0)/(m[3]+mt0)
mt1 = m[3]+mt0
x0 = a0*(1+e0)
x1 = a1*(1+e1)
v0 = sqrt(GNEWT*mt0/a0*(1-e0)/(1+e0))
v1 = sqrt(GNEWT*mt1/a1*(1-e1)/(1+e1))
fa0 = m[1]/mt0
fb0 = m[2]/mt0
fa1 = mt0/mt1
fb1 = m[3]/mt1
x[1,1] = -fb1*x1-fb0*x0
x[1,2] = -fb1*x1+fa0*x0
x[1,3] = fa1*x1
v[2,1] = -fb1*v1-fb0*v0
v[2,2] = -fb1*v1+fa0*v0
v[2,3] = fa1*v1
return
end

# Computes conserved quantities:
function consq!(m::Array{Float64,1},x::Array{Float64,2},v::Array{Float64,2},n::Int64,
              L::Array{Float64,1},p::Array{Float64,1},xcm::Array{Float64,1})
fill!(p,0.0) # =zeros(NDIM)
fill!(xcm,0.0) # =zeros(NDIM)
fill!(L,0.0) #  =zeros(NDIM)
E = 0.0
for i=1:n
  if m[i] != 0
    for j=1:NDIM
      p[j] += m[i]*v[j,i];
      xcm[j] += m[i]*x[j,i];
    end
    L[1] += m[i]*(x[2,i]*v[3,i]-x[3,i]*v[2,i]);
    L[2] += m[i]*(x[3,i]*v[1,i]-x[1,i]*v[3,i]);
    L[3] += m[i]*(x[1,i]*v[2,i]-x[2,i]*v[1,i]);
  end
  for j=1:NDIM
    E += 1.0/2.0*m[i]*v[j,i]*v[j,i];
  end
  for j=i+1:n
    if m[j] != 0
      rij = 0.;
      for k=1:NDIM
        rij += (x[k,i]-x[k,j])*(x[k,i]-x[k,j]);
      end
      rij = sqrt(rij);
      E -= GNEWT*m[i]*m[j]/rij;
    end
  end
end
mt=sum(m)
if mt != 0
  for i=1:NDIM
    xcm[i] /= mt;
    p[i] /= mt;
  end
end
println("L: ",L)
println("p: ",p)
println("xcm: ",xcm)
return E
end


text:
s1 = "The quick brown fox jumps over the lazy dog α,β,γ";
println(s1);
print("this is ");
print("one line\n");

c1 = 'c';
println(c1);
println("this is c1 -> ",c1); # ok
# println("this is c1 -> "+c1); # not ok

num = 7;
println(7);
println("num -> ", num);

s1_caps = uppercase(s1)
s1_lower = lowercase(s1)
println(s1_caps, "\n", s1_lower)

println(s1[11]);
println(s1[1:11]);
show(s1[end-10:end]);
println();

s3 = string("this", " and", " that")
println(s3)

ext: jl

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# integer program formulation by Oddo Zhang, for move-up
# From thesis: "Simulation optimisation and Markov models for dynamic ambulance redeployment"

# initialise data relevant to move up
function initZhangIp!(sim::Simulation, zhangIpData::ZhangIpData)
	zid = sim.moveUpData.zhangIpData = zhangIpData # shorthand

	# check number of stations
	numStations = sim.numStations
	@assert(length(zid.marginalBenefits) == numStations)
	@assert(length(zid.stationCapacities) == numStations)

	# check station capacities
	for i = 1:numStations
		@assert(zid.stationCapacities[i] <= sim.stations[i].capacity)
	end

	zid.stationSlots = []
	zid.benefitSlots = []
	for i = 1:numStations
		for j = 1:zid.stationCapacities[i]
			push!(zid.stationSlots, i)
			push!(zid.benefitSlots, zid.marginalBenefits[i][j])
		end
	end

	zid.marginalBenefitsDecreasing = all(i -> issorted(zid.marginalBenefits[i], lt=<=, rev=true), 1:numStations)
	if !zid.marginalBenefitsDecreasing
		p1 = vcat([findall(zid.stationSlots .== j)[1:end-1] for j = 1:numStations]...)
		p2 = vcat([findall(zid.stationSlots .== j)[2:end] for j = 1:numStations]...)
		zid.stationSlotsOrderPairs = hcat(p1, p2)
	end

	@assert(sim.travel.numSets == 1) # otherwise, need to be careful in calculating the regret travel time for move up of on-road ambulances
end
function initZhangIp!(sim::Simulation; paramsFilename::String = "")
	zid = readZhangIpParamsFile(paramsFilename)
	initZhangIp!(sim, zid)
end

# determine move ups to make at current time
# returns list of ambulances to be moved, and list of their destinations (stations)
function zhangIpMoveUp(sim::Simulation)::Tuple{Vector{Ambulance}, Vector{Station}}

	# shorthand names:
	zid = sim.moveUpData.zhangIpData
	@unpack stationSlots, benefitSlots = zid
	@unpack ambulances, stations, numStations = sim

	# get currently movable ambulances, and at-hospital ambulances
	movableAmbs = filter(a -> isAmbMovable(a), ambulances)
	atHospitalAmbs = filter(a -> a.status == ambAtHospital, ambulances)
	@assert(!any(a -> isAmbMovable(a) && a.status == ambAtHospital, ambulances)) # @assert(isempty(intersect(movableAmbs, atHospitalAmbs))) is much slower
	if isempty(movableAmbs) return moveUpNull() end

	# let moveUpAmbs be the ambulances that can be moved now, and those that can be moved later (currently at-hospital)
	moveUpAmbs = vcat(movableAmbs, atHospitalAmbs)
	numMoveUpAmbs = length(moveUpAmbs)

	# calculate travel time for each move-up ambulance to reach every station
	ambToStationTimes = Array{Float,2}(undef, numMoveUpAmbs, numStations)
	for i = 1:numMoveUpAmbs
		ambToStationTimes[i,:] = ambMoveUpTravelTimes!(sim, moveUpAmbs[i])
	end

	# calculate adjustedAmbToStationTimes from ambToStationTimes, according to Zhang's thesis
	# for at-station or newly-freed ambulances: adjusted time = original time
	# for on-road ambulances: adjusted time = original time * ( time spent on route + time to drive to station - time from route origin to station <= regret-travel-time threshold ? discount factor : 1 )
	# for at-hospital ambulances: adjusted time = original time + expected remaining handover duration
	adjustedAmbToStationTimes = copy(ambToStationTimes)
	travelMode = getTravelMode!(sim.travel, lowPriority, sim.time)
	for i = 1:numMoveUpAmbs
		ambulance = moveUpAmbs[i]
		if ambulance.status == ambIdleAtStation || ambulance.status == ambFreeAfterCall
			# no change, adjustedAmbToStationTimes[i,:] = ambToStationTimes[i,:]
		elseif isGoingToStation(ambulance.status)
			for j = 1:numStations
				# calculate "regret" travel time
				if ambulance.stationIndex == j
					regretTravelTime = 0.0
				else
					# get travel time from ambulance.route.startLoc to station j
					time1 = ambulance.route.startFNodeTime - ambulance.route.startTime
					pathTravelTime = shortestPathTravelTime(sim.net, travelMode.index, ambulance.route.startFNode, stations[j].nearestNodeIndex)
					time2 = offRoadTravelTime(travelMode, stations[j].nearestNodeDist)
					tob = time1 + pathTravelTime + time2 # on and off road travel times

					tox = sim.time - ambulance.route.startTime # time spent on current route
					txb = ambToStationTimes[i,j]
					regretTravelTime = tox + txb - tob
				end

				# apply discount if regret travel time is within threshold
				if regretTravelTime <= zid.regretTravelTimeThreshold
					adjustedAmbToStationTimes[i,j] *= zid.onRoadMoveUpDiscountFactor
				end
			end
		elseif ambulance.status == ambAtHospital
			adjustedAmbToStationTimes[i,:] .+= zid.expectedHospitalHandoverDuration
		else
			error()
		end
	end

	ambToStationCosts = adjustedAmbToStationTimes * zid.travelTimeCost

	# solve
	if zid.marginalBenefitsDecreasing
		# solve as assignment problem
		(movableAmbStations, atHospitalAmbStations) = solveZhangIpAssignmentProblem(stationSlots, benefitSlots, ambToStationCosts, length(movableAmbs), length(atHospitalAmbs))
	else
		# solve as an integer program
		(movableAmbStations, atHospitalAmbStations) = solveZhangIp(stationSlots, benefitSlots, ambToStationCosts, length(movableAmbs), length(atHospitalAmbs), zid.marginalBenefitsDecreasing;
			stationSlotsOrderPairs = zid.stationSlotsOrderPairs)
	end

	ambStations = stations[movableAmbStations]

	return movableAmbs, ambStations
end

# solve as an integer program
function solveZhangIp(stationSlots::Vector{Int}, benefitSlots::Vector{Float}, ambToStationCosts::Array{Float,2},
	numMovableAmbs::Int, numAtHospitalAmbs::Int,
	marginalBenefitsDecreasing::Bool; stationSlotsOrderPairs::Array{Int,2} = Array{Int,2}(undef,0,0))::Tuple{Vector{Int}, Vector{Int}}
	# ambToStationCosts[:,j] gives costs to move different ambulances to station j
	# ambToStationCosts[1:numMovableAmbs,:] are for ambs that can be moved now
	# ambToStationCosts[(1:numAtHospitalAmbs).+numMovableAmbs,:] are for ambs that are at hospital
	# stationSlotsOrderPairs is not needed if marginalBenefitsDecreasing = true

	(numMoveUpAmbs, numStations) = size(ambToStationCosts)
	if numMoveUpAmbs == 0 return (Int[], Int[]) end
	numStationSlots = length(stationSlots)
	@assert(numMoveUpAmbs == numMovableAmbs + numAtHospitalAmbs)
	@assert(numStationSlots >= numMoveUpAmbs)
	@assert(numStationSlots == length(benefitSlots))
	@assert(all(in(1:numStations), stationSlots))

	if checkMode
		if marginalBenefitsDecreasing
			for j = 1:numStations
				# check that benefit of ambulance k at station j is > benefit of ambulance k+1, otherwise marginalBenefitsDecreasing should be false
				@assert(issorted(benefitSlots[findall(stationSlots .== j)], lt=<=, rev=true))
			end
		end
	end

	# shorthand:
	a = numMoveUpAmbs
	ai = 1:numMovableAmbs # row indices of ambToStationCosts that are for movable ambs
	aj = (1:numAtHospitalAmbs).+numMovableAmbs # row indices of ambToStationCosts that are for at-hospital ambs
	s = numStations
	n = length(stationSlots) # = length(benefitSlots)

	# using JuMP
	model = Model()
	set_optimizer(model, with_optimizer(GLPK.Optimizer)) # faster without presolve; have not compared with other solvers
	# set_optimizer(model, with_optimizer(Cbc.Optimizer, logLevel=0))

	@variable(model, x[i=1:a,j=1:s], Bin)
	marginalBenefitsDecreasing ? @variable(model, 0 <= y[k=1:n] <= 1) : @variable(model, y[k=1:n], Bin)
	# x[i,j] = 1 if ambulance moveUpAmbs[i] should be moved to station j
	# sum(y[stationSlots .== j]) = number of ambulances assigned to station j

	@constraints(model, begin
		(movableAmbAtOneLocation[i=ai], sum(x[i,:]) == 1) # each movable ambulance must be assigned to one station
		(atHospitalAmbAtMostOneLocation[i=aj], sum(x[i,:]) <= 1) # each at-hospital ambulance may be assigned to zero or one station
		(stationAmbCounts[j=1:s], sum(x[:,j]) == sum(y[k] for k=findall(stationSlots .== j)))
		# (stationAmbCounts[j=1:s], sum(x[:,j]) >= sum(y[k] for k=findall(stationSlots .== j))) # should have same effect as "==" constraint (instead of ">="), but may be faster?
	end)

	if !marginalBenefitsDecreasing
		# need to enforce station slot filling order
		p = stationSlotsOrderPairs # shorthand
		@constraint(model, stationSlotsFillingOrder[k=1:size(p,1)], y[p[k,1]] >= y[p[k,2]])
	end

	@expressions(model, begin
		totalBenefitAtStations, sum(y .* benefitSlots)
		totalAmbTravelCosts, sum(x .* ambToStationCosts)
	end)

	# solve
	@objective(model, Max, totalBenefitAtStations - totalAmbTravelCosts)
	@stdout_silent optimize!(model)
	@assert(termination_status(model) == MOI.OPTIMAL)

	# get solution
	vals = Dict()
	vals[:x] = JuMP.value.(x)
	vals[:y] = JuMP.value.(y)

	# solution
	sol = convert(Array{Bool,2}, round.(vals[:x]))
	movableAmbStations = [findfirst(sol[i,:]) for i = ai]
	atHospitalAmbStations = [findfirst(sol[i,:]) for i = aj]
	atHospitalAmbStations = convert(Vector{Int}, replace(atHospitalAmbStations, nothing => nullIndex))

	if checkMode
		@assert(all(sum(sol, dims=2) .<= 1)) # each ambulance can be used in move up at most once

		# check that y values are ordered correctly
		stationSlotsFilled = convert(Vector{Bool}, round.(vals[:y]))
		for j = 1:numStations
			@assert(issorted(stationSlotsFilled[stationSlots .== j], rev=true))
		end
	end

	return movableAmbStations, atHospitalAmbStations # atHospitalAmbStations[i] == nullIndex if amb i (of those at hospital) is not moved
end

# solve as an assignment problem
# this assumes that the benefit of adding amb k to a station is > benefit of adding amb k+1
function solveZhangIpAssignmentProblem(stationSlots::Vector{Int}, benefitSlots::Vector{Float}, ambToStationCosts::Array{Float,2},
	numMovableAmbs::Int, numAtHospitalAmbs::Int)::Tuple{Vector{Int}, Vector{Int}}
	# ambToStationCosts[:,j] gives costs to move different ambulances to station j
	# ambToStationCosts[1:numMovableAmbs,:] are for ambs that can be moved now
	# ambToStationCosts[(1:numAtHospitalAmbs).+numMovableAmbs,:] are for ambs that are at hospital

	(numMoveUpAmbs, numStations) = size(ambToStationCosts)
	if numMoveUpAmbs == 0 return (Int[], Int[]) end
	@assert(numMoveUpAmbs == numMovableAmbs + numAtHospitalAmbs)
	numStationSlots = length(stationSlots)
	@assert(numStationSlots >= numMoveUpAmbs)
	@assert(numStationSlots == length(benefitSlots))
	@assert(all(in(1:numStations), stationSlots))

	if checkMode && false # skip this check, it is slow
		for j = 1:numStations
			# check that benefit of ambulance k at station j is > benefit of ambulance k+1, otherwise cannot solve as assignment problem
			@assert(issorted(benefitSlots[findall(stationSlots .== j)], lt=<=, rev=true))
		end
	end

	# shorthand
	a = numMoveUpAmbs
	ai = 1:numMovableAmbs # row indices of ambToStationCosts that are for movable ambs
	aj = (1:numAtHospitalAmbs).+numMovableAmbs # row indices of ambToStationCosts that are for at-hospital ambs
	n = numStationSlots + numAtHospitalAmbs # total number of slots, need one dummy per each amb at hospital
	si = 1:numStationSlots # column indices of ambToStationCosts that are for real stations
	sj = numStationSlots+1:n # column indices of ambToStationCosts that are for dummy stations

	# formulate assignment problem and solve
	@assert(a <= n) # needed for Hungarian algorithm
	weights = zeros(Float, a, n)
	for i = 1:a, j = si
		weights[i,j] = -(benefitSlots[j] - ambToStationCosts[i,stationSlots[j]]) # note that Hungarian algorithm looks for min cost assignment
	end
	weights[ai, sj] .= Inf # movable ambs cannot be assigned to dummy stations
	matching = Hungarian.munkres(weights) # returns a sparse matrix

	# extract solution
	I = findall(matching .== Hungarian.STAR)
	(ambIndices, stationSlotIndices) = (getindex.(I,1), getindex.(I,2))
	ambStations = zeros(Int,a) # will be the station (real or dummy) of each ambulance
	stationSlotsExtended = vcat(stationSlots, fill(nullIndex, numAtHospitalAmbs)) # add dummy station indices
	ambStations[ambIndices] = stationSlotsExtended[stationSlotIndices]
	movableAmbStations = ambStations[ai]
	atHospitalAmbStations = ambStations[aj]

	if checkMode
		@assert(!any(isequal(nullIndex), movableAmbStations)) # make sure that no movable ambs are assigned to dummy stations

		if false # skip this check, it is slow
			# check that station slots are filled in correct order
			stationSlotsFilled = fill(false, numStationSlots)
			for i in stationSlotIndices
				if i <= numStationSlots stationSlotsFilled[i] = true end
			end
			for j = 1:numStations
				@assert(issorted(stationSlotsFilled[findall(stationSlots .== j)], rev=true))
			end
		end
	end

	return movableAmbStations, atHospitalAmbStations # atHospitalAmbStations[i] == nullIndex if amb i (of those at hospital) is not moved
end


text:
#Function to minimize Bregman divergences h() with respect to a given
#vector y over the submodular base polytope of a cardinality-based submodular function

# using Distributions

export card_inc_fix

function card_inc_fix(F::CardBasedAtom,
	                    w₀::AbstractArray,
										  divergence::String = "euclidean")
  # here F(S) = g(|S|) for some non-decreasing concave function g, g(0)=0.
  # currently implemented for divergences: Euclidean
  # check if g is concave, non-decreasing, g is an Array for now

  n = length(w₀)
  g = zeros(n)
  for i = 1:n
    g[i] = evaluate(F, collect(1:i))[1]
  end

	sorted_indices = sortperm(w₀, rev=true)
	w = w₀[sorted_indices]
	if(w[1] <= 0 && divergence != "euclidean")
   	error("The input is out of range for itakura-saito, logistic, entropy")
	elseif(w[1] >= 1&&divergence == "logistic")
  	error("The input is out of range for logistic")
	end

	counter = 0                             # maintains the tight set
	epsilon = zeros(n,1)

	if(divergence == "euclidean")
    x_opt = zeros(n,1)
  	# x_opt = round.(x_opt, 4)
  	feas_x_opt = (g[n]/n) * ones(n,1)
  	while(counter < n)
      for j = (counter + 1) : n
        if(counter == 0)
          epsilon[j] = (g[j] - sum(w[1:j]))/j
        else
          epsilon[j] = (g[j] - g[counter] - sum(w[counter+1:j]))/(j-counter)
        end
      end
      (e_k, I) = findmin(epsilon[counter+1:n])
      k = I + counter
    	x_opt[counter+1:k] = w[counter+1: k] + e_k
      counter = k
    	temp = counter+1
      while(temp<=n && x_opt[temp] - w[temp] <= e_k)
        x_opt[temp] = e_k+ w[temp]
        temp = temp+1
      end
    end
	end
	temp = zeros(n,1)
	for i =1:n
		temp[sorted_indices[i]] = x_opt[i]
	end
	return temp
end

card_inc_fix(p::AssocPoly, w::AbstractArray, divergence::String = "euclidean") = card_inc_fix(p.F, w, divergence)


text:
using SCIP
using MathProgBase
using FactCheck

include("more_tests.jl")
include("mixintprog.jl")

include(joinpath(Pkg.dir("JuMP"),"test","solvers.jl"))

using JuMP
solver = SCIPSolver("display/verblevel", 0)
heursolver = SCIPSolver("display/verblevel", 0,
                        "presolving/maxrounds", 0,
                        "separating/maxrounds", 0,
                        "separating/maxroundsroot", 0)
nlsolver = SCIPSolver("display/verblevel", 0)
# "numerics/feastol", 1e-7) # doesn't work, no convert from Ptr{Void} to Float64
if length(lp_solvers) < 1
    error("Need at least one LP solver, eg, Clp")
end
lp_solvers = lp_solvers[1:1] # just need one
ip_solvers = [solver]
quad_solvers = []
quad_mip_solvers = []
quad_mip_solvers = [solver]
lazy_solvers = [solver]
lazylocal_solvers = []
cut_solvers = []
cutlocal_solvers = []
heur_solvers = [heursolver]
info_solvers = []
minlp_solvers = [nlsolver]

# nlp_solvers fails becuase:
# - some dual information: fixed by deleting getconstrduals
# - some strange exprresions: can be fixed by adding SCIPSolver to the excluded
#   list in "Test ifelse". This is on JuMP's side
# - numerical tolerances: I tried setting SCIP feastol to 1e-9 and still got
#   numerical issues. However, the solution is feasible for SCIP.
#nlp_solvers = [nlsolver]

include(joinpath(Pkg.dir("JuMP"),"test","model.jl"))
include(joinpath(Pkg.dir("JuMP"),"test","probmod.jl"))
include(joinpath(Pkg.dir("JuMP"),"test","qcqpmodel.jl"))
include(joinpath(Pkg.dir("JuMP"),"test","callback.jl"))
include(joinpath(Pkg.dir("JuMP"),"test","nonlinear.jl"))
#include(joinpath(Pkg.dir("JuMP"),"test","sdp.jl"))

# return proper code (e.g. for TravisCI)
FactCheck.exitstatus()


text:
#
# Correctness Tests
#

using Base.Test

my_tests = ["naive.jl", "kdtree.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
@testset "Test composition order" begin
  Qx, x = FlintQQ["x"]
  K, a = NumberField(x^6 + 108, "a")
  A = automorphisms(K)
  for f in A
    for g in A
      @test g(f(a)) == compose(f, g)(a)
      @test g(f(a)) == (f * g)(a)
    end
  end
end

@testset "Induce image" begin
  Qx, x = FlintQQ["x"]
  K, a = NumberField(x^6 + 108, "a")
  A = automorphisms(K)
  OK = maximal_order(K)
  lP = prime_decomposition(OK, 7)
  lP1 = [x[1] for x in lP]
  for i = 1:length(lP)
    I_new = A[2](lP1[1])
    id = findfirst(isequal(I_new), lP1)
    @test id != nothing
  end
  f = hom(K, K, a^4//12+a//2)
  E = EquationOrder(K)
  I = ideal(E, E(a))
  @test_throws ErrorException Hecke.induce_image(f, I)

end


text:
## test the main entry point in Celeste: the `infer` function
import JLD

import Celeste: Config, BoundingBox
using Celeste.SDSSIO
using Celeste.ParallelRun
import Celeste.ParallelRun: OptimizedSource
using Celeste.DeterministicVI

function compute_obj_value(images::Vector{<:Image},
                           catalog::Vector{CatalogEntry},
                           box::BoundingBox,
                           results::Vector{OptimizedSource})

    # TODO: This stuff is duplicated from ParallelRun.infer_box.
    # We should refactor infer_box to return the objective value in some way!
    patches = Model.get_sky_patches(images, catalog)
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_ids = find(entry_in_range, catalog)

    # There must be a vp for every patch in the call to elbo().
    # So, here we must limit patches to just the targets we optimized
    # and pass [1, 2, 3, ...] as the target indexes.
    ea = ElboArgs(images, patches[target_ids, :],
                  collect(1:length(target_ids)); include_kl=false)
    vp = [r.vs for r in results]
    DeterministicVI.elbo(ea, vp).v[]
end

@testset "infer" begin
    @testset "infer_box(...; method = :single_vi) runs" begin
        # very small patch of sky that turns out to have 4 sources.
        # We checked that this patch is in the given field.
        box = BoundingBox(164.39, 164.41, 39.11, 39.13)
        images = SampleData.get_sdss_images(3900, 6, 269)
        result = ParallelRun.infer_box(images, box; method=:single_vi)
    end

    @testset "infer_box(...; method = :mcmc) runs" begin
        # very small patch of sky that turns out to have 4 sources.
        # We checked that this patch is in the given field.
        box = BoundingBox(164.39, 164.41, 39.11, 39.13)
        images = SampleData.get_sdss_images(3900, 6, 269)
        result = ParallelRun.infer_box(images, box; method=:mcmc,
                                       config=Config(2.0, 3, 2, 3))
    end

    @testset "joint vs single on overlapping sources" begin
        images = SampleData.get_sdss_images(4263, 5, 119)
        catalog = SampleData.get_sdss_catalog(4263, 5, 119)

        # This box has 3 overlapping objects in it.
        box = BoundingBox(0.467582, 0.473275, 0.588383, 0.595095)

        results_single = ParallelRun.infer_box(images, box;
                                               catalog = catalog,
                                               method=:single_vi)
        results_joint = ParallelRun.infer_box(images, box;
                                              catalog = catalog,
                                              method=:joint_vi)

        @test length(results_single) == 3
        @test length(results_joint) == 3

        score_single = compute_obj_value(images, catalog, box, results_single)
        score_joint = compute_obj_value(images, catalog, box, results_joint)

        @test score_joint > score_single
    end
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames: DataFrame

include("../helper/gpu.jl")

@testset "Likelihood" begin
    # Uses a PC with 4 variables, and tests 3 of the configurations to
    # match with python. Also tests all probabilities sum up to 1.

    prob_circuit = zoo_psdd("little_4var.psdd");
    @test prob_circuit isa ProbCircuit;

    # Step 1. Check Probabilities for 3 samples
    data = DataFrame(BitArray([0 0 0 0; 0 1 1 0; 0 0 1 1]));
    true_prob = [0.07; 0.03; 0.13999999999999999]

    calc_prob = EVI(prob_circuit, data)
    calc_prob = exp.(calc_prob)

    @test true_prob ≈ calc_prob atol=1e-7;

    # Step 2. Add up all probabilities and see if they add up to one
    N = 4;
    data_all = generate_data_all(N)

    calc_prob_all = EVI(prob_circuit, data_all)
    calc_prob_all = exp.(calc_prob_all)
    sum_prob_all = sum(calc_prob_all)

    @test 1 ≈ sum_prob_all atol = 1e-7;

    cpu_gpu_agree_approx(data_all) do d
        EVI(prob_circuit, d)
    end

    # make sure log-likelihoods are -Inf when the input is not satisfied
    alltrue = multiply(pos_literals(ProbCircuit,4))
    @test all(EVI(alltrue, data) .== -Inf)

    cpu_gpu_agree(data) do d
        EVI(alltrue, d)
    end

end

@testset "Bagging models' likelihood" begin
    dfb = DataFrame(BitMatrix([true true; true true; true true; true true]))
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))
    # bag_dfb = bagging_dataset(dfb; num_bags = 2, frac_examples = 1.0)
    bag_dfb = Array{DataFrame}(undef, 2)
    bag_dfb[1] = dfb[[2, 1, 3, 4], :]
    bag_dfb[2] = dfb[[4, 3, 2, 1], :]

    r = compile(SharedProbCircuit, r, 2)

    estimate_parameters(r, bag_dfb; pseudocount = 1.0)

    ll = log_likelihood_per_instance(r, dfb)
    @test ll[1] ≈ -0.2107210 atol = 1e-6
    @test ll[2] ≈ -0.2107210 atol = 1e-6
    @test ll[3] ≈ -0.2107210 atol = 1e-6
    @test ll[4] ≈ -0.2107210 atol = 1e-6

    @test log_likelihood(r, dfb) ≈ -0.8428841 atol = 1e-6

    @test log_likelihood_avg(r, dfb) ≈ -0.2107210 atol = 1e-6
end


text:
using Documenter, AbstractPlotting
using Markdown, Pkg, Random, FileIO
using MakieGallery, Makie
import AbstractPlotting: _help, to_string, to_func, to_type
using MakieGallery: eval_examples, generate_thumbnail, master_url, print_table
using MakieGallery: @cell, @block, @substep

cd(@__DIR__)
database = MakieGallery.load_database()

pathroot = normpath(@__DIR__, "..")
docspath = joinpath(pathroot, "docs")
srcpath = joinpath(docspath, "src")
buildpath = joinpath(docspath, "build")
mediapath = joinpath(homedir(), "ReferenceImages", "gallery")

# =============================================
# automatically generate an overview of the atomic functions, using a source md file
@info("Generating functions overview")
path = joinpath(srcpath, "functions-overview.md")
srcdocpath = joinpath(srcpath, "src-functions.md")

plotting_functions = (
    AbstractPlotting.atomic_functions..., contour, arrows,
    barplot, poly, band, slider, vbox
)

open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plotting functions overview")
    src = read(srcdocpath, String)
    println(io, src, "\n")
    for func in plotting_functions
        fname = to_string(func)
        println(io, "## `$fname`\n")
        println(io, "```@docs")
        println(io, "$fname")
        println(io, "```\n")
        help_attributes(io, func; extended = true)
        # add previews of all tags related to function
        for example in database
            fname in example.tags || continue
            base_path = joinpath(mediapath, string(example.unique_name))
            thumb = master_url(mediapath, joinpath(base_path, "media", "thumb.jpg"))
            code = master_url(mediapath, joinpath(base_path, "index.html"))
            src_lines = example.file_range
            println(io, """[![library lines $src_lines]($thumb)]($code)""")
        end
        println(io, "\n")
    end
end

cd(docspath)

# =============================================
# automatically generate an overview of the plot attributes (keyword arguments), using a source md file
@info("Generating attributes page")
include("../src/plot_attr_desc.jl")
path = joinpath(srcpath, "plot-attributes.md")
srcdocpath = joinpath(srcpath, "src-plot-attributes.md")
open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plot attributes")
    src = read(srcdocpath, String)
    println(io, src)
    print(io, "\n")
    println(io, "## List of attributes")
    print_table(io, plot_attr_desc)
end

# =============================================
# automatically generate an overview of the axis attributes, using a source md file
@info("Generating axis page")
path = joinpath(srcpath, "axis.md")
srcdocpath = joinpath(srcpath, "src-axis.md")
include("../src/Axis2D_attr_desc.jl")
include("../src/Axis3D_attr_desc.jl")

open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    src = read(srcdocpath, String)
    println(io, src)
    print(io)
    # Axis2D section
    println(io, "## `Axis2D`")
    println(io, "### `Axis2D` attributes groups")
    print_table(io, Axis2D_attr_desc)
    print(io)
    for (k, v) in Axis2D_attr_groups
        println(io, "#### `:$k`\n")
        print_table(io, v)
        println(io)
    end
    # Axis3D section
    println(io, "## `Axis3D`")
    println(io, "### `Axis3D` attributes groups")
    print_table(io, Axis3D_attr_desc)
    print(io)
    for (k, v) in Axis3D_attr_groups
        println(io, "#### `:$k`\n")
        print_table(io, v)
        println(io)
    end
end

# automatically generate an overview of the function signatures, using a source md file
@info("Generating signatures page")
path = joinpath(srcpath, "signatures.md")
srcdocpath = joinpath(srcpath, "src-signatures.md")
open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plot function signatures")
    src = read(srcdocpath, String)
    println(io, src)
    print(io, "\n")
    println(io, "```@docs")
    println(io, "convert_arguments")
    println(io, "```\n")

    println(io, "See [Plot attributes](@ref) for the available plot attributes.")
end

# build docs with Documenter
@info("Running `makedocs` with Documenter.")

using Documenter, MakieGallery, Makie, AbstractPlotting

makedocs(
    modules = [Makie, AbstractPlotting],
    doctest = false, clean = true,
    format = Documenter.HTML(prettyurls = false),
    sitename = "Makie.jl",
    pages = Any[
        "Home" => "index.md",
        "Basics" => [
            "basic-tutorials.md",
            "statsmakie.md",
            "help_functions.md",
            "functions-overview.md",
            "signatures.md",
            "plot-attributes.md",
            "colors.md",
            "theming.md",
            # "extending.md",
            "axis.md",
            "interaction.md",
            "recipes.md",
            "output.md",
            # "layout.md",
            "troubleshooting.md"
        ],
        "Developer Documentation" => [
            "why-makie.md",
        #     "devdocs.md",
            "abstractplotting_api.md"
        ],
    ]
)

using Conda, Documenter
using Base64
# deploy
ENV["DOCUMENTER_DEBUG"] = "true"
if !haskey(ENV, "DOCUMENTER_KEY")
    # Workaround for when deploying locally and silly Windows truncating the env variable
    # on the CI these should be set!
    ENV["TRAVIS_BRANCH"] = "latest"
    ENV["TRAVIS_PULL_REQUEST"] = "false"
    ENV["TRAVIS_REPO_SLUG"] = "github.com/JuliaPlots/MakieGallery.jl.git"
    ENV["TRAVIS_TAG"] = "v7"
    ENV["TRAVIS_OS_NAME"] = ""
    ENV["TRAVIS_JULIA_VERSION"] = ""
    ENV["PATH"] = string(ENV["PATH"], Sys.iswindows() ? ";" : ":", Conda.SCRIPTDIR)
    ENV["DOCUMENTER_KEY"] = readchomp(joinpath(homedir(), "documenter.key"))
end

run(`pip install --upgrade pip`)
cd(@__DIR__)
deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-cinder"),
    repo = "github.com/JuliaPlots/MakieGallery.jl.git",
)


text:
module AuraLighting

using ZMQ

export AuraMbControl, AuraControlClient, getcolor, setcolor, setmode
export rbgtoi, itorbg, startserver, iscorrectcontroller, sendexit
export AuraGPUControl, AuraKeyboardControl, AuraMouseControl, EnumerateControls

const Handle = Ptr{Cvoid}
const Hptr = Ptr{Ptr{Cvoid}}
const Bptr = Ptr{UInt8}

"""
NOTE: THIS DLL FILE MUST BE EITHER BE LOCATED IN THE DIRECTORY WHERE PROGRAM IS RUN,
OR ELSE YOU SHOULD PLACE A COPY IN A DIRECTORY IN YOUR PATH.
"""
const DLLNAME = "AURA_SDK.dll"

""" rgb integer to UInt8 (r, g, b) """
itorbg(i) = [UInt8((i >> 16) & 0xff), UInt8((i >> 8) & 0xff), UInt8(i & 0xff)]

""" UInt8 (r, g, b) to rgb integer """
rbgtoi(r, g, b) = ((UInt32(r) << 16) | (UInt32(g) << 8) | UInt32(b)) & 0xffffff

abstract type AuraControl end

struct AuraMbControl <: AuraControl
    controllernumber::Int
    LEDcount::Int
    handle::Handle
    colorbuf::Vector{UInt8}
    buflen::Int
    port::Int
    AuraMbControl(c, n, h, p) = new(c, n, h, zeros(UInt, n * 3), n * 3, p)
end

"""
    function AuraMbControl(cont=1; asservice=false, port=5555)

Constructor for an AuraMBControl.
cont: controller number, defaults to 1 (first or only controller found)
port: port number of ZMQ service, defaults to 5555
"""
function AuraMbControl(cont=1, port=5555)
    GC.enable(false)
    hcount = ccall((:EnumerateMbController, DLLNAME), Cint, (Hptr, Cint), C_NULL, 0)
    handles = [C_NULL for _ in 1:hcount]
    GC.@preserve handles begin
        ccall((:EnumerateMbController, DLLNAME), Cint, (Hptr, Cint), handles, hcount)
        handle = handles[hcount]
        LEDcount = ccall((:GetMbLedCount, DLLNAME), Cint, (Handle,), handle)
        return AuraMbControl(cont, LEDcount, handle, port)
    end
end

struct AuraGPUControl <: AuraControl
    controllernumber::Int
    LEDcount::Int
    handle::Handle
    colorbuf::Vector{UInt8}
    buflen::Int
    port::Int
    AuraGPUControl(c, n, h, p) = new(c, n, h, zeros(UInt, n * 3), n * 3, p)
end

"""
    function AuraGPUControl(cont=1, port=5556)

Constructor for an AuraGPUControl.
cont: controller number, defaults to 1 (first or only controller found)
port: port number of ZMQ service, defaults to 5556
"""
function AuraGPUControl(cont=1, port=5556)
    GC.enable(false)
    hcount = ccall((:EnumerateGPU, DLLNAME), Cint, (Hptr, Cint), C_NULL, 0)
    handles = [C_NULL for _ in 1:hcount]
    GC.@preserve handles begin
        ccall((:EnumerateGPU, DLLNAME), Cint, (Hptr, Cint), handles, hcount)
        handle = handles[max(min(hcount, cont), 1)]
        LEDcount = ccall((:GetGPULedCount, DLLNAME), Cint, (Handle,), handle)
        return AuraGPUControl(cont, LEDcount, handle, port)
    end
end

struct AuraKeyboardControl <: AuraControl
    LEDcount::Int
    handle::Handle
    colorbuf::Vector{UInt8}
    buflen::Int
    port::Int
    AuraKeyboardControl(n, h, p) = new(n, h, zeros(UInt, n * 3), n * 3, p)
end

"""
    function AuraKeyboardControl(cont=1; asservice=false, port=5557)
Constructor for an AuraKeyboardControl.
port: port number of ZMQ service, defaults to 5557
"""
function AuraKeyboardControl(port=5557)
    GC.enable(false)
    handles = [C_NULL]
    GC.@preserve handles begin
        ccall((:CreateClaymoreKeyboard, DLLNAME), Cint, (Hptr,), handles)
        handle = handles[1]
        LEDcount = ccall((:GetClaymoreKeyboardLedCount, DLLNAME), Cint, (Handle,), handle)
        return AuraKeyboardControl(LEDcount, handle, port)
    end
end

struct AuraMouseControl <: AuraControl
    LEDcount::Int
    handle::Handle
    colorbuf::Vector{UInt8}
    buflen::Int
    port::Int
    AuraMouseControl(n, h, p) = new(n, h, zeros(UInt, n * 3), n * 3, p)
end

"""
    function AuraMouseControl(cont=1; asservice=false, port=5558)
Constructor for an AuraMouseControl.
port: port number of ZMQ service, defaults to 5558
"""
function AuraMouseControl(port=5558)
    GC.enable(false)
    handles = [C_NULL]
    GC.@preserve handles begin
        ccall((:CreateRogMouse, DLLNAME), Cint, (Hptr,), handles)
        handle = handles[1]
        LEDcount = ccall((:GetRogLedCount, DLLNAME), Cint, (Handle,), handle)
        return AuraMouseControl(LEDcount, handle, port)
    end
end

function startserver(au::AuraControl)
    @async begin ZMQservice(au) end
end

""" only motherboard and GPU have controller number, all others are just 1 """
controllernumber(au::AuraMbControl) = au.controllernumber
controllernumber(au::AuraGPUControl) = au.controllernumber
controllernumber(au::AuraControl) = 1

"""
   function setmode(au::AuraMbControl, setting::Integer=0)

Set mode of motherboard Aura controller from software control to an auto mode

A setting of 0 will change to the auto mode. A setting of 1 is software control.
"""
function setmode(au::AuraMbControl, setting::Integer)
    0 <= setting <= 1 || return false
    h = au.handle
    return ccall((:SetMbMode, DLLNAME), Cint, (Handle, Cint), h, setting)
end

function setmode(au::AuraGPUControl, setting::Integer)
    0 <= setting <= 1 || return false
    h = au.handle
    return ccall((:SetGPUMode, DLLNAME), Cint, (Handle, Cint), h, setting)
end

function setmode(au::AuraKeyboardControl, setting::Integer)
    0 <= setting <= 1 || return false
    h = au.handle
    return ccall((:SetClaymoreKeyboardMode, DLLNAME), Cint, (Handle, Cint), h, setting)
end

function setmode(au::AuraMouseControl, setting::Integer)
    0 <= setting <= 1 || return false
    h = au.handle
    return ccall((:SetRogMouseMode, DLLNAME), Cint, (Handle, Cint), h, setting)
end

"""
    function getcolor(auramb::AuraMbControl)

Get RGB color as a tuple of red, green, and blue values (0 to 255 each).
"""
function getcolor(au::AuraMbControl)
    for i in 1:au.buflen
        au.colorbuf[i] = 0x0
    end
    buf = au.colorbuf
    GC.@preserve buf begin
        ccall((:GetMbColor, DLLNAME), Cint, (Handle, Bptr, Cint),
            au.handle, buf, au.buflen)
        return Int(buf[1]), Int(buf[2]), Int(buf[3])
    end
end

""" only motherboard currently can get color, though all types can set color """
function getcolor(au::AuraControl)
    return (-1, -1, -1)
end

"""
    function setcolor(au::AuraMbControl, red, green, blue)

Set RGB color via setting color with separate red, green, and blue values
"""
function setcolor(au::AuraMbControl, red, green, blue)
    for i in 1:3:au.buflen-1
        au.colorbuf[i], au.colorbuf[i+1], au.colorbuf[i+2] = red, green, blue
    end
    buf = au.colorbuf
    GC.@preserve buf begin
        success = ccall((:SetMbColor, DLLNAME), Cint, (Handle, Bptr, Cint),
            au.handle, buf, au.buflen)
        return success == 1
    end
end

"""
    function setcolor(au::AuraGPUControl, red, green, blue)

Set RGB color via setting color with separate red, green, and blue values
"""
function setcolor(au::AuraGPUControl, red, green, blue)
    for i in 1:3:au.buflen-1
        au.colorbuf[i], au.colorbuf[i+1], au.colorbuf[i+2] = red, green, blue
    end
    buf = au.colorbuf
    GC.@preserve buf begin
        success = ccall((:SetGPUColor, DLLNAME), Cint, (Handle, Bptr, Cint),
            au.handle, buf, au.buflen)
        return success == 1
    end
end

"""
    function setcolor(au::AuraKeyboardControl, red, green, blue)

Set RGB color via setting color with separate red, green, and blue values
"""
function setcolor(au::AuraKeyboardControl, red, green, blue)
    for i in 1:3:au.buflen-1
        au.colorbuf[i], au.colorbuf[i+1], au.colorbuf[i+2] = red, green, blue
    end
    buf = au.colorbuf
    GC.@preserve buf begin
        success = ccall((:SetClaymoreKeyboardColor, DLLNAME), Cint, (Handle, Bptr, Cint),
            au.handle, buf, au.buflen)
        return success == 1
    end
end

"""
    function setcolor(au::AuraControl, red, green, blue)

Set RGB color via setting color with separate red, green, and blue values
"""
function setcolor(au::AuraMouseControl, red, green, blue)
    for i in 1:3:au.buflen-1
        au.colorbuf[i], au.colorbuf[i+1], au.colorbuf[i+2] = red, green, blue
    end
    buf = au.colorbuf
    GC.@preserve buf begin
        success = ccall((:SetRogMouseColor, DLLNAME), Cint, (Handle, Bptr, Cint),
            au.handle, buf, au.buflen)
        return success == 1
    end
end

"""
    function setcolor(au::AuraControl, rgb::Integer)

Set color of the LED light as an integer in 24-bit RGB format.
The color is of form hex 0xRRGGBB, where RR is the red component, GG green,
    and BB the blue values of an 24-bit RGB coded color.
Black is 0, white is 0x00ffffff, red 0xff0000, green 0x00ff00, blue 0x0000ff
"""
function setcolor(au::AuraControl, rgb)
    r, g, b = itorbg(rgb)
    setcolor(au, r, g, b)
end

"""
    function ZMQservice(au::AuraControl)

Serve requests via ZMQ to control the Aura lighting controller on the motherboard.
This must be run in 32-bit Windows mode with admin privileges.
"""
function ZMQservice(au::AuraControl)
    sock = Socket(REP)
    bind(sock, "tcp://*:$(au.port)")
    try
        while true
            message = recv(sock, String)
            cmd = message * " 0 0 "
            words = split(cmd, r"\s+")
            if words[1] == "getcolor"
                r, g, b = getcolor(au)
                if r == g == b == -1
                    send("ERROR not supported")
                else
                    send(sock, "OK $(rbgtoi(r, g, b))")
                end
            elseif words[1] == "setcolor" && (c = tryparse(Int, words[2])) != nothing
                setcolor(au, c)
                send(sock, "OK")
            elseif words[1] == "setmode" && (n = tryparse(Int, words[2])) != nothing
                setmode(au, n)
                send(sock, "OK")
            elseif words[1] == "getcontroller"
                send(sock, "OK $(controllernumber(au))")
            elseif words[1] == "exit"
                send(sock, "OK")
                break
            else
                @warn("Unknown command received: $message")
                send(sock, "ERROR in message received: $message")
            end
        end
    catch y
        @warn("ZMQ server fatal error: $y with message $message")
    finally
        close(sock)
    end
end

"""
    EnumerateControls()
Find any accessible Aura lighting AuraControl objects and return a set of such,
or an empty vector if none are found. Note will also return empty if the
Windows 32-bit DLL AURA_SDK.dll is not found (in current directory or in a directory
previously set as a valid DLL directory) or if it is not loadable (eg, not 32-bit Windows Julia).
"""
function EnumerateControls(verbose=true)
    found = AuraLighting.AuraControl[]
    begin
        try
            auramb = AuraMbControl()
            push!(found, auramb)
            verbose && @warn("Found a motherboard Aura controller")
        catch y
            verbose && @warn("No motherboard controller found: $y")
        end
    end
    begin
        try
            auragpu = AuraGPUControl()
            push!(found, auragpu)
            verbose && @warn("Found a GPU Aura controller")
        catch y
            verbose && @warn("No GPU controller found: $y")
        end
    end
    begin
        try
            aurakb = AuraKeyboardControl()
            push!(found, aurakb)
            verbose && @warn("Found a keyboard Aura controller")
        catch y
            verbose && @warn("No keyboard Aura controller found: $y")
        end
    end
    begin
        try
            auramouse = AuraMouseControl()
            push!(found, auramouse)
            verbose && @warn("Found a mouse Aura controller")
        catch y
            verbose && @warn("No mouse Aura controller found: $y")
        end
    end
    return found
end

"""
    struct AuraControlClient

Client for a service to change Aura lighting via the 32-bit Win32 Aura lighting SDK.

This client can be any application and can be 64-bit even though the server must be
32-bit because ASUStek only provided a 32-bit DLL for Windows in its AuraSDK library.
"""
struct AuraControlClient
    sock::Socket
    controllernumber::Int
    function AuraControlClient(port=5555, server="localhost", cont=1)
        sock = Socket(REQ)
        connect(sock, "tcp://$server:$port")
        new(sock, cont)
    end
end

"""
    function iscorrectcontroller(client::AuraControlClient)

Check if the client's controller number matches the server's controller number.
It is not actually necessary these two match, but checking this may help avoid
sending commands to the wrong controller.
"""
function iscorrectcontroller(client::AuraControlClient)
    send(client.sock, "getcontroller")
    message = ZMQ.recv(client.sock, String)
    try
        s = split(message, r"\s+")
        n = tryparse(Int, s[2])
        return n == client.controllernumber
    catch
        return false
    end
end

"""
    function getcolor(client::AuraMBControlClient)

Get Aura lighting color as a tuple of red, green, and blue.
"""
function getcolor(client::AuraControlClient)
    try
        send(client.sock, "getcolor")
        message = ZMQ.recv(client.sock, String)
        s = split(message, r"\s+")
        if s[1] != "OK"
            @warn("Error, getcolor may not be supported")
        return (0, 0, 0)
        else
            c = parse(Int, s[2])
            return Tuple(itorbg(c))
        end
    catch y
        @warn("Error getting color: $y")
        return (0, 0, 0)
    end
end

"""
    function setcolor(client::AuraControlClient, color::Integer)

Set Aura lighting color to an RGB integer of form 0xrrggbb.
Return: true on success, false on failure
"""
function setcolor(client::AuraControlClient, color)
    try
        send(client.sock, "setcolor $(color & 0xffffff)")
        message = ZMQ.recv(client.sock, String)
        return message[1:2] == "OK"
    catch y
        @warn("Error setting color to $color: $y")
        return false
    end
end

"""
    function setcolor(client::AuraControlClient, r, g, b)

Set Aura lighting color to RGB color with components r red, g green, b blue.
Return: true on success, false on failure
"""
setcolor(client::AuraControlClient, r, g, b) = setcolor(client, rbgtoi(r, g, b))

"""
   function setmode(client::AuraControlClient, mode::Integer)

Set the mode of the controller to 0 for auto mode, 1 for software controlled.
Return true on success
"""
function setmode(client::AuraControlClient, mode)
    0 <= mode <= 1 || return false
    try
        send(client.sock, "setmode $mode")
        message = recv(client.sock, String)
        return message[1:2] == "OK"
    catch y
        @warn("Error setting mode")
        return false
    end
end

function sendexit(client::AuraControlClient)
    try
        send(client.sock, "exit")
        message = recv(client.sock, String)
        if message[1:2] == "OK"
            @info("Sent exit command to server, terminating socket")
            close(client.sock)
            return true
        end
    catch y
        @warn("Error sending exit command to server")
        return false
    end
end

end # module


text:
# using Pkg
# Pkg.add("DataFramesMeta")

using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using CSV
using DataFramesMeta

function q3g()
# Question 1
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdataBeta0.csv"
df = CSV.read(HTTP.get(url).body)

df = @transform(df, bus_id = 1:size(df,1))

dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

# Question 2
sigma_hat = glm(@formula(Y ~ Odometer + Branded ), df_long, Binomial(), LogitLink())
println(sigma_hat)

# Question 3a
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body)

Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
X1 = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
X2 = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])

# Question 3b
function create_grids()

    function xgrid(theta,xval)
        N      = length(xval)
        xub    = vcat(xval[2:N],Inf)
        xtran1  = zeros(N,N)
        xtran1c = zeros(N,N)
        lcdf   = zeros(N)
        for i=1:length(xval)
            xtran1[:,i]   = (xub[i] .>= xval).*(1 .- exp.(-theta*(xub[i] .- xval)) .- lcdf)
            lcdf        .+= xtran1[:,i]
            xtran1c[:,i] .+= lcdf
        end
        return xtran1, xtran1c
    end

    zval   = collect([.25:.01:1.25]...)
    zbin   = length(zval)
    xval   = collect([0:.125:25]...)
    xbin   = length(xval)
    tbin   = xbin*zbin
    xtran  = zeros(tbin,xbin)
    xtranc = zeros(xbin,xbin,xbin)
    for z=1:zbin
        xtran[(z-1)*xbin+1:z*xbin,:],xtranc[:,:,z] = xgrid(zval[z],xval)
    end

    return zval,zbin,xval,xbin,xtran
end
zval,zbin,xval,xbin,xtran = create_grids()

@views @inbounds function q3e()
# Question 3c
T=20
future_value_array=zeros(size(xtran,1),2,T+1)
β=0.9

for t=T+1:1
  for b=1:size(df,1)
   for z=1:zbin
    for x=1:xbin
       row=xval[x] + (zval[z]-1)*xbin
       v1[t]= 1.92596-0.148154*xval[x]+1.05919*Branded[b]+xtran[row,:]'*FV[(zval[z]-1)*xbin+1:zval[z]*xbin,b+1,t+1]
       v0[t]= xtran[1+(zval[z]-1)*xbin,:]'*FV[(zval[z]-1)*xbin+1:zval[z]*xbin,b+1,t+1]
       futurevaluearray[:,:,t]=β*log(exp(v0[t])+exp(v1[t]))
end
end
end
end

# Q3d
log_likelihood_value=zeros(T+1)
for t=T+1:1
row0=1+(:Zst-1)*xbin
row1=:Xst+(:Zst-1)*xbin
v1t-v0t
(xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,B[i]+1,t+1]

P_i1t =exp(v1t −v0t)/(1+exp(v1t −v0t))
P_i0t = 1−P_i1t

log_likelihood_value[t]=d_ijt*log(P_ijt)
loglikehood=sum(log_likelihood_value)
end

return -loglikehood
end

end

# Question 3e
println(q3e())

# Question 3g
println(q3g())

max_stars_repo_path: ProblemSets/PS5-ddc/PS5_Dayuan Xie.jl

max_issues_repo_path: ProblemSets/PS5-ddc/PS5_Dayuan Xie.jl

max_forks_repo_path: ProblemSets/PS5-ddc/PS5_Dayuan Xie.jl


text:
module Grids

export Center, Face
export AbstractTopology, Periodic, Bounded, Flat, Connected, topology

export AbstractGrid, AbstractUnderlyingGrid, halo_size, total_size
export AbstractRectilinearGrid, RectilinearGrid
export XRegRectilinearGrid, YRegRectilinearGrid, ZRegRectilinearGrid, HRegRectilinearGrid, RegRectilinearGrid
export AbstractCurvilinearGrid, AbstractHorizontallyCurvilinearGrid
export LatitudeLongitudeGrid, XRegLatLonGrid, YRegLatLonGrid, ZRegLatLonGrid
export ConformalCubedSphereFaceGrid, ConformalCubedSphereGrid
export node, xnode, ynode, znode, xnodes, ynodes, znodes, nodes
export offset_data, new_data
export on_architecture

using CUDA
using Adapt
using OffsetArrays

using Oceananigans
using Oceananigans.Architectures

import Base: size, length, eltype, show
import Oceananigans: short_show
import Oceananigans.Architectures: architecture

#####
##### Abstract types
#####

"""
    Center

A type describing the location at the center of a grid cell.
"""
struct Center end

"""
	Face

A type describing the location at the face of a grid cell.
"""
struct Face end

"""
    AbstractTopology

Abstract supertype for grid topologies.
"""
abstract type AbstractTopology end

"""
    Periodic

Grid topology for periodic dimensions.
"""
struct Periodic <: AbstractTopology end

"""
    Bounded

Grid topology for bounded dimensions, e.g., wall-bounded dimensions.
"""
struct Bounded <: AbstractTopology end

"""
    Flat

Grid topology for flat dimensions, generally with one grid point, along which the solution
is uniform and does not vary.
"""
struct Flat <: AbstractTopology end

"""
    Connected

Grid topology for dimensions that are connected to other models or domains on both sides.
"""
const Connected = Periodic  # Right now we just need them to behave like Periodic dimensions except we change the boundary conditions.

"""
    AbstractGrid{FT, TX, TY, TZ}

Abstract supertype for grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractUnderlyingGrid{FT, TX, TY, TZ}

Abstract supertype for "primary" grids (as opposed to grids with immersed boundaries)
with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractUnderlyingGrid{FT, TX, TY, TZ, Arch} <: AbstractGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractRectilinearGrid{FT, TX, TY, TZ}

Abstract supertype for rectilinear grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractRectilinearGrid{FT, TX, TY, TZ, Arch} <: AbstractUnderlyingGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractCurvilinearGrid{FT, TX, TY, TZ}

Abstract supertype for curvilinear grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractCurvilinearGrid{FT, TX, TY, TZ, Arch} <: AbstractUnderlyingGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractHorizontallyCurvilinearGrid{FT, TX, TY, TZ}

Abstract supertype for horizontally-curvilinear grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractHorizontallyCurvilinearGrid{FT, TX, TY, TZ, Arch} <: AbstractCurvilinearGrid{FT, TX, TY, TZ, Arch} end

include("grid_utils.jl")
include("zeros.jl")
include("new_data.jl")
include("automatic_halo_sizing.jl")
include("input_validation.jl")
include("grid_generation.jl")
include("rectilinear_grid.jl")
include("conformal_cubed_sphere_face_grid.jl")
include("latitude_longitude_grid.jl")

end


text:
module Probably

using Base.Iterators: filter, take, repeated
import LinearAlgebra

export State, QBool, QUInt, QBuffer, epr,
  measure, hadamard, not, cnot, swap,
  classical, control

eye(m::AbstractMatrix) = Matrix{eltype(m)}(LinearAlgebra.I, size(m))
a ⊗ b = kron(a, b)

include("operators.jl")
include("constants.jl")
include("state.jl")
include("bits.jl")

include("types/buffer.jl")
include("types/bool.jl")
include("types/int.jl")

end # module


text:
"""
    AbstractDataGraph{T,VL,VD,ED,GD} <: AbstractGraph{T}

General template for graphs with metadata.

- `T<:Integer` is the type of vertices
- `VL` is the type of vertex labels (cannot be a subtype of `Integer`)
- `VD` is the type of vertex data objects
- `ED` is the type of edge data objects
- `GD` is the type of the graph data object
"""
abstract type AbstractDataGraph{T<:Integer,VL,VD,ED,GD} <: AbstractGraph{T} end

Base.copy(g::AbstractDataGraph) = deepcopy(g)

function Base.show(io::IO, g::AbstractDataGraph{T,VL,VD,ED,GD}) where {T,VL,VD,ED,GD}
    n, m = nv(g), ne(g)
    return print(
        io,
        "DataGraph of size ($n, $m) with vertex labels of type $VL, vertex data of type $VD, edge data of type $ED and graph data of type $GD.",
    )
end


text:
"""
Family of nonlinear transformation parametrized by λ that defaults to Box Cox
"""
struct nonlinearTransform{T1<:Function, T2<:Function, T3<:Function, T4<:Function, T5<:Function, T6<:Function, T7<:Function}
    f::T1
    df::T2
    d2f::T3
    df_hyp::T4
    d2f_hyp::T5
    d_mixed::T6
    inv::T7
end

@doc raw"""
    BoxCox(λ::Real)

The Box-Cox power transformation function with parameter λ.

```math
g_\lambda(x) = \begin{cases}
    \frac{x^\lambda - 1}{\lambda} & \lambda \neq 0\\
    \ln(x) & \lambda = 0
\end{cases}
```

External links
* [Power Transform on Wikipedia](https://en.wikipedia.org/wiki/Power_transform)

...
# Arguments
* `x::float64`: evaluation point
* `lambda::float=1.0`: lambda hyperparameter

"""
function boxCox(x, lambda=1)
    lambda == 0 ? Base.log.(x) : (float(x).^lambda.-1)./lambda
end

"""
Derivative of Box-Cox power transformation w.r.t x
"""
function boxCoxPrime(x, lambda=1)
    lambda==0 ? float(x).^(-1) : float(x).^(lambda .-1)
end

"""
Second derivative of Box-Cox power transformation w.r.t x
"""
function boxCoxPrime2(x, lambda=1)
    lambda==0 ? -float(x).^(-2) : (lambda-1)*float(x).^(lambda-2)
end

"""
Derivative of Box-Cox power transformation w.r.t lambda
"""
function boxCoxPrime_lambda(x, lambda=1)
    lambda==0 ? 0 : (lambda * float(x).^lambda .* log.(x) .- float(x).^lambda .+ 1)/lambda^2
end

"""
Second derivative of Box-Cox power transformation w.r.t lambda
"""
function boxCoxPrime_lambda2(x, lambda=1)
    lambda==0 ? 0 : float(x).^lambda .* (Base.log.(x)^2)/lambda -
    (2*x.^lambda .* Base.log.(x))/lambda^2 + 2*(float(x).^lambda-1)/lambda^3
end

"""
Mixed derivative of Box-Cox power transformation w.r.t lambda and x
"""
function boxCoxMixed_lambda_z(x, lambda=1)
    lambda==0 ? 0 : float(x).^(lambda-1) .* Base.log.(x)
end

"""
Inverse Box Cox power transformation
"""
function invBoxCox(y, lambda=1)
    lambda==0 ? Base.exp.(y) : Base.exp.(Base.log.(lambda.*y.+1)./lambda)
end

#define BoxCox Object
boxCoxObj = nonlinearTransform(boxCox, boxCoxPrime, boxCoxPrime2, boxCoxPrime_lambda, boxCoxPrime_lambda2, boxCoxMixed_lambda_z, invBoxCox)

@doc raw"""
    TODO Unimplimented
"""
struct YeoJohnson{T} <: Function end

@doc raw"""
    TODO Unimplimented
"""
struct ArandaOrdaz{T} end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
$(FUNCTIONNAME)

Adds linear objective cut constraint to the `x.relaxed_optimizer`.
"""
function objective_cut!(m::GlobalOptimizer, check_safe::Bool)
    f = m._working_problem._objective_saf
    u = m._global_upper_bound
    if  u < Inf
        b = f.constant
        f.constant = 0.0
        if check_safe && is_safe_cut!(m, f)
            s = LT(u - b + _constraint_tol(m))
            m._affine_objective_cut_ci = MOI.add_constraint(_relaxed_optimizer(m), f, s)
        end
        f.constant = b
        m._new_eval_objective = false
    end
    return
end

"""
    RelaxResultStatus

Status code used internally to determine how to interpret the results from the
solution of a relaxed problem.
"""
@enum(RelaxResultStatus, RRS_OPTIMAL, RRS_DUAL_FEASIBLE, RRS_INFEASIBLE, RRS_INVALID)

"""
$(SIGNATURES)

Takes an `MOI.TerminationStatusCode` and a `MOI.ResultStatusCode` and returns
the tuple `(valid_result::Bool, feasible::Bool)`. The value `valid_result` is
`true` if the pair of codes prove that either the subproblem solution was solved
to global optimality or the subproblem solution is infeasible. The value of
`feasible` is true if the problem is feasible and false if the problem is infeasible.
"""
function relaxed_problem_status(t::MOI.TerminationStatusCode,
                                p::MOI.ResultStatusCode,
                                d::MOI.ResultStatusCode)

    if (t == MOI.OPTIMAL) && (p == MOI.FEASIBLE_POINT)
        return RRS_OPTIMAL
    elseif t == MOI.INFEASIBLE
        if (p == MOI.INFEASIBILITY_CERTIFICATE) ||
           (p == MOI.NO_SOLUTION) || (p == MOI.UNKNOWN_RESULT_STATUS)
            return RRS_INFEASIBLE
        end
    elseif (t == MOI.INFEASIBLE_OR_UNBOUNDED && p == MOI.NO_SOLUTION)
        return RRS_INFEASIBLE
    end
    (d == MOI.FEASIBLE_POINT) && return RRS_DUAL_FEASIBLE
    return RRS_INVALID
end

"""
$(SIGNATURES)

Updates the relaxed constraint by setting the constraint set of `v == x*`` ,
`xL_i <= x_i`, and `x_i <= xU_i` for each such constraint added to the relaxed
optimizer. Resets integral valued constraints to either `EqualTo` or `Interval`
constraints.
"""
function update_relaxed_problem_box!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    d = _relaxed_optimizer(m)
    for i = 1:_variable_num(BranchVar(), m)
        l = _lower_bound(BranchVar(), m, i)
        u = _upper_bound(BranchVar(), m, i)
        v = VI(_bvi(m, i))
        if l == u
            ci_vi_et = MOI.add_constraint(d, v, ET(l))
            push!(m._relaxed_variable_et, (ci_vi_et,i))
        else
            ci_vi_lt = MOI.add_constraint(d, v, LT(u))
            ci_vi_gt = MOI.add_constraint(d, v, GT(l))
            m._node_to_sv_leq_ci[i] = ci_vi_lt
            m._node_to_sv_geq_ci[i] = ci_vi_gt
            push!(m._relaxed_variable_lt, (ci_vi_lt,i))
            push!(m._relaxed_variable_gt, (ci_vi_gt,i))
        end
    end
    return
end

const SOLUTION_EPS = 0.05
function store_lower_solution!(m::GlobalOptimizer{R,S,Q}, d::T) where {R,S,Q<:ExtensionType,T}
    for i = 1:_variable_num(FullVar(), m)
        l = _lower_bound(FullVar(), m, i)
        u = _upper_bound(FullVar(), m, i)
        ladj = l + SOLUTION_EPS*(u - l)
        uadj = u - SOLUTION_EPS*(u - l)
        x = MOI.get(d, MOI.VariablePrimal(), m._relaxed_variable_index[i])
        (x < ladj) && (x = ladj)
        (x > uadj) && (x = uadj)
        m._lower_solution[i] = x
    end
    return
end

function reset_relaxation!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    d = _relaxed_optimizer(m)
    m._cut_iterations = 1
    m._obbt_performed_flag = false
    m._working_problem._relaxed_evaluator.is_first_eval = true
    fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)

    m._new_eval_objective = true
    m._new_eval_constraint = true

    # delete added affine constraints
    foreach(c -> MOI.delete(d, c), m._affine_relax_ci)
    empty!(m._affine_relax_ci)

    # delete variable
    foreach(c -> MOI.delete(d, c[1]), m._relaxed_variable_et)
    foreach(c -> MOI.delete(d, c[1]), m._relaxed_variable_lt)
    foreach(c -> MOI.delete(d, c[1]), m._relaxed_variable_gt)
    foreach(c -> MOI.delete(d, c), m._relaxed_variable_integer)
    empty!(m._relaxed_variable_et)
    empty!(m._relaxed_variable_lt)
    empty!(m._relaxed_variable_gt)
    empty!(m._relaxed_variable_integer)

    # delete objective cut
    !isnothing(m._affine_objective_cut_ci) && MOI.delete(d, m._affine_objective_cut_ci)
    return
end

"""
$(FUNCTIONNAME)

"""
function set_first_relax_point!(m::GlobalOptimizer)
    if m._cut_iterations == 1
        m._working_problem._relaxed_evaluator.is_first_eval = true
        m._new_eval_constraint = true
        m._new_eval_objective = true
        for i = 1:_variable_num(FullVar(), m)
            l = _lower_bound(FullVar(), m, i)
            u = _upper_bound(FullVar(), m, i)
            if isfinite(l) && isfinite(u)
                x = 0.5*(l + u)
            elseif isfinite(l)
                x = min(0.0, u)
            elseif isfinite(u)
                x = max(0.0, l)
            else
                x = 0.0
            end
            _set_lower_solution!(FullVar(), m, x, i)
        end
    end
    return
end

"""
$(TYPEDSIGNATURES)

A routine that adds relaxations for all nonlinear constraints and quadratic constraints
corresponding to the current node to the relaxed problem. This adds an objective cut
(if specified by `objective_cut_on`) and then sets the `_new_eval_constraint` flag
to false indicating that an initial evaluation of the constraints has occurred. If
the `objective_cut_on` flag is `true` then the `_new_eval_objective` flag is also
set to `false` indicating that the objective expression was evaluated.
"""
function relax_all_constraints!(t::ExtensionType, m::GlobalOptimizer, k::Int)
    check_safe = (k == 1) ? false : m._parameters.cut_safe_on
    wp = m._working_problem
    wp._relaxed_evaluator.is_first_eval = m._new_eval_constraint
    foreach(f -> relax!(m, f, k, check_safe), wp._sqf_leq)
    foreach(f -> relax!(m, f, k, check_safe), wp._sqf_eq)
    valid_relax_flag = true
    if valid_relax_flag
        for nl in wp._nonlinear_constr
            valid_relax_flag &= relax!(m, nl, k, check_safe)
        end
    end
    m._new_eval_constraint = false
    (k == 1) && objective_cut!(m, check_safe)
    return valid_relax_flag
end
relax_constraints!(t::ExtensionType, m::GlobalOptimizer, k::Int) = relax_all_constraints!(t, m, k)
relax_constraints!(m::GlobalOptimizer{R,S,Q}, k::Int) where {R,S,Q<:ExtensionType} = relax_constraints!(_ext_typ(m), m, k)

function relax_problem!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    wp = m._working_problem
    if m._cut_iterations == 1
        reset_relaxation!(m)
        if m._nonlinear_evaluator_created
            set_node!(wp._relaxed_evaluator, m._current_node)
            set_reference_point!(m)
            fill!(wp._relaxed_evaluator.subexpressions_eval, false)
        end
        wp._relaxed_evaluator.is_post = m._parameters.subgrad_tighten
        _set_has_value!(wp._objective, false)
        wp._relaxed_evaluator.interval_intersect = false
        update_relaxed_problem_box!(m)
        set_first_relax_point!(m)
    else
        set_reference_point!(m)
    end
    valid_relax_flag = relax_constraints!(m, m._cut_iterations)
    MOI.set(_relaxed_optimizer(m), MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return valid_relax_flag
end

"""
$(SIGNATURES)

Retrieves the lower and upper duals for variable bounds from the
`relaxed_optimizer` and sets the appropriate values in the
`_lower_lvd` and `_lower_uvd` storage fields.
"""
function set_dual!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    d = _relaxed_optimizer(m)
    if MOI.get(d, MOI.DualStatus()) === MOI.FEASIBLE_POINT
        for (c, i) in m._relaxed_variable_lt
            m._lower_uvd[i] = MOI.get(d, MOI.ConstraintDual(), c)
        end
        for (c, i) in m._relaxed_variable_gt
            m._lower_lvd[i] = MOI.get(d, MOI.ConstraintDual(), c)
        end
    else
        fill!(m._lower_lvd, 0.0)
        fill!(m._lower_uvd, 0.0)
    end
    return
end

interval_objective_bound!(m::GlobalOptimizer, f::Nothing) = nothing
function interval_objective_bound!(m::GlobalOptimizer, f::AffineFunctionIneq)
    m._working_problem._relaxed_evaluator.is_first_eval = true
    fL, fU = bound_objective(m)
    if fL > m._lower_objective_value
        m._lower_objective_value = fL
        fill!(m._lower_lvd, 0.0)
        fill!(m._lower_uvd, 0.0)
        m._cut_add_flag = false
    end
end
function interval_objective_bound!(m::GlobalOptimizer, f)
    m._working_problem._relaxed_evaluator.is_first_eval = true
    fL, fU = bound_objective(m)
    fv = _is_input_min(m) ? fL : -fU
    if fv > m._lower_objective_value
        m._lower_objective_value = fv
        fill!(m._lower_lvd, 0.0)
        fill!(m._lower_uvd, 0.0)
        m._cut_add_flag = false
    end
    return
end
interval_objective_bound!(m::GlobalOptimizer) = interval_objective_bound!(m, m._working_problem._objective)

"""
$(SIGNATURES)

Runs interval, linear, quadratic contractor methods followed by obbt and a
constraint programming walk up to tolerances specified in
`EAGO.Optimizer` object.
"""
function preprocess!(t::ExtensionType, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}

    feasible_flag = true
    reset_relaxation!(m)
    if _fbbt_lp_depth(m) >= _iteration_count(m)
        load_fbbt_buffer!(m)
        for _ = 1:_fbbt_lp_repetitions(m)
            ns = NodeBB(_current_node(m))
            for f in m._working_problem._saf_leq
                !(feasible_flag = feasible_flag && fbbt!(m, f)) && break
            end
            !feasible_flag && break
            for f in m._working_problem._saf_eq
                !(feasible_flag = feasible_flag && fbbt!(m, f)) && break
            end
            (same_box(ns,_current_node(m),0.0) || !feasible_flag) && break
        end
        unpack_fbbt_buffer!(m)
    end
    #@show "FBBT Feasibility", feasible_flag

    # done after cp to prevent using cp specific flags in cut generation
    set_first_relax_point!(m)
    # nonlinear CP can detect infeasibility and bound objective even if
    # the relaxation is ill-posed, so one is always used to mitigate numerical issues
    cp_reps = _cp_depth(m) >= _iteration_count(m) ? _cp_repetitions(m) : 1
    for _ = 1:_cp_repetitions(m)
        ns = NodeBB(_current_node(m))
        feasible_flag = feasible_flag && set_constraint_propagation_fbbt!(m)
        (same_box(ns,_current_node(m),0.0) || !feasible_flag) && break
    end
    #@show "CP Feasibility", feasible_flag

    if _obbt_depth(m) >= _iteration_count(m)
        for k = 1:_obbt_repetitions(m)
            ns = NodeBB(_current_node(m))
            feasible_flag = feasible_flag && obbt!(m)
            m._obbt_performed_flag = true
            (same_box(ns,_current_node(m),0.0) || !feasible_flag) && break
        end
    end
    #@show "OBBT Feasibility", feasible_flag
    m._preprocess_feasibility = feasible_flag
    return
end
preprocess!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = preprocess!(_ext_typ(m), m)

"""
$(SIGNATURES)

Checks if a cut should be added and computes a new reference point to add the
cut at. If no cut should be added the constraints not modified in place are
deleted from the relaxed optimizer and the solution is compared with the
interval lower bound. The best lower bound is then used.
"""
function cut_condition(t::ExtensionType, m::GlobalOptimizer)
    obj_old = m._last_cut_objective
    obj_new = m._lower_objective_value
    flag = m._cut_iterations < _cut_max_iterations(m)
    flag &= obj_new - obj_old > _cut_ϵ_abs(m)
    flag &= obj_new - obj_old > _cut_ϵ_rel(m)*abs(obj_new)
    return flag
end
cut_condition(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = cut_condition(_ext_typ(m), m)

is_integer_subproblem(m) = !continuous(_current_node(m))
"""
$(SIGNATURES)

Constructs and solves the relaxation using the default EAGO relaxation scheme
and optimizer on node `y`.
"""
function lower_problem!(t::ExtensionType, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}

    d = _relaxed_optimizer(m)
    m._last_cut_objective = typemin(Float64)
    m._lower_objective_value = typemin(Float64)

    t_status = MOI.OPTIMIZE_NOT_CALLED
    p_status = MOI.OTHER_RESULT_STATUS
    d_status = MOI.OTHER_RESULT_STATUS
    status = RRS_INVALID

    set_first_relax_point!(m)
    MOI.set(d, MOI.ObjectiveFunction{SAF}(), m._working_problem._objective_saf)
    while true
        relax_problem!(m)
        m._last_cut_objective = m._lower_objective_value
        MOI.optimize!(d)
        t_status = MOI.get(d, MOI.TerminationStatus())
        p_status = MOI.get(d, MOI.PrimalStatus())
        d_status = MOI.get(d, MOI.DualStatus())
        status = relaxed_problem_status(t_status, p_status, d_status)
        if status != RRS_OPTIMAL
            break
        end
        m._lower_objective_value = MOI.get(d, MOI.ObjectiveValue())
        #@show m._lower_objective_value
        if cut_condition(m)
            store_lower_solution!(m, d)
            m._cut_iterations += 1
        else
            break
        end
    end

    # activate integrality conditions for MIP & solve MIP subproblem
    if is_integer_subproblem(m) && (status !== RRS_INFEASIBLE)
        m._last_cut_objective = m._lower_objective_value
        for i = 1:_variable_num(BranchVar(), m)
            l = _lower_bound(BranchVar(), m, i)
            u = _upper_bound(BranchVar(), m, i)
            if is_integer(BranchVar(), m, i) && (l != u)
                c_integer = MOI.add_constraint(d, VI(_bvi(m, i)), MOI.Integer())
                push!(m._relaxed_variable_integer, c_integer)
            end
        end
        MOI.optimize!(d)
        t_status = MOI.get(d, MOI.TerminationStatus())
        p_status = MOI.get(d, MOI.PrimalStatus())
        d_status = MOI.get(d, MOI.DualStatus())
        status = relaxed_problem_status(t_status, p_status, d_status)
    end

    # check status -- if not feasible/infeasible then fallback to interval bounds
    if status == RRS_OPTIMAL
        m._lower_objective_value = MOI.get(d, MOI.ObjectiveValue())
    end
    #@show m._lower_objective_value

    m._lower_termination_status = t_status
    m._lower_primal_status = p_status
    m._lower_dual_status = d_status
    status = relaxed_problem_status(t_status, p_status, d_status)
    if status == RRS_INFEASIBLE
        m._lower_feasibility  = false
        m._lower_objective_value = -Inf
        return
    end

    # set dual values
    set_dual!(m)
    m._lower_feasibility = true
    store_lower_solution!(m, d)
    if status == RRS_DUAL_FEASIBLE
        m._lower_objective_value = MOI.get(d, MOI.DualObjectiveValue())
    end
    #@show m._lower_objective_value
    interval_objective_bound!(m)
    return
end
lower_problem!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = lower_problem!(_ext_typ(m), m)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Testing the hypothesis that evolution under a varying environment increases evolvability.

# The purpose of the function run_run_pop_evolvability() is to explore settings of popsize,
# mutrate, and ngens for the population-based evolutionary algorithm.  The function allows
# for different settings of popsize and mutrate via range parameters.  The parameter
# ngens_popsize is the product of ngens and popsize.  Instead of using ngens as a parameter,
# ngens is computed by the formula  ngens = Int(round(ngens_popsize/popsize))  which
# specifies a number of function evaluations for the run.
function run_run_pop_evolvability( nreps::Int64, p::Parameters, popsize_rng::CGP.IntRange, gl::GoalList, ngens_popsize::Int64,
      start_phmut_gen::Int64, phmut_interval::Int64, mutrate_rng::CGP.FloatRange;
    csvfile::String="" )
  sumdf = DataFrame()
  sumdf.gen = Int64[]
  sumdf.popsize = Int64[]
  sumdf.mutrate = Float64[]
  sumdf.max_fitness = Float64[]
  sumdf.evolvability = Float64[]
  sumdf.robustness= Float64[]
  run_tuples = Tuple[]
  for popsize in popsize_rng
    ngens = Int(round(ngens_popsize/popsize))
    for mutrate in mutrate_rng
      push!(run_tuples, (popsize,ngens,mutrate))
    end
  end
  results=pmap(x->run_pop_evolvability(nreps,p,x[1],gl,x[2],start_phmut_gen, phmut_interval,x[3],use_pmap=false),run_tuples)
  for res in results
    push!(sumdf,res)
  end
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# goallist: ",gl)
      println(f,"# nreps: ",nreps)
      println(f,"# start_phmut_gen: ",start_phmut_gen)
      println(f,"# phmut_interval: ",phmut_interval)
      CSV.write( f, sumdf, append=true, writeheader=true )
    end
  end
  sumdf
end

# Averages results over reps
# use_pmap must be false when called by run_run_pop_evolvability
# If length(gl) == 0, generate goallist within each call to pop_evolvability
function run_pop_evolvability( nreps::Int64, p::Parameters, popsize::Int64, gl::GoalList, ngens::Int64,
    start_phmut_gen::Int64, phmut_interval::Int64,  mutrate::Float64=1.0; numgoals::Int64=1,
    csvfile::String="", use_pmap::Bool=false, all_gens::Bool=false )
  #@assert p.numoutputs == length(gl[1])
  println("mutrate: ",mutrate,"  gl: ",gl,"  length(gl): ",length(gl))
  avgdf = DataFrame()
  avgdf.gen = collect(1:ngens)
  avgdf.popsize = fill(popsize,ngens)
  avgdf.mutrate = fill(mutrate,ngens)
  avgdf.max_fitness = zeros(Float64,ngens)
  avgdf.evolvability = zeros(Float64,ngens)
  avgdf.robustness = zeros(Float64,ngens)
  if use_pmap
    dflist = pmap(x->pop_evolvability( p, popsize, gl, ngens, start_phmut_gen, phmut_interval, mutrate,numgoals=numgoals), collect(1:nreps))
  else
    dflist = map(x->pop_evolvability( p, popsize, gl, ngens, start_phmut_gen, phmut_interval, mutrate,numgoals=numgoals), collect(1:nreps))
  end
  for df in dflist
    avgdf.max_fitness += df.max_fitness
    avgdf.evolvability += df.evolvability
    avgdf.robustness += df.robustness
  end
  avgdf.max_fitness /= nreps
  avgdf.evolvability /= nreps
  avgdf.robustness /= nreps
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# goallist: ",gl)
      println(f,"# nreps: ",nreps)
      CSV.write( f, avgdf, append=true, writeheader=true )
    end
  end
  if all_gens
    avgdf
  else
    avgdf[end,:]
  end
end

# If length(gl)==0, generate goallist within each call to pop_evolvability
function pop_evolvability( p::Parameters, popsize::Int64, gl::GoalList, ngens::Int64,
    start_phmut_gen::Int64, phmut_interval::Int64,  mutrate::Float64=1.0;
    csvfile::String="", numgoals::Int64=1 )
  if length(gl) > 0
    glmut = deepcopy(gl)
  else
    glmut = randgoallist( numgoals, p )
  end
  #println("glmut: ",glmut)
  df = DataFrame()
  df.gen = Int64[]
  df.max_fitness = Float64[]
  df.evolvability = Float64[]
  df.robustness = Float64[]
  funcs = default_funcs( p.numinputs )
  pop = [ random_chromosome(p,funcs) for i = 1: popsize ]
  count_mutations = 0
  fitness_vector = zeros(Float64,popsize)
  maxfit = 0.0
  for i in 1:popsize
    c = pop[i]
    c.fitness = fitness_funct( p, c, glmut )
    fitness_vector[i] = pop[i].fitness
    maxfit = fitness_vector[i] > maxfit ? fitness_vector[i] : maxfit
  end
  ( evolvability, robustness ) = pop_evolvability_robustness( pop )
  push!( df.gen, 1 )
  push!( df.max_fitness, maxfit )
  push!( df.evolvability, evolvability )
  push!( df.robustness, robustness )
  for gen = 2:ngens
    if gen >= start_phmut_gen && gen % phmut_interval == 0
      phmut!( glmut, p )
      println("phmut gen: ",gen,"  glmut: ",glmut)
    end
    for i in 1:popsize
      c = pop[i]
      if rand() <= mutrate
        pop[i] = c = mutate_chromosome!(deepcopy(c),funcs)[1]
        c.fitness = fitness_funct( p, c, glmut )
        #println("mutation count: ",count_mutations,"  phenotype: ",@sprintf("0x%x",(output_values(c)[1])),"  fitness: ",c.fitness)
        count_mutations += 1
      end
    end
    fitness_vector = ([ pop[i].fitness for i = 1:popsize ])
    maxfit = findmax(fitness_vector)[1]
    propsel!( pop, fitness_vector, maxfit=maxfit )
    ( evolvability, robustness ) = pop_evolvability_robustness( pop )
    push!( df.gen, gen )
    push!( df.max_fitness, maxfit )
    push!( df.evolvability, evolvability )
    push!( df.robustness, robustness )
  end
  println("glmut: ",glmut,"  mutrate: ",mutrate,"  max_fit: ",maxfit)
  df
end

function phmut!( glmut::GoalList, p::Parameters )
  gl_component = length(glmut) > 1 ? rand(1:length(glmut)) : 1
  goal_index = length(glmut[1]) > 1 ? rand(1:length(glmut[1])) : 1
  max_shift = 2^p.numinputs - 1
  shift = rand(0:max_shift)
  println("gl_component: ",gl_component,"  goal_index: ",goal_index,"  shift: ",shift)
  mut_mask = MyInt(1) << shift
  #@printf("mut_mask: 0x%x\n",mut_mask)
  # A more elegant version of the next line is:
  # glmut[gl_component][goal_index] ⊻= mut_mask
  glmut[gl_component][goal_index] = xor(glmut[gl_component][goal_index], mut_mask)
  return glmut
end

function fitness_funct( p::Parameters, c::Chromosome, gl::GoalList )
  maximum( 1.0-hamming_distance( g, output_values(c), p.numinputs ) for g in gl )
end


text:
using Images
using ImageView
using MAT

vars = matread("../../DATA/allFaces.mat")

n = Int(vars["n"])
m = Int(vars["m"])
faces = vars["faces"]
nfaces = Int.(vars["nfaces"])

allPersons = zeros(n*6, m*6)

count = 1
for i = 1:6
    for j = 1:6
        allPersons[1+(i-1)*n:i*n, 1+(j-1)*m:j*m] = reshape(faces[:,1+sum(nfaces[1:count-1])], n, m)
        global count += 1
    end
end
imshow(allPersons)

for person = 1:length(nfaces)
    subset = faces[:, 1+sum(nfaces[1:person-1]):sum(nfaces[1:person])]
    allFaces = zeros(n*8, m*8);

    global count = 1
    for i = 1:8
        for j = 1:8
            if (count <= nfaces[person])
                allFaces[1+(i-1)*n:i*n, 1+(j-1)*m:j*m] = reshape(subset[:,count], n, m)
                global count += 1
            end
        end
    end
    imshow(allFaces)
end


text:
# ---
# title: 21. Merge Two Sorted Lists
# id: problem21
# author: zhwang
# date: 2022-01-20
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/merge-two-sorted-lists/description/>
# hidden: true
# ---
#
# Merge two sorted linked lists and return it as a new **sorted** list. The new
# list should be made by splicing together the nodes of the first two lists.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/03/merge_ex1.jpg)
#
#
#
#     Input: l1 = [1,2,4], l2 = [1,3,4]
#     Output: [1,1,2,3,4,4]
#
#
# **Example 2:**
#
#
#
#     Input: l1 = [], l2 = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: l1 = [], l2 = [0]
#     Output: [0]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in both lists is in the range `[0, 50]`.
#   * `-100 <= Node.val <= 100`
#   * Both `l1` and `l2` are sorted in **non-decreasing** order.
#
#
## @lc code=start
using LeetCode

function merge_two_lists(
    list1::Union{ListNode,Nothing}, list2::Union{ListNode,Nothing}
)::Union{ListNode,Nothing}
    isnothing(list1) && isnothing(list2) && return nothing
    node = fake_head = ListNode()
    while !isnothing(list1) && !isnothing(list2)
        node.next = ListNode()
        node = next(node)
        if val(list1) < val(list2)
            node.val = list1.val
            list1 = next(list1)
        else
            node.val = list2.val
            list2 = next(list2)
        end
    end
    l = isnothing(list2) ? list1 : list2
    while !isnothing(l)
        node.next = ListNode()
        node = next(node)
        node.val = l.val
        l = next(l)
    end
    return next(fake_head)
end

## @lc code=end


text:
#
# Miscellaneous functions
#

# Redefine unsafe_load, unsafe_store!, read, and write so that they pack the type
function define_packed(ty::DataType)
    @assert isbitstype(ty)
    packed_offsets = cumsum([jlsizeof(x) for x in ty.types])
    sz = pop!(packed_offsets)
    pushfirst!(packed_offsets, 0)

    if sz != jlsizeof(ty)
        @eval begin
            function jlunsafe_store!(p::Ptr{$ty}, x::$ty)
                $([:(jlunsafe_store!(convert(Ptr{$(ty.types[i])}, p+$(packed_offsets[i])), getfield(x, $i)))
                   for i = 1:length(packed_offsets)]...)
            end
            function jlunsafe_load(p::Ptr{$ty})
                $(Expr(:new, ty, [:(jlunsafe_load(convert(Ptr{$(ty.types[i])}, p+$(packed_offsets[i]))))
                                   for i = 1:length(packed_offsets)]...))
            end
            jlsizeof(::Union{$ty,Type{$ty}}) = $sz
        end
    end

    @eval begin
        @inline jlwrite(io::Union{MmapIO,BufferedWriter}, x::$ty) = _write(io, x)
        @inline jlread(io::Union{MmapIO,BufferedReader}, x::Type{$ty}) = _read(io, x)
        function jlread(io::IO, ::Type{$ty})
            $(Expr(:new, ty, [:(jlread(io, $(ty.types[i]))) for i = 1:length(packed_offsets)]...))
        end
        function jlwrite(io::IO, x::$ty)
            $([:(jlwrite(io, getfield(x, $i))) for i = 1:length(packed_offsets)]...)
            nothing
        end
    end
    nothing
end

"""
    read_size(io::IO, flags::UInt8)

Loads a variable-length size according to flags
Expects that the first two bits of flags mean:
0   The size of the Length of Link Name field is 1 byte.
1   The size of the Length of Link Name field is 2 bytes.
2   The size of the Length of Link Name field is 4 bytes.
3   The size of the Length of Link Name field is 8 bytes.
Returns the size as an Int
"""
function read_size(io::IO, flags::UInt8)
    if (flags & 3) == 0
        Int(jlread(io, UInt8))
    elseif (flags & 3) == 1
        Int(jlread(io, UInt16))
    elseif (flags & 3) == 2
        Int(jlread(io, UInt32))
    else
        Int(jlread(io, UInt64))
    end
end

# Determine what the size flag should be
# Same rules as above
function size_flag(sz::Integer)
    if sz <= typemax(UInt8)
        UInt8(0)
    elseif sz <= typemax(UInt16)
        UInt8(1)
    elseif sz <= typemax(UInt32)
        UInt8(2)
    else
        UInt8(3)
    end
end

# Store a size
function write_size(io::IO, sz::Integer)
    if sz <= typemax(UInt8)
        jlwrite(io, UInt8(sz))
    elseif sz <= typemax(UInt16)
        jlwrite(io, UInt16(sz))
    elseif sz <= typemax(UInt32)
        jlwrite(io, UInt32(sz))
    else
        jlwrite(io, UInt64(sz))
    end
end

# Get the size of the size
function size_size(sz::Integer)
    if sz <= typemax(UInt8)
        1
    elseif sz <= typemax(UInt16)
        2
    elseif sz <= typemax(UInt32)
        4
    else
        8
    end
end

"""
    symbol_length(x::Symbol)

Returns the length of the string represented by `x`.
"""
symbol_length(x::Symbol) = ccall(:strlen, Int, (Cstring,), x)


text:
export GuildEmbed

"""
A [`Guild`](@ref) embed.
More details [here](https://discordapp.com/developers/docs/resources/guild#guild-embed-object).
"""
struct GuildEmbed
    enabled::Bool
    channel_id::Nullable{Snowflake}
end
@boilerplate GuildEmbed :constructors :docs :lower :merge :mock


text:
using SparseArrays, LinearAlgebra

include("read.jl")
include("write.jl")


text:
using PolygonArea
using Test

@testset "PolygonArea.jl" begin
    include("test_points.jl")
    include("test_halfplanes.jl")
    include("test_polygons.jl")
end


text:
"""
```
solve(m::AbstractModel)
```

Driver to compute the model solution and augment transition matrices.

### Arguments
- `m`: the model object

### Return values
 - TTT, RRR, and CCC matrices of the state transition equation:
    ```
    S_t = TTT*S_{t-1} + RRR*ϵ_t + CCC
    ```
"""
function solve(m::AbstractModel)

    # Get equilibrium condition matrices
    Γ0, Γ1, C, Ψ, Π  = eqcond(m)

    # Solve model
    TTT_gensys, CCC_gensys, RRR_gensys, fmat, fwt, ywt, gev, eu, loose = gensys(Γ0, Γ1, C, Ψ, Π, 1+1e-6)
    if !((eu[1] == 1) & (eu[2] == 1))
        throw(GensysError("Gensys does not give existence"))
    end
    TTT_gensys = real(TTT_gensys)
    RRR_gensys = real(RRR_gensys)
    CCC_gensys = reshape(CCC_gensys, length(CCC_gensys), 1)

    # Augment states
    TTT, RRR, CCC = augment_states(m, TTT_gensys, RRR_gensys, CCC_gensys)

    return TTT, RRR, CCC
end

"""
```
GensysError <: Exception
```
A `GensysError` is thrown when Gensys does not give a unique solution, or no solution
exists. If a `GensysError`is thrown during Metropolis-Hastings, it is caught by `posterior`.
`posterior` then returns a value of `-Inf`, which Metropolis-Hastings always rejects.
### Fields
* `msg::ASCIIString`: Info message. Default = "Error in gensys."
"""
type GensysError <: Exception
    msg::ASCIIString
end
GensysError() = GensysError("Error in gensys.")
Base.showerror(io::IO, ex::GensysError) = print(io, ex.msg)


text:
function calFLDOPartition!(
    ocn :: Ocean
)

    @loop_hor ocn i j let

        FLDO = ocn.FLDO[i, j]

        if FLDO == -1
            ocn.FLDO_ratio_top[i, j] = 0.0
            ocn.FLDO_ratio_bot[i, j] = 1.0
        else
            ocn.FLDO_ratio_top[i, j] = (ocn.h_ML[i, j] + ocn.zs[FLDO, i, j]) / ocn.hs[FLDO, i, j]
            ocn.FLDO_ratio_bot[i, j] = 1.0 - ocn.FLDO_ratio_top[i, j]
        end
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM