text:
using AnalogDiscovery

function caps(devIndex=0)
  cDev = enumDevices()
  @show ("Available",cDev)
  myDev = Int32(devIndex)
  if cDev > 0
    if enumDeviceIsOpened(myDev)
      error("Error opening device, already in use.")
    end
    hdwf = deviceOpen(myDev)
    if hdwf == Int32(0)
      error("Failed to open the requested device.")
    end
  else
    error("No devices available.")
  end

  bmin,bmax = analogInBufferSizeInfo(hdwf)
  show(("Buffer Size, Min:",bmin,"Max:",bmax))

  chVoltMin,chVoltMax,chSteps = analogInChannelRangeInfo(hdwf)
  show(("chVoltMin:",chVoltMin,"chVoltMax:",chVoltMax,"chSteps:",chSteps))
  chCount = analogInChannelCount(hdwf)
  show(("chCount:",chCount))

  deviceClose(hdwf)
  nothing
end

caps()


text:
@testset "Cut Type" begin
    p = Pattern(strengthened_benders=1, lagrangian=0)
    @test SDDiP.getcuttype(1, p) == :strengthened_benders
    @test SDDiP.getcuttype(2, p) == :strengthened_benders

    p = Pattern(benders=1, lagrangian=0)
    @test SDDiP.getcuttype(1, p) == :benders
    @test SDDiP.getcuttype(2, p) == :benders

    p = Pattern(benders=1, strengthened_benders=0, lagrangian=1)
    @test SDDiP.getcuttype(1, p) == :lagrangian
    @test SDDiP.getcuttype(2, p) == :benders
    @test SDDiP.getcuttype(3, p) == :lagrangian
    @test SDDiP.getcuttype(4, p) == :benders

    p = Pattern(benders=2, strengthened_benders=1, lagrangian=0)
    @test SDDiP.getcuttype(1, p) == :benders
    @test SDDiP.getcuttype(2, p) == :strengthened_benders
    @test SDDiP.getcuttype(3, p) == :benders
    @test SDDiP.getcuttype(4, p) == :benders
    @test SDDiP.getcuttype(5, p) == :strengthened_benders

    p = Pattern(benders=1, strengthened_benders=1, lagrangian=1)
    @test SDDiP.getcuttype(1, p) == :strengthened_benders
    @test SDDiP.getcuttype(2, p) == :lagrangian
    @test SDDiP.getcuttype(3, p) == :benders
    @test SDDiP.getcuttype(4, p) == :strengthened_benders
end

@testset "Binary Expansion" begin
    @test_throws Exception binexpand(0)
    @test binexpand(1) == [1]
    @test binexpand(2) == [0, 1]
    @test binexpand(3) == [1, 1]
    @test binexpand(4) == [0, 0, 1]
    @test binexpand(5) == [1, 0, 1]
    @test binexpand(6) == [0, 1, 1]
    @test binexpand(7) == [1, 1, 1]

    @test binexpand(typemax(Int)) == ones(Int, SDDiP._2i_L)

    @test 0 == bincontract([0])
    @test 1 == bincontract([1])
    @test 0 == bincontract([0, 0])
    @test 1 == bincontract([1, 0])
    @test 2 == bincontract([0, 1])
    @test 3 == bincontract([1, 1])
    @test 2 == bincontract([0, 1, 0])
    @test 3 == bincontract([1, 1, 0])
    @test 4 == bincontract([0, 0, 1])
    @test 5 == bincontract([1, 0, 1])
    @test 6 == bincontract([0, 1, 1])
    @test 7 == bincontract([1, 1, 1])
    @test typemax(Int) == bincontract(ones(Int, SDDiP._2i_L))

    @test binexpand(1, length=3) == [1, 0, 0]
    @test binexpand(2, length=3) == [0, 1, 0]
    @test binexpand(4, length=3) == [0, 0, 1]
    @test_throws Exception binexpand(8, length=3)

    @test binexpand(1, maximum=7) == [1, 0, 0]
    @test binexpand(2, maximum=7) == [0, 1, 0]
    @test binexpand(4, maximum=7) == [0, 0, 1]
    @test binexpand(4, maximum=7, length=2) == [0, 0, 1]
    @test_throws Exception binexpand(8, maximum=7)

    @test_throws(Exception, bincontract(Float64, [1, 0], -0.1))
    @test isapprox(bincontract(Float64, [0],       0.1), 0.0, atol=1e-4)
    @test isapprox(bincontract(Float64, [1],       0.1), 0.1, atol=1e-4)
    @test isapprox(bincontract(Float64, [0, 1],    0.1), 0.2, atol=1e-4)
    @test isapprox(bincontract(Float64, [1, 1],    0.1), 0.3, atol=1e-4)
    @test isapprox(bincontract(Float64, [0, 1, 0], 0.1), 0.2, atol=1e-4)
    @test isapprox(bincontract(Float64, [1, 1, 0], 0.1), 0.3, atol=1e-4)
    @test isapprox(bincontract(Float64, [1, 0, 1], 0.1), 0.5, atol=1e-4)
    @test binexpand(0.5)        == binexpand(5)
    @test binexpand(0.54)       == binexpand(5)
    @test binexpand(0.56, 0.1)  == binexpand(6)
    @test binexpand(0.5, 0.01)  == binexpand(50)
    @test binexpand(0.54, 0.01) == binexpand(54)
    @test binexpand(0.56, 0.01) == binexpand(56)

    @test binexpand(0.5, length=5)        == binexpand(5, length=5)
    @test binexpand(0.56, 0.01, length=8) == binexpand(56, length=8)

    @test binexpand(0.5, maximum=3.1)         == binexpand(5, length=5)
    @test binexpand(0.56, 0.01, maximum=2.55) == binexpand(56, length=8)
end

@testset "SDDiP Examples" begin
for example in [
        "SDDiPnewsvendor.jl",
        "all_blacks.jl",
        "generation_expansion.jl",
        "stochastic_all_blacks.jl",
        "booking_management.jl",
        "airconditioning.jl",
        "vehicle_location.jl",
        "integer_optimality.jl"
    ]
    @testset "$example" begin
        include(joinpath(examples_dir, example))
    end
end
end

ext: jl


text:
module NurseAllocation

using JuMP
using Gurobi
using LinearAlgebra
using MathOptInterface

export nurse_allocation

function nurse_allocation(
		initial_nurses::Array{<:Real,1},
		demand::Array{<:Real,2},
		adj_matrix::BitArray{2};

		sent_penalty::Real=0,
		smoothness_penalty::Real=0,

		no_artificial_shortage::Bool=false,
		no_worse_shortage::Bool=false,
		fully_connected::Bool=false,

		sendreceive_switch_time::Int=0,
		min_send_amt::Real=0,
		setup_cost::Real=0,

		verbose::Bool=false,
)

	###############
	#### Setup ####
	###############

	N, T = size(demand)
	@assert(size(initial_nurses) == (N,))
	@assert(size(adj_matrix) == (N,N))

	###############
	#### Model ####
	###############

	model = Model(Gurobi.Optimizer)
	if !verbose set_silent(model) end

	###############
	## Variables ##
	###############

	@variable(model, sent[1:N,1:N,1:T] >= 0)
	@variable(model, obj_dummy[1:N,1:T] >= 0)

	#################
	## Expressions ##
	#################

	# compute active nurses
	@expression(model, active_nurses[i=1:N,t=0:T],
		initial_nurses[i]
		- sum(sent[i,:,1:t])
		+ sum(sent[:,i,1:t])
	)

	objective = @expression(model, sum(obj_dummy))

	######################
	## Hard Constraints ##
	######################

	# sent nurses ≦ active nurses
	@constraint(model, [i=1:N,t=1:T], sum(sent[i,:,t]) <= active_nurses[i,t-1])

	# objective
	@constraint(model, [i=1:N,t=1:T], obj_dummy[i,t] >= demand[i,t] - active_nurses[i,t])

	for i in 1:N, t in 1:T
		fix(sent[i,i,t], 0, force=true)
	end
	if !fully_connected
		for i in 1:N, j in 1:N
			if ~adj_matrix[i,j]
				for t in 1:T
					fix(sent[i,j,t], 0, force=true)
				end
			end
		end
	end

	##########################
	## Optional Constraints ##
	##########################

	if no_artificial_shortage
		for i in 1:N, t in 1:T
			if initial_nurses[i] >= demand[i,t]
				@constraint(model, active_nurses[i,t] >= demand[i,t])
			end
		end
	end

	if no_worse_shortage
		for i in 1:N, t in 1:T
			if initial_nurses[i] < demand[i,t]
				@constraint(model, active_nurses[i,t] >= initial_nurses[i])
			end
		end
	end

	if min_send_amt <= 0
		@constraint(model, sent .>= 0)
	else
		@constraint(model, [i=1:N,j=1:N,t=1:T], sent[i,j,t] in MOI.Semicontinuous(Float64(min_send_amt), Inf))
	end

	if sendreceive_switch_time > 0
		@constraint(model, [i=1:N,t=1:T-1],
			[sum(sent[:,i,t]), sum(sent[i,:,t:min(t+sendreceive_switch_time,T)])] in MOI.SOS1([1.0, 1.0])
		)
		@constraint(model, [i=1:N,t=1:T-1],
			[sum(sent[:,i,t:min(t+sendreceive_switch_time,T)]), sum(sent[i,:,t])] in MOI.SOS1([1.0, 1.0])
		)
	end

	########################
	## Optional Penalties ##
	########################

	if sent_penalty > 0
		add_to_expression!(objective, sent_penalty*sum(sent))
	end

	if smoothness_penalty > 0
		@variable(model, smoothness_dummy[i=1:N,j=1:N,t=1:T-1] >= 0)
		@constraint(model, [t=1:T-1],  (sent[:,:,t] - sent[:,:,t+1]) .<= smoothness_dummy[:,:,t])
		@constraint(model, [t=1:T-1], -(sent[:,:,t] - sent[:,:,t+1]) .<= smoothness_dummy[:,:,t])

		add_to_expression!(objective, smoothness_penalty * sum(smoothness_dummy))
		add_to_expression!(objective, smoothness_penalty * sum(sent[:,:,1]))
	end

	if setup_cost > 0
		@variable(model, setup_dummy[i=1:N,j=i+1:N], Bin)
		@constraint(model, [i=1:N,j=i+1:N], [1-setup_dummy[i,j], sum(sent[i,j,:])+sum(sent[j,i,:])] in MOI.SOS1([1.0, 1.0]))
		add_to_expression!(objective, setup_cost*sum(setup_dummy))
	end

	###############
	#### Solve ####
	###############

	@objective(model, Min, objective)
	optimize!(model)

	return model
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
type QHistory{I,V} <: UnivalueHistory{I}
    lastiter::I
    storage::Deque{Tuple{I,V}}

    function QHistory(::Type{V}, ::Type{I})
        new(typemin(I), Deque{Tuple{I,V}}())
    end
end

function QHistory{I,V}(v::Type{V}, i::Type{I} = Int)
    QHistory{I,V}(v, i)
end

# ==========================================================================
#

Base.length(history::QHistory) = length(history.storage)
Base.enumerate(history::QHistory) = history.storage
Base.first(history::QHistory) = front(history.storage)
Base.last(history::QHistory) = back(history.storage)

function Base.push!{I,V}(
        history::QHistory{I,V},
        iteration::I,
        value::V)
    lastiter = history.lastiter
    iteration > lastiter || throw(ArgumentError("Iterations must increase over time"))
    history.lastiter = iteration
    push!(history.storage, (iteration, value))
    value
end

function Base.push!{I,V}(
        history::QHistory{I,V},
        value::V)
    lastiter = history.lastiter == typemin(I) ? zero(I) : history.lastiter
    iteration = lastiter + one(history.lastiter)
    history.lastiter = iteration
    push!(history.storage, (iteration, value))
    value
end

function Base.get{I,V}(history::QHistory{I,V})
    l = length(history)
    k, v = front(history.storage)
    karray = zeros(I, l)
    varray = Array(V, l)
    i = 1
    for (k, v) in enumerate(history)
        karray[i] = k
        varray[i] = v
        i += 1
    end
    karray, varray
end

Base.print{I,V}(io::IO, history::QHistory{I,V}) = print(io, "$(length(history)) elements {$I,$V}")

function Base.show{I,V}(io::IO, history::QHistory{I,V})
    println(io, "QHistory")
    println(io, "    types: $I, $V")
    print(io,   "    length: $(length(history))")
end


text:
using Optim

function readAndProcCoinData(filename)
  data = readdlm(filename)
  ndata = length(data)
  nheads = int(zeros(ndata))
  procdata = int(zeros(ndata))
  for i=1:ndata
    nheads[i] = length(matchall(r"H", data[i]))
    procdata[i] = binomial(50, nheads[i])
  end
  return nheads, procdata
end

function sigmoid(x)
  return 1.0./(1.0.+exp(-x))
end

function bin50(k, p)
  # d here is (50 choose d)
  # println(exp(log(d) + d*log(p)  ))
  return procdata[k]*(p^nheads[k])*((1.0-p)^(50-nheads[k]))
end

function loglik(x)
  dists = Any[]
  theta = sigmoid(x[1:4])
  cs = sigmoid(x[5:end])
  ntheta = length(theta)

  cs = reshape(cs, nprocdata, 4)
  cs = cs ./ sum(cs,2)

  ll = 0
  for k=1:length(procdata)
    for i=1:ntheta
      ll += bin50(k, theta[i])*cs[k,i]
    end
  end
  ll *= -1
end

function getMeSomeMLE(procdata, theta)
  optimize(loglik, theta, show_trace=true, autodiff=true,
            method = :cg, ftol = 1e-12, iterations=5000)
end

# read in and process the data
nheads, procdata = readAndProcCoinData("coindata.txt")
nprocdata = length(procdata)

# optimize and return the maximum likelihood estimates
theta = rand(4)
cs = rand(nprocdata*4)
x = [theta; cs]
results = getMeSomeMLE(procdata, x)
x_opt = results.minimum
theta_opt = x_opt[1:4]
cs_opt = x_opt[5:end]
println("Maximum Likelihood Estimates: ", sigmoid(theta_opt))

max_stars_repo_name: haroldsoh/puzzles

max_issues_repo_name: haroldsoh/puzzles

max_forks_repo_name: haroldsoh/puzzles


text:
module TestIndexing

using Test
using RLEVectors

@testset begin

# ind2run
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),5) == 3
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),6) == 3
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),7) == 4
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),7:9) == 4:5
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),[3, 5, 6]) == [2, 3, 3]

# ind2runcontext
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),5) == (3,1,1)
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),6) == (3,2,0)
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),7) == (4,1,1)
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),7:9) == (4,5,1,1)

## getindex for single position
x = RLEVector([1,2,3,4],[2,4,6,8])
@test x[1] == 1
@test x[3] == 2
@test x[8] == 4
@test_throws BoundsError x[9]

## getindex on multiple positions
x = RLEVector([1,2,3,4],[2,4,6,8])
@test x[ [1,2,3] ] == RLEVector([1,1,2])
@test x[ [8,2,4] ] == RLEVector([4,1,2])
@test x[ 2:5 ] == RLEVector([1,2,2,3])
@test x[ 5:-1:2 ] == RLEVector([3,2,2,1])
@test x[ 1:end ] == x

## getindex with logical
x = RLEVector([1,2],[2,4])
@test x[ [true, false, false, true] ] == RLEVector([1, 2])

## setindex! for single position
y = RLEVector([3,4,4,5,5])
y[1] = 20
@test collect(y) == [20,4,4,5,5]

# No change
x = RLEVector([1,2,3,4],[2,4,6,8])
x[4] = 2
@test collect(x) == [1,1,2,2,3,3,4,4]
@test x.runvalues == [1,2,3,4]
@test x.runends == [2,4,6,8]

# left of run
x = RLEVector([1,2,3,4],[2,4,6,8])
x[1] = 9
@test collect(x) == [9,1,2,2,3,3,4,4]
@test x.runvalues == [9,1,2,3,4]
@test x.runends == [1,2,4,6,8]

# right of run
x = RLEVector([1,2,3,4],[2,4,6,8])
x[2] = 7
@test collect(x) == [1,7,2,2,3,3,4,4]
@test x.runvalues == [1,7,2,3,4]
@test x.runends == [1,2,4,6,8]

# left of run making repeat
x = RLEVector([1,2,3,4],[2,4,6,8])
x[3] = 1
@test collect(x) == [1,1,1,2,3,3,4,4]
@test x.runvalues == [1,2,3,4]
@test x.runends == [3,4,6,8]

# right of run making repeat
x = RLEVector([1,2,3,4],[2,4,6,8])
x[4] = 3
@test collect(x) == [1,1,2,3,3,3,4,4]
@test x.runvalues == [1,2,3,4]
@test x.runends == [2,3,6,8]

# middle of run
x = RLEVector([1,2,3],[3,6,9])
x[5] = 5
@test collect(x) == [1,1,1,2,5,2,3,3,3]
@test x.runvalues == [1,2,5,2,3]
@test x.runends == [3,4,5,6,9]

# ends
x = RLEVector([1,2,3],[3,6,9])
x[1] = 10
@test collect(x) == [10,1,1,2,2,2,3,3,3]
x = RLEVector([1,2,3],[3,6,9])
x[9] = 10
@test collect(x) == [1,1,1,2,2,2,3,3,10]

# singleton ranges
x = RLEVector([1,2,2,3,3])
x[1] = 5
@test collect(x) == [5,2,2,3,3]
x = RLEVector([1,2,2,3,3])
x[1] = 2
@test collect(x) == [2,2,2,3,3]
x = RLEVector([1,1,2,3,3])
x[3] = 3
@test collect(x) == [1,1,3,3,3]
x = RLEVector([1,1,2,3,3])
x[3] = 1
@test collect(x) == [1,1,1,3,3]
x = RLEVector([3,1,2,1,3])
x[3] = 1
@test collect(x) == [3,1,1,1,3]

# ## range with scalar
# # in middle, no match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[10:11] .= 5
# @test x.runvalues == [1,2,3,5,3,4]
# @test x.runends == [4,8,9,11,12,16]
#
# # just left, no match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[9:10] .= 5
# @test x.runvalues == [1,2,5,3,4]
# @test x.runends == [4,8,10,12,16]
#
# # just right, no match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[11:12] .= 5
# @test x.runvalues == [1,2,3,5,4]
# @test x.runends == [4,8,10,12,16]
#
# # just left, match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[9:10] .= 2
# @test x.runvalues == [1,2,3,4]
# @test x.runends == [4,10,12,16]
#
# # just right, match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[11:12] .= 4
# @test x.runvalues == [1,2,3,4]
# @test x.runends == [4,8,10,16]
#
# # span, no match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[9:12] .= 5
# @test x.runvalues == [1,2,5,4]
# @test x.runends == [4,8,12,16]
#
# # span, left match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[9:12] .= 2
# @test x.runvalues == [1,2,4]
# @test x.runends == [4,12,16]
#
# # span, right match
# x = RLEVector([1,2,3,4],[4,8,12,16])
# x[9:12] .= 4
# @test x.runvalues == [1,2,4]
# @test x.runends == [4,8,16]
#
# # span, both match
# x = RLEVector([1,2,3,2],[4,8,12,16])
# x[9:12] .= 2
# @test x.runvalues == [1,2]
# @test x.runends == [4,16]
#
# # reverse range with scalar
# x = RLEVector([1,2,3,4],[2,4,6,8])
# x[4:-1:2] .= 5
# @test collect(x) == [1,5,5,5,3,3,4,4]
# @test x.runvalues == [1,5,3,4]
# @test x.runends == [1,4,6,8]
#
## range with vector
x = RLEVector([1,2,3,4],[2,4,6,8])
x[2:4] = [5,6,7]
@test collect(x) == [1,5,6,7,3,3,4,4]

## range with range
x = RLEVector([1,2,3,4],[2,4,6,8])
x[2:4] = 5:7
@test collect(x) == [1,5,6,7,3,3,4,4]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[1:2] = 5:6
@test collect(x) == [5,6,2,2,3,3,4,4]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[1:3] = 5:7
@test collect(x) == [5,6,7,2,3,3,4,4]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[7:8] = 5:6
@test collect(x) == [1,1,2,2,3,3,5,6]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[6:8] = 5:7
@test collect(x) == [1,1,2,2,3,5,6,7]

# range with vector and match
x = RLEVector([1,2,3,4],[2,4,6,8])
x[1:2] = [2,2]
@test collect(x) == [2,2,2,2,3,3,4,4]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[2:3] = [1,1]
@test collect(x) == [1,1,1,2,3,3,4,4]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[2:3] = [2,2]
@test collect(x) == [1,2,2,2,3,3,4,4]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[6:7] = [3,3]
@test collect(x) == [1,1,2,2,3,3,3,4]
x = RLEVector([1,2,3,4],[2,4,6,8])
x[7:8] = [3,3]
@test collect(x) == [1,1,2,2,3,3,3,3]

# reverse range with vector
x = RLEVector([1,2,3,4],[2,4,6,8])
x[4:-1:2] = [5,6,7]
@test x[4:-1:2] == [5,6,7]

# Colon
x = RLEVector([1,2,3,4],[2,4,6,8])
@test x[:] == x
#x[:] .= 4
#@test x == RLEVector([4 for i in 1:8])

# Logical
x = RLEVector([1,2],[2,4])
@test x[ [ true,true,true,false ] ] == RLEVector([1,2],[2,3])
#x[ [true,true,true,false] ] .= 4
#@test x == RLEVector([4,4,4,2])
x = RLEVector([1,2],[2,4])
x[ [true,true,true,false] ] = [4,5,6]
@test x == RLEVector([4,5,6,2])

# eachrange iterator
x = RLEVector([1, 1, 2, 2, 7, 12])
y = collect(eachrange(x))
@test y[1] == (1,1:2)
@test y[2] == (2,3:4)
@test y[3] == (7,5:5)
@test y[4] == (12,6:6)
@test length(eachrange(x)) == 4

# iterate
y = [1, 1, 2, 2, 7, 12]
x = RLEVector([1, 1, 2, 2, 7, 12])
out = Vector{eltype(x)}()
global next = iterate(x)
while next !== nothing
    (i,state) = next
    push!(out, i)
    global next = iterate(x, state)
end
@test out == y

# tapply
factor = repeat( ["a","b","c","d","e"], inner=5 )
rle = RLEVector( factor )
x = collect(1:25)
tapply_res = Dict( "a" => mean(x[1:5]), "b" => mean(x[6:10]), "c" => mean(x[11:15]), "d" => mean(x[16:20]), "e" => mean(x[21:25]) )
@test tapply_res == tapply( x, factor, mean )
@test tapply_res == tapply( x, rle, mean )

# tapply, non-unique RLE values
factor = repeat( ["a","b","c","d","e","b"], inner=5 )
rle = RLEVector( factor )
x = collect(1:30)
tapply_res = Dict( "a" => mean(x[1:5]), "b" => mean(x[vcat(6:10,26:30)]), "c" => mean(x[11:15]), "d" => mean(x[16:20]), "e" => mean(x[21:25]) )
@test_throws ArgumentError tapply( x, rle, mean )

end # testset

end # module


text:
using NBodyIPs, JuLIP, Test, IPFitting, DataFrames
using JuLIP.Potentials: evaluate_d
using IPFitting: Dat, LsqDB
using NBodyIPs: BondLengthDesc, BondAngleDesc, bodyorder
Lsq = IPFitting.Lsq
Err = IPFitting.Errors
using LinearAlgebra: norm

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

r0 = rnn(:Si)
calc = StillingerWeber()
data = generate_data(:Si, 2, 0.33*r0, 100, calc)

TRANSFORM = "exp( - 2 * (r/$r0 - 1.5) )"

rcut2 = cutoff(calc)*2.0
CUTOFF2 = PolyCutSym(2, rcut2)
D2 = BondLengthDesc(TRANSFORM, CUTOFF2)
rcut3 = cutoff(calc)*2.0
CUTOFF3 = PolyCut(2, rcut3)
D3 = BondLengthDesc(TRANSFORM, CUTOFF3)

B = [nbpolys(2, D2, 6); nbpolys(3, D3, 6)]

db = LsqDB("", B, data)

Ib2 = findall(b -> bodyorder(b) == 2, db.basis)
Ib3 = findall(b -> bodyorder(b) == 3, db.basis)

IP2, _ = Lsq.lsqfit( db,
                     E0 = 0.0,
                     configweights = Dict("rand" => 1.0),
                     obsweights   = Dict("E" => 30.0, "F" => 1.0),
                     combineIP = NBodyIP,
                     Ibasis = Ib2 )

# IP3a is the standard fit with Vref which is evaluated during Lsq assembly
IP3a, _ = Lsq.lsqfit( db,
                      Vref = IP2,
                      configweights = Dict("rand" => 1.0),
                      obsweights   = Dict("E" => 30.0, "F" => 1.0),
                      combineIP = NBodyIP,
                      Ibasis = Ib3 )

# precompute the observations for Vref in parallel
# (use JULIA_NUM_THREADS to use multi-threading for this!)
IP2f = fast(IP2)
Err.add_fits!(IP2f, db.configs, fitkey = "Vref")

# in the following comment Vref = IP2 is used only to construct the
# final potential but it will ot be evaluated, instead the precomputed
# observations from the `add_fits!` line are used.
IP3b, _ = Lsq.lsqfit( db,
                      Vref = IP2,
                      configweights = Dict("rand" => 1.0),
                      obsweights   = Dict("E" => 30.0, "F" => 1.0),
                      combineIP = NBodyIP,
                      Ibasis = Ib3 )

println(@test IP3a.components[1] == IP3b.components[1])
println(@test IP3a.components[2] === IP3b.components[2])
println(@test IP3a.components[3].c ≈ IP3b.components[3].c)
println(@test !(IP3a.components[3] === IP3b.components[3]))


text:
using PartialFunctions
using Test

a(x) = x^2
greet(greeting, name, punctuation) = "$(greeting), $(name)$(punctuation)"

@testset "Partial functions" begin
    @test map((+)$2, [1,2,3]) == [3, 4, 5]
    @test repr(map $ a) == "map(a, ...)"
    @test (map $ a)([1, 2, 3]) == [1, 4, 9]

    @test greet("Hello", "Bob", "!") == "Hello, Bob!"
    sayhello = greet $ "Hello"
    @test repr(sayhello) == "greet(\"Hello\", ...)"
    @test repr("text/plain", sayhello) == repr(sayhello)

    @test sayhello("Bob", "!") == "Hello, Bob!"
    hi_bob = greet $ "Hi" $ "Bob" $ "!"
    @test hi_bob isa PartialFunctions.PartialFunction{typeof(greet), Tuple{String, String, String}, NamedTuple{(), Tuple{}}}
    @test hi_bob <| () == "Hi, Bob!"
    @test sayhello <| ("Jimmy", "?")... == "Hello, Jimmy?"

    @test greet $ ("Hi", "Bob") <| "!" == "Hi, Bob!"
end

@testset "Reversed functions" begin
    revmap = flip(map)
    @test flip(revmap) == map
    @test revmap([1,2,3], sin) == map(sin, [1,2,3])

    func(x, y) = x - y
    func(x, y, z) = x - y - z
    @test func(1, 2) == -1
    @test func(1, 3, 6) == -8
    flipped = flip(func)
    @test flipped(2, 1) == -1
    @test flipped(3, 2, 1) == -4
end

@testset "Keyword Arguments" begin
    a = [[1,2,3], [1,2]]
    sort_by_length = sort $ (; by = length)
    @test sort(a, by = length) == sort_by_length(a)

    sort_a_by_length = sort $ (a, (;by = length))
    @test sort(a, by = length) == sort_a_by_length()

    sort_a_by_length_2 = sort $ ((a,), (;by = length))
    @test sort_a_by_length == sort_a_by_length_2

    @test repr(sort_a_by_length) == "sort([[1, 2, 3], [1, 2]], ...; by = length, ...)"
end


text:
clear()
@testset "check ScalePath op's gradient" begin
    # [1] prepare input data and its label
    input = randn(19,17,7);
    label = randn(19,17,7);
    x = Variable(input, type=Array{Float64});
    l = Variable(label, type=Array{Float64});

    # [2] prepare a learnabel scalar multiplier
    k = ScalePath(1.0, ndims=3, type=Array{Float64});

    # [3] forward and get gradient by backpropagation
    y  = forward(k, x);
    C₁ = mseCost(y, l);
    backward();
    ∇ = Delta.δ(k.scale)[1];  # backpropagation gradient

    # [4] with a samll change of weight, forward again
    Δ = 1e-6;
    Delta.ᵛ(k.scale)[1] += Δ;

    y  = forward(k, x);
    C₂ = mseCost(y, l);
    backward();

    # [5] check if the auto-grad is true or not
    ϵ    = eps(Float64)
    ∂C∂k = (C₂ - C₁)/Δ;            # numerical gradient
    err  = (∂C∂k - ∇)/(∇+ϵ);       # relative error
    err  = abs(err) * 100;         # relative error in %
    err  = err < 1e-1 ? 0.0 : err; # zero-error if under 0.1%
    @test err < 1.0
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using ImageDraw, Plots

function draw_boxes_and_landmarks(img, boxes, lms; conf=nothing, len=640)
    to_draw = deepcopy(RGB{N0f8}.(img[1,:,:], img[2,:,:], img[3,:,:]))
    colors = [
        RGB{N0f8}(1.0,0.0,0.0), # red
        RGB{N0f8}(0.0,0.0,1.0), # blue
        RGB{N0f8}(0.0,1.0,0.0), # green
        RGB{N0f8}(1.0,1.0,0.0), # yellow
        RGB{N0f8}(0.0,1.0,1.0), # turquoise
        RGB{N0f8}(1.0,0.0,1.0) # magenta
    ]

    for person in 1:size(boxes, 2)
        rect = Int.(ceil.(boxes[:, person]))
        if any(x->x >= 0, boxes[1:4, person])
            to_draw = draw!(to_draw, Polygon([
                        (rect[1, 1], rect[2, 1]),
                        (rect[3, 1], rect[2, 1]),
                        (rect[3, 1], rect[4, 1]),
                        (rect[1, 1], rect[4, 1])]), colors[1])
        end
        colorid = 2
        if lms !== nothing
            pts = Int.(ceil.(lms[:, person]))
            for lm in 1:2:10
                if pts[lm,1] >= 0 && pts[lm+1,1] >= 0
                   to_draw = draw!(to_draw, CirclePointRadius(pts[lm,1],pts[lm+1,1],2), colors[colorid])
                end
                colorid += 1
            end
        end
    end

    if conf !== nothing
        plt = plot(to_draw, size=(len,len))
        for person in 1:size(conf)[2]
            pt_x = boxes[1, person] + (boxes[3, person] - boxes[1, person]) / 2
            pt_y = boxes[2, person] - 10
            annotate!(pt_x, pt_y, text(string(round.(value(conf[2,person]); digits=4)), :orange, 6))
        end
        return plt
    else
        return to_draw
    end

end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
"""

Executable path corresponding to a source file, or a model.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
function executable_path(source_path::AbstractString)
    if Sys.iswindows()
        source_path =  source_path * ".exe"
    end
    source_path = source_path
end

executable_path(model::T) where T <: CmdStanModels =
    executable_path(model.tmpdir)

"""

Error thrown when a Stan model fails to compile.

$(TYPEDEF)

# Extended help

Accessing fields directly is part of the API.

$(FIELDS)
"""
struct StanModelError <: Exception
    name::String
    message::String
end

function Base.showerror(io::IO, e::StanModelError)
    print(io, "\nError when compiling SampleModel ", e.name, ":\n",
          e.message)
end

"""

Ensure that a compiled model executable exists, and return its path.

$(SIGNATURES)

# Extended help

If compilation fails, a `StanModelError` is returned instead.

Internal, not exported.
"""
function ensure_executable(model::T) where T <: CmdStanModels
    @unpack cmdstan_home, exec_path = model
    error_output = IOBuffer()
    is_ok = cd(cmdstan_home) do
        success(pipeline(`make -f $(cmdstan_home)/makefile -C $(cmdstan_home) $(exec_path)`;
                         stderr = error_output))
    end
    if is_ok
        exec_path
    else
        throw(StanModelError(model.name, String(take!(error_output))))
    end
end

"""

Compile a model, throwing an error if it failed.

$(SIGNATURES)
"""
function stan_compile(model::T) where T <: CmdStanModels
    ensure_executable(model)
    nothing
end


text:
Base.isapprox(x::Tuple, y::Tuple; kws...) = isapprox(collect(x), collect(y); kws...)

@testset "Analysis struct" begin
@testset "Uniformly distributed values" begin
    n = 20
    data = rand(Uniform(0,10), n)
    curve = CalculatedABC.ABCcurve(data)

    analysis1 = CalculatedABC.ABCanalysis(curve)
    @test isapprox(analysis1.pareto,      (0.4, 0.6192318028485929))
    @test isapprox(analysis1.break_even,  (0.5, 0.7210686469145638))
    @test isapprox(analysis1.demark_AB,   (0.4, 0.6192318028485929))
    @test isapprox(analysis1.submarginal, (0.6, 0.8125961014234772))

    analysis2 = CalculatedABC.ABCanalysis(data)
    @test isapprox(analysis2.pareto,      (0.4, 0.6192318028485929))
    @test isapprox(analysis2.break_even,  (0.5, 0.7210686469145638))
    @test isapprox(analysis2.demark_AB,   (0.4, 0.6192318028485929))
    @test isapprox(analysis2.submarginal, (0.6, 0.8125961014234772))
    @test analysis2.a_indices == [2, 13, 1, 3, 6, 14, 7, 19]
    @test analysis2.b_indices == [12, 16, 18, 4]
    @test analysis2.c_indices == [5, 15, 17, 9, 20, 11, 10, 8]
end

@testset "Exponentially distributed values" begin
    n = 20
    data = rand(Exponential(1), n)
    curve = CalculatedABC.ABCcurve(data)

    analysis1 = CalculatedABC.ABCanalysis(curve)
    @test isapprox(analysis1.pareto,      (0.35, 0.6503414525690906))
    @test isapprox(analysis1.break_even,  (0.35, 0.6503414525690906))
    @test isapprox(analysis1.demark_AB,   (0.35, 0.6503414525690906))
    @test isapprox(analysis1.submarginal, (0.5, 0.7704643555492489))

    analysis2 = CalculatedABC.ABCanalysis(data)
    @test isapprox(analysis2.pareto,      (0.35, 0.6503414525690906))
    @test isapprox(analysis2.break_even,  (0.35, 0.6503414525690906))
    @test isapprox(analysis2.demark_AB,   (0.35, 0.6503414525690906))
    @test isapprox(analysis2.submarginal, (0.5, 0.7704643555492489))
    @test analysis2.a_indices == [2, 19, 13, 20, 6, 4, 5]
    @test analysis2.b_indices == [15, 3, 12]
    @test analysis2.c_indices == [16, 7, 17, 9, 11, 14, 18, 8, 1, 10]
end

@testset "SwissInhabitants data" begin
    data = readdlm("swissinhabitants")[:,1]
    curve = CalculatedABC.ABCcurve(data)

    analysis1 = CalculatedABC.ABCanalysis(curve)
    @test isapprox(analysis1.pareto,      (0.25379834254143646, 0.7219586074342286))
    @test isapprox(analysis1.break_even,  (0.22168508287292818, 0.6913764432243525))
    @test isapprox(analysis1.demark_AB,   (0.22168508287292818, 0.6913764432243525))
    @test isapprox(analysis1.submarginal, (0.3715469613259669, 0.810278191608252))

    analysis2 = CalculatedABC.ABCanalysis(data)
    @test isapprox(analysis2.pareto,      (0.25379834254143646, 0.7219586074342286))
    @test isapprox(analysis2.break_even,  (0.22168508287292818, 0.6913764432243525))
    @test isapprox(analysis2.demark_AB,   (0.22168508287292818, 0.6913764432243525))
    @test isapprox(analysis2.submarginal, (0.3715469613259669, 0.810278191608252))
    # too many to test directly
    @test (analysis2.a_indices |> length) == 642
    @test (analysis2.b_indices |> length) == 434
    @test (analysis2.c_indices |> length) == 1820
    @test minimum(data[analysis2.a_indices]) >= maximum(data[analysis2.b_indices])
    @test minimum(data[analysis2.b_indices]) >= maximum(data[analysis2.c_indices])
end
end


text:
module TestPaths

using FilePathsBase

import Base: ==

__init__() = FilePathsBase.register(TestPath)

# We'll make a semicolon separate test path
struct TestPath <: AbstractPath
    parts::Tuple{Vararg{String}}
end

TestPath() = TestPath(tuple())

function TestPath(str::AbstractString)
    str = string(str)

    if isempty(str)
        return TestPath(tuple("."))
    end

    tokenized = split(str, ";")
    if isempty(tokenized[1])
        tokenized[1] = ";"
    end
    return TestPath(tuple(map(String, tokenized)...))
end

==(a::TestPath, b::TestPath) = a.parts == b.parts
Base.string(path::TestPath) = join([path.parts...], ";")
FilePathsBase.ispathtype(::Type{TestPath}, str::AbstractString) = startswith(str, "test|;;")
Base.show(io::IO, path::TestPath) = print(io, "p\"$path\"")

end


text:
using fci
using JLD2

@load "_testdata_h4_triplet_ymatrix.jld2"

norbs = 4      #for the h4 example
configs = [1,3,4], [4,2,3], [2], [1]
         #alpha, alpha, beta, beta
idxs = [3,4,2,1]

@testset "get index" begin
    for i in 1:4
        config = configs[i]
        idx = idxs[i]
        if i <= 2
            y = ya
            string = zeros(Int8, norbs)
            string[config] .= 1

            index = 1
            start = [1,1]

            for value in string
                if value == 0
                    #move down but dont add to index
                    start[1] = start[1]+1
                end
                if value == 1
                    #move right and add value to index
                    start[2] = start[2] + 1
                    index += y[start[1], start[2]]
                end
            end
            @test idx == index

        else
            y = yb
            string = zeros(Int8, norbs)
            string[config] .= 1

            index = 1
            start = [1,1]

            for value in string
                if value == 0
                    #move down but dont add to index
                    start[1] = start[1]+1
                end
                if value == 1
                    #move right and add value to index
                    start[2] = start[2] + 1
                    index += y[start[1], start[2]]
                end
            end
            @test idx == index
        end
    end
end

max_stars_repo_name: nbraunsc/fci

max_issues_repo_name: nbraunsc/fci

max_forks_repo_name: nbraunsc/fci


text:
import Base.LinAlg.BLAS: axpy!
import Base.broadcast!

abstract AbstractCuArray{T,N}
export AbstractCuArray

function Base.copy!{T,N}(y::AbstractCuArray{T,N}, x::AbstractCuArray{T,N})
    t = ctype(T)
    f = @nvrtc """
    $array_h
    __global__ void f(Array<$t,$N> y, Array<$t,$N> x) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        int len_y = y.length();
        if (idx < len_y) {
            if (x.length() == len_y) y(idx) = x(idx);
            else {
                int subs[$N];
                y.idx2sub(idx, subs);
                y(subs) = x(subs);
            }
        }
    } """
    f(y, x, dx=length(y))
    y
end


text:
using HyperGraphs
using LinearAlgebra: Diagonal

x1 = HyperGraph([HyperEdge([1, 2, 3]), HyperEdge([2, 3, 4])])
x2 = HyperGraph([HyperEdge(["1", "2", "3"]), HyperEdge(["2", "3", "4"])])
x3 = ChemicalHyperGraph([ChemicalHyperEdge(), ChemicalHyperEdge()])

# degree matrix
@test degree_matrix(x1) == degree_matrix(x2)

# incidence matrix
@test incidence_matrix(x1) == incidence_matrix(x2)

# cardinality matrix
@test cardinality_matrix(x1) == cardinality_matrix(x2) == Diagonal([3, 3])

# weight matrix
@test weight_matrix(x3) == Diagonal([1, 1])


text:
function elastancemodel!(system::CVSystem,n::Int64)
    # ventricular elastance
    CVModule.elastancefn!(system,n)

    # update chamber pressures
    system.heart.lv.P[n+1] = system.heart.lv.E[n+1]*(
        system.heart.lv.V[n+1] - system.heart.lv.V0);
    system.heart.rv.P[n+1] = system.heart.rv.E[n+1]*(
        system.heart.rv.V[n+1] - system.heart.rv.V0);

    system.heart.la.P[n+1] = system.heart.la.E*(
        system.heart.la.V[n+1] - system.heart.la.V0);
    system.heart.ra.P[n+1] = system.heart.ra.E*(
        system.heart.ra.V[n+1] - system.heart.ra.V0);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
function N=sphere2cart(lam,phi,r)
x=cos(lam)*cos(phi)*r;
y=sin(lam)*cos(phi)*r;
z=sin(phi)*r;
N=[x;y;z];
end


text:
# This file is a part of LegendGeSim.jl, licensed under the MIT License (MIT).

import Test
Test.@testset "Package LegendGeSim" begin

include("test_utils.jl")
# include("test_another_source_file.jl")
# ... more includes ...

end # testset


text:
const DIR = @__DIR__

@info "Installing Jupyter kernel..."
using IJulia
installkernel("Julia_SSPOC_demo", "--project=" * DIR * "/../")
@info "Build complete."


text:
###############################################################################
"""
Provides utilities to transform between physical (x,y,z) and pixel (u,v) coordinate systems.
"""
module ImageRectification
###############################################################################

###############################################################################
using StaticArrays              # provides SMatrix, SVector, for faster linear algebra
using ProgressMeter
###############################################################################

###############################################################################
"""
FindXYFromUVZ(u, v, z, DLT)

    Returns XY coordinates of a point.

    Calculated using the DLT formulation from eqn 3 in Holland et al. z is
    known a-priori (usually the water surface):
"""
function FindXYFromUVZ end
FindXYFromUVZ(u::Number, v::Number, z::Number, DLT::NamedTuple) = uvz2xy(u, v, z, DLT)
FindXYFromUVZ(uv::Array, z::Number, DLT::NamedTuple) = FindXYFromUVZ(view(uv,:,1), view(uv,:,2), z::Number, DLT)
FindXYFromUVZ(uv::Array, z::Array, DLT::NamedTuple) = FindXYFromUVZ(view(uv,:,1), view(uv,:,2), z::Array, DLT)
#
function FindXYFromUVZ(u, v, z::Number, DLT::NamedTuple)
    x = Array{Float64}(undef, size(u))
    y = similar(x)
    for i in eachindex(u); x[i],y[i] = uvz2xy(u[i], v[i], z, DLT); end
    return x,y
end
function FindXYFromUVZ(u, v, z::Array, DLT::NamedTuple)
    x = Array{Float64}(undef, size(u))
    y = similar(x)
    for i in eachindex(u); x[i],y[i] = uvz2xy(u[i], v[i], z[i], DLT); end
    return x,y
end
# Helper functions for FindXYFromUVZ():
Amatrix(u, v, z, DLT::NamedTuple) = SMatrix{2,2}( u*DLT.L₉-DLT.L₁,  v*DLT.L₉-DLT.L₅, u*DLT.L₁₀-DLT.L₂, v*DLT.L₁₀-DLT.L₆ )
bvector(u, v, z, DLT::NamedTuple) = SVector( DLT.L₄-u + (DLT.L₃-u*DLT.L₁₁)*z, DLT.L₈-v + (DLT.L₇-v*DLT.L₁₁)*z )
uvz2xy(u, v, z, DLT) = Amatrix(u, v, z, DLT) \ bvector(u, v, z, DLT)
###############################################################################

###############################################################################
""" Calculate the uv -> xyz coordinate conversion using each of the sets of DLTs: """
function CalculateModifiedxyCoordinates(u,v, WaterSurfaceElev, DLTs;
                                        wseNoiseFactor = 0.0)
    Perturb_wse = randn(size(DLTs,2)) .* wseNoiseFactor
    Modified_x, Modified_y = Rect.FindXYFromUVZ_DLT_z_sensitivity(u,v,
                                                                  # WaterSurfaceElev .+ randn(size(DLTs,2)) .* wseNoiseFactor,
                                                                  WaterSurfaceElev .+ Perturb_wse,
                                                                  DLTs)
    return Modified_x, Modified_y, Perturb_wse
end

###############################################################################

###############################################################################
# This option is for sensitivity analysis of water surface elevation (and a single DLT set):
function FindXYFromUVZ_WSE_sensitivity(u::Array, v::Array, z::Array, DLT::NamedTuple)
    x = Array{Float64}(undef, size(u)..., length(z))
    y = Array{Float64}(undef, size(v)..., length(z))
    @showprogress 1 "Calculating $(length(u)) xy coordinates using $(size(DLT,2)) DLT sets:" for d in 1:length(z)
        for j in 1:size(u,2), i in 1:size(u,1)
            x[i,j,d],y[i,j,d] = FindXYFromUVZ(u[i,j], v[i,j], z[d], DLT)
        end
    end
    return x,y
end
# This option is for sensitivity analysis of multiple DLT sets:
function FindXYFromUVZ_DLT_sensitivity(u::Array, v::Array, z::Float64, DLTs::Array)
    x = Array{Float64}(undef, size(u)..., size(DLT,2))
    y = Array{Float64}(undef, size(v)..., size(DLT,2))
    @showprogress 1 "Calculating $(length(u)) xy coordinates using $(size(DLT,2)) DLT sets:" for d in 1:size(DLT,2)
        for j in 1:size(u,2), i in 1:size(u,1)
            x[i,j,d],y[i,j,d] = FindXYFromUVZ(u[i,j], v[i,j], z, DLT[d])
        end
    end
    return x,y
end
# This option is for sensitivity analysis of multiple DLT sets, along with
# multiple z values, for analysis of sensitivity to, e.g., water surface
# elevation:
function FindXYFromUVZ_DLT_z_sensitivity(u::Array, v::Array, z::Array, DLT::Array)
    (length(z) != size(DLT,2)) ? @error("lengths of z and DLT do not match!") : nothing
    x = Array{Float64}(undef, size(u)..., size(DLT,2))
    y = Array{Float64}(undef, size(v)..., size(DLT,2))
    @showprogress 1 "Calculating $(length(u)) xy coordinates using $(size(DLT,2)) DLT sets:" for d in 1:size(DLT,2) for j in 1:size(u,2), i in 1:size(u,1)
            x[i,j,d],y[i,j,d] = FindXYFromUVZ(u[i,j], v[i,j], z[d], DLT[d])
        end
    end
    return x,y
end
###############################################################################

###############################################################################
###############################################################################
###############################################################################
end                             # end of module
###############################################################################
###############################################################################
###############################################################################

max_stars_repo_name: saschweitzer/camera-calibration

max_issues_repo_name: saschweitzer/camera-calibration

max_forks_repo_name: saschweitzer/camera-calibration


text:
# Illustrate the concept of mutual consistency
function illustrate_mc()
    dgp = illustration_dgp()
    assumptions = Dict{Symbol, Any}(:lb => 0, :ub => 1)
    basis1 = (interacted_bernstein_basis(1, notℓ = 2),
              interacted_bernstein_basis(1, notℓ = 2))
    basis2 = (interacted_bernstein_basis(3, notℓ = 1),
              interacted_bernstein_basis(3, notℓ = 1))
    bases = [basis1, basis2]
    tp = ey1(dgp, ℓ = 2)

    # for plotting
    ev = DataFrame(u = 0:.01:1, d = 1)
    results = DataFrame(u = ev.u)

    # Impose mutual consistency
    assumptions[:mutually_consistent] = true
    r = compute_bounds(tp, bases, assumptions, dgp)
    ev.z = fill([0,0], nrow(ev))
    results[:, "MTR (l=1 d=1 z2=0)"] = evaluate_mtr(r[:mtr_ub][1][2], ev)
    ev.z = fill([0,1], nrow(ev))
    results[:, "MTR (l=1 d=1 z2=1)"] = evaluate_mtr(r[:mtr_ub][1][2], ev)
    ev.z = fill([0,0], nrow(ev))
    results[:, "MTR (l=2 d=1 z1=0) con"] = evaluate_mtr(r[:mtr_ub][2][2], ev)
    ev.z = fill([1,0], nrow(ev))
    results[:, "MTR (l=2 d=1 z1=1) con"] = evaluate_mtr(r[:mtr_ub][2][2], ev)

    # Do not impose mutual consistency. Keep target parameter defined in terms
    # of the second model. Then the max MTR will be a second model MTR that is
    # inconsistent with the first model MTR.
    assumptions[:mutually_consistent] = false
    r = compute_bounds(tp, bases, assumptions, dgp)
    ev.z = fill([0,0], nrow(ev))
    results[:, "MTR (l=2 d=1 z1=0) incon"] = evaluate_mtr(r[:mtr_ub][2][2], ev)
    ev.z = fill([1,0], nrow(ev))
    results[:, "MTR (l=2 d=1 z1=1) incon"] = evaluate_mtr(r[:mtr_ub][2][2], ev)
    return results
end
export illustrate_mc

# Show the impact of mutual consistency for an instrument-invariant parameter
# like the ATT
function simulation_att()
    dgp = simulation_dgp()

    results = DataFrame(degree = 1:1:9)
    results[:, "LB (l = 1; MC off)"] .= NaN
    results[:, "UB (l = 1; MC off)"] .= NaN
    results[:, "LB (l = 2; MC off)"] .= NaN
    results[:, "UB (l = 2; MC off)"] .= NaN
    results[:, "LB (MC on)"] .= NaN
    results[:, "UB (MC on)"] .= NaN
    assumptions = Dict{Symbol, Any}(:lb => 0, :ub => 1)

    for k in 1:nrow(results)
        bases = [(interacted_bernstein_basis(results[k, :degree], notℓ = 2),
                  interacted_bernstein_basis(results[k, :degree], notℓ = 2)),
                 (interacted_bernstein_basis(results[k, :degree], notℓ = 1),
                  interacted_bernstein_basis(results[k, :degree], notℓ = 1))]

        assumptions[:mutually_consistent] = false
        tp = att(dgp, ℓ = 1)
        r = compute_bounds(tp, bases, assumptions, dgp)
        results[k, 2:3] = [r[:lb], r[:ub]]

        tp = att(dgp, ℓ = 2)
        r = compute_bounds(tp, bases, assumptions, dgp)
        results[k, 4:5] = [r[:lb], r[:ub]]

        assumptions[:mutually_consistent] = true
        r = compute_bounds(tp, bases, assumptions, dgp)
        results[k, 6:7] = [r[:lb], r[:ub]]
    end

    knots = vcat(0, dgp.pscore, 1)
    bases = [(interacted_constantspline_basis(knots, notℓ = 2),
              interacted_constantspline_basis(knots, notℓ = 2)),
             (interacted_constantspline_basis(knots, notℓ = 1),
              interacted_constantspline_basis(knots, notℓ = 1))]

    assumptions[:mutually_consistent] = false
    tp = att(dgp, ℓ = 1)
    r = compute_bounds(tp, bases, assumptions, dgp)
    results[:, "NP LB (l = 1; MC off)"] .= r[:lb]
    results[:, "NP UB (l = 1; MC off)"] .= r[:ub]

    tp = att(dgp, ℓ = 2)
    r = compute_bounds(tp, bases, assumptions, dgp)
    results[:, "NP LB (l = 2; MC off)"] .= r[:lb]
    results[:, "NP UB (l = 2; MC off)"] .= r[:ub]

    assumptions[:mutually_consistent] = true
    r = compute_bounds(tp, bases, assumptions, dgp)
    results[:, "NP LB (MC on)"] .= r[:lb]
    results[:, "NP UB (MC on)"] .= r[:ub]

    return results
end
export simulation_att

# Show that mutual consistency can also help tighten inference on instrument-
# dependent parameters.
function simulation_prte()
    dgp = simulation_dgp()
    δ = .2
    tp = prte_plusδpercent(dgp, δ, ℓ = 1)
    knots = vcat(0, dgp.pscore, (1 + δ) .* dgp.pscore, 1)

    results = DataFrame(degree = 1:1:9)
    results[:, "LB (MC off)"] .= NaN
    results[:, "UB (MC off)"] .= NaN
    results[:, "LB (l = 2; NP)"] .= NaN
    results[:, "UB (l = 2; NP)"] .= NaN
    results[:, "LB (l = 2; NP; Decr)"] .= NaN
    results[:, "UB (l = 2; NP; Decr)"] .= NaN
    results[:, "LB (l = 2; NP; Linear)"] .= NaN
    results[:, "UB (l = 2; NP; Linear)"] .= NaN

    for k in 1:nrow(results)
        bases = [(interacted_bernstein_basis(results[k, :degree], notℓ = 2),
                  interacted_bernstein_basis(results[k, :degree], notℓ = 2)),
                 (interacted_constantspline_basis(knots, notℓ = 1),
                  interacted_constantspline_basis(knots, notℓ = 1))]

        assumptions = Dict{Symbol, Any}(:lb => 0, :ub => 1,
                                        :mutually_consistent => false)
        r = compute_bounds(tp, bases, assumptions, dgp)
        results[k, 2:3] = [r[:lb], r[:ub]]

        assumptions[:mutually_consistent] = true
        r = compute_bounds(tp, bases, assumptions, dgp)
        results[k, 4:5] = [r[:lb], r[:ub]]

        assumptions[:decreasing_level] = [(2, 0), (2,1)]
        assumptions[:decreasing_difference] = [2]
        r = compute_bounds(tp, bases, assumptions, dgp)
        results[k, 6:7] = [r[:lb], r[:ub]]

        bases[2] = (interacted_bernstein_basis(1, notℓ = 1),
                    interacted_bernstein_basis(1, notℓ = 1))
        r = compute_bounds(tp, bases, assumptions, dgp)
        results[k, 8:9] = [r[:lb], r[:ub]]
    end

    return results
end
export simulation_prte

# Illustrate misspecification caused for a PRTE by using the wrong choice model
function prte_misspecification(; ℓ_gen = 1)
    dgp = prte_dgp(ℓ_gen = ℓ_gen) # ℓ_gen shouldn't matter here (and doesn't)
    knots = vcat(0, dgp.pscore, 1)

    # Set up PRTE target parameter
    dgp_new = DGP(suppZ = dgp.suppZ,
                  pscore = dgp.pscore,
                  densZ = fill(.25, 4),
                  mtrs = dgp.mtrs)
    tp = prte_newz(dgp, dgp_new, ℓ = 1)

    nrows = 5
    results = DataFrame(name = fill("", nrows),
                        lb = fill(+Inf, nrows), ub = fill(+Inf, nrows))

    # True value
    truth = eval_tp(tp, [dgp.mtrs], dgp)
    results[1, :name] = "True value"
    results[1, :lb], results[1, :ub] = truth, truth

    # IAM model using both instruments
    # --> should get bounds that are non-positive
    bases = [(constantspline_basis(knots), constantspline_basis(knots))]
    assumptions = Dict{Symbol, Any}(:lb => 0, :ub => 1)
    r = compute_bounds(tp, bases, assumptions, dgp)
    results[2, :name] = "Both instruments with IAM"
    results[2, :lb] = r[:lb]
    results[2, :ub] = r[:ub]

    # Use the first instrument only
    bases = [(interacted_constantspline_basis(knots, notℓ = 2),
              interacted_constantspline_basis(knots, notℓ = 2))]
    r = compute_bounds(tp, bases, assumptions, dgp)
    results[3, :name] = "Instrument 1 with IAM"
    results[3, :lb] = r[:lb]
    results[3, :ub] = r[:ub]

    # Use the second instrument only
    bases = [(interacted_constantspline_basis(knots, notℓ = 1),
              interacted_constantspline_basis(knots, notℓ = 1))]
    r = compute_bounds(tp, bases, assumptions, dgp)
    results[4, :name] = "Instrument 2 with IAM"
    results[4, :lb] = r[:lb]
    results[4, :ub] = r[:ub]

    # Use both instruments together and impose MC
    assumptions = Dict{Symbol, Any}(:lb => 0, :ub => 1,
                                    :mutually_consistent => true)
    bases = [(interacted_constantspline_basis(knots, notℓ = 2),
              interacted_constantspline_basis(knots, notℓ = 2)),
             (interacted_constantspline_basis(knots, notℓ = 1),
              interacted_constantspline_basis(knots, notℓ = 1))]
    r = compute_bounds(tp, bases, assumptions, dgp)
    results[5, :name] = "Both instruments with PM"
    results[5, :lb] = r[:lb]
    results[5, :ub] = r[:ub]
    return results
end
export prte_misspecification


text:
import HypothesisTests.LjungBoxTest
import HypothesisTests.BoxPierceTest

#
#
# """
#     LjungBoxTest(d::UncertainDataset; n::Int = 10000,
#         lag::Int = 1, dof::Int = 0) -> LjungBoxTest
#
# Compute the Ljung-Box `Q` statistic to test the null hypothesis of
# independence in a data series represented by the uncertain dataset `d`.
#
# `lag` specifies the number of lags used in the construction of `Q`. When
# testing the residuals of an estimated model, `dof` has to be set to the number
# of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
# model, set `dof=p+q`.
# """
# function LjungBoxTest(d::UncertainDataset; lag::Int = 1, dof::Int = 0)
#     LjungBoxTest(resample(d, 1)[1], lag, dof)
# end
#
# """
#     LjungBoxTest(d::UncertainDataset, n_tests::Int = 100; lag::Int = 1, dof::Int = 0) -> LjungBoxTest
#
# Compute the Ljung-Box `Q` statistic to test the null hypothesis of
# independence in a data series represented by the uncertain dataset `d`
# by performing the test on `n` independent draws of the dataset.
#
# `lag` specifies the number of lags used in the construction of `Q`. When
# testing the residuals of an estimated model, `dof` has to be set to the number
# of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
# model, set `dof=p+q`.
# """
# function LjungBoxTest(d::UncertainDataset, n_tests::Int = 100; lag::Int = 1, dof::Int = 0)
#     [LjungBoxTest(resample(d)[1], lag, dof) for i = 1:n_tests]
# end
#
# """
#     BoxPierceTest(d::UncertainDataset;
#         lag::Int = 1, dof::Int = 0) -> Vector{BoxPierceTest}
#
# Compute the Box-Pierce `Q` statistic to test the null hypothesis of
# independence in a data series represented by the uncertain dataset `d`.
#
# `lag` specifies the number of lags used in the construction of `Q`. When
# testing the residuals of an estimated model, `dof` has to be set to the number
# of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
# model, set `dof=p+q`.
# """
# function BoxPierceTest(d::UncertainDataset; lag::Int = 1, dof::Int = 0)
#     BoxPierceTest(resample(d, 1)[1], lag, dof)
# end
#
#
# """
#     BoxPierceTest(d::UncertainDataset, n_tests::Int = 100;
#         lag::Int = 1, dof::Int = 0) -> BoxPierceTest
#
# Compute the BoxPierceTest `Q` statistic to test the null hypothesis of
# independence in a data series represented by the uncertain dataset `d`
# by performing the test on `n` independent draws of the dataset.
#
# `lag` specifies the number of lags used in the construction of `Q`. When
# testing the residuals of an estimated model, `dof` has to be set to the number
# of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
# model, set `dof=p+q`.
# """
# function BoxPierceTest(d::UncertainDataset, n_tests::Int = 100;
#         lag::Int = 1, dof::Int = 0)
#     [BoxPierceTest(resample(d, 1)[1], lag, dof) for i = 1:n_tests]
# end
#
#
# export
# LjungBoxTest,
# BoxPierceTest


text:
using MiniGNN
using JLD, HDF5
# Beware! https://github.com/JuliaIO/JLD.jl/issues/234
using Flux.Data: DataLoader
using CUDA
using MiniFastai
using Flux: logitcrossentropy, onehotbatch, Descent, Dense, Dropout, relu, Chain
using SparseArrays
using MiniGNN: create_batch

@load "train_data.jld" train_data
#size Nsamples x V
@load "train_labels.jld" train_labels
@load "test_data.jld" test_data
@load "test_labels.jld" test_labels

#small = 1000
#train_data = train_data[1:small,:]
#train_labels = train_labels[1:small]
train_labels = onehotbatch(train_labels, 0:9)
#size Nlabels x Nsamples
test_labels = onehotbatch(test_labels,0:9)
d = load("List_graphs.jld")
L_g = d["L_g"]
sizes = d["sizes"]

bs = 100

train_loader = DataLoader((permutedims(train_data),train_labels), batchsize=bs, shuffle=true)
val_loader = DataLoader((permutedims(test_data),test_labels), batchsize=bs, shuffle=false)
databunch = Databunch(train_loader,val_loader)

batch, labels = create_batch(100,train_loader)
input = convert(Array{Float32,3},process_data(batch))

p=0.5
mat1 = MiniGNN.make_rescaled_mat(L_g[1])
layer1 = MiniGNN.GraphConv(mat1,1=>32,5)
mp = MiniGNN.GraphMaxPool(4)
mat2 = MiniGNN.make_rescaled_mat(L_g[3])
layer2 = MiniGNN.GraphConv(mat2,32=>64,5)
dim_fc = Int(64*size(mat2)[1]/4)
model = Chain(layer1 , mp, layer2, mp, MiniGNN.graphflatten, Dense(dim_fc,512, relu), Dropout(p), Dense(512,10))

learning_rate = 0.05
SGD = Descent(learning_rate)
loss(pred,y) = logitcrossentropy(pred,y)

learner_gnn = Learner(model,SGD,loss,databunch,cbs=(AvgStatsCallback(),))
transform = x -> convert(Array{Float32,3},process_data(x))
fit!(learner_gnn,30,transform)

model(input)
#one batch
using Flux
using Zygote: pullback, update!
ps = params(model)
            #source https://fluxml.ai/Flux.jl/stable/training/training/
            # https://fluxml.ai/Zygote.jl/latest/adjoints/#Pullbacks-1
y_pred, back_net = pullback(() -> model(input), ps)
current_loss, back_loss = pullback(y -> loss(y,labels),y_pred)
            #current_loss, back = pullback(() -> los_f(mod(xb),yb), ps)
gs = back_net(back_loss(1)[1])#back(one(current_loss))
update!(SGD,ps,gs)


text:
export TiledIndex, SampleSet, ExplicitSampleSet, TiledSampleSet
export zero_lattice_component, shift_lattice_component, indextype

# Sample Set Indices
struct TiledIndex{D}
    on_lattice::SVector{D,Int}
    off_lattice::Int
end
TiledIndex{D}(i::Int) where {D} = TiledIndex(zero(SVector{D,Int}), i)
Base.zero(::Type{TiledIndex{D}}) where {D} = TiledIndex{D}(0)
function Base.hash(x::TiledIndex, h::UInt=zero(UInt))
    hash(reduce((r, x) -> (1 << 19 - 1)*r + x, x.on_lattice .% UInt), x.off_lattice + h)
    # hash((x.on_lattice.data..., x.off_lattice), h)
    # reduce((r, x) -> 31*r + x, SVector(h, (x.on_lattice.data .% UInt)..., x.off_lattice % UInt))
end
zero_lattice_component(i::TiledIndex{D}) where {D} = TiledIndex(zero(SVector{D,Int}), i.off_lattice)
shift_lattice_component(i::TiledIndex, shift) = TiledIndex(i.on_lattice + shift, i.off_lattice)
const SampleIndex = Union{Int,TiledIndex}

# Sample Sets
abstract type SampleSet end
## ExplicitSampleSet
mutable struct ExplicitSampleSet{S<:State} <: SampleSet
    init::S
    goal_samples::Vector{S}
    V::Vector{S}
end
Base.eltype(nodes::ExplicitSampleSet{S}) where {S} = S
Base.length(nodes::ExplicitSampleSet) = length(nodes.V)
Base.getindex(nodes::ExplicitSampleSet, i::Int) = i > 0 ? nodes.V[i] : (i == 0 ? nodes.init : nodes.goal_samples[-i])
Base.eachindex(nodes::ExplicitSampleSet) = eachindex(nodes.V)
indextype(nodes::ExplicitSampleSet) = Int
# node_info_datastructure(nodes::ExplicitSampleSet, f) = [f(nodes[i]) for i in eachindex(nodes)]
@inline function node_info_datastructure(nodes::ExplicitSampleSet, ::Type{T}, f) where {T}
    DefaultDict{indextype(nodes),T}(i -> f(nodes[i]), passkey=true)
end

## TiledSampleSet
mutable struct TiledSampleSet{S<:State,M} <: SampleSet
    init::S
    goal_samples::Vector{S}
    A::M
    b::S
    Voff::Vector{S}
end
Base.eltype(nodes::TiledSampleSet{S}) where {S} = S
Base.length(nodes::TiledSampleSet) = length(nodes.Voff)
function Base.getindex(nodes::TiledSampleSet, i::TiledIndex)
    if i.off_lattice > 0
        nodes.A*i.on_lattice + nodes.b + nodes.Voff[i.off_lattice]
    elseif i.off_lattice == 0
        nodes.init
    else
        nodes.goal_samples[-i.off_lattice]
    end
end
function Base.getindex(nodes::TiledSampleSet, i::Int)
    @assert i <= 0 "::Int TiledSampleSet indices must be <= 0 (corresponding to init or goal_samples)"
    i == 0 ? nodes.init : nodes.goal_samples[-i]
end
indextype(nodes::TiledSampleSet) = TiledIndex{size(nodes.A, 2)}
@inline function node_info_datastructure(nodes::TiledSampleSet, ::Type{T}, f) where {T}
    DefaultDict{indextype(nodes),T}(i -> f(nodes[i]), passkey=true)
end

max_stars_repo_name: pbpf/clean

max_issues_repo_name: pbpf/clean

max_forks_repo_name: pbpf/clean


text:
function LutherFieldIntegrator(arg0::Field, arg1::RealFieldElement)
    return LutherFieldIntegrator((Field, RealFieldElement), arg0, arg1)
end

function get_a(obj::LutherFieldIntegrator)
    return jcall(obj, "getA", Vector{Vector{RealFieldElement}}, ())
end

function get_b(obj::LutherFieldIntegrator)
    return jcall(obj, "getB", Vector{RealFieldElement}, ())
end

function get_c(obj::LutherFieldIntegrator)
    return jcall(obj, "getC", Vector{RealFieldElement}, ())
end


text:
using Random
using JSON

mutable struct Network
    layers::Vector{Int64}
    neurons::Vector{Vector{Float64}}
    weights::Vector{Matrix{Float64}}
    biases::Vector{Vector{Float64}}
    z::Vector{Vector{Float64}}
    batch::Int64
    epochs::Int64
    η::Float64
    act
    act_p
    C
    ∇C_a

    function Network(layers,batch,epochs,η,act,act_p,C,∇C_a)
        new(layers::Vector{Int64},
        [[0.0 for j=1:layers[i]] for i=1:length(layers)],
        [randn(Float64,layers[i+1],layers[i]) for i=1:(length(layers)-1)],
        [rand!(MersenneTwister(1234), zeros(layers[i])) for i=1:(length(layers)-1)],
        [[0.0 for j=1:layers[i]] for i=1:length(layers)],
        batch::Int64,epochs::Int64,η::Float64,act,act_p,C,∇C_a)
    end
end

function feedforward!(net::Network,input::Vector{Float64})::Vector{Float64}
    net.layers[1] = input

    # a^(l+1) = σ(w^l*a^l + b^1)
    for i=1:(length(net.layers)-1)
        net.z[i] = net.weights[i]*net.neurons[i] + net.biases[i]
        net.neurons[i+1] = net.z[i].net.act()
    end

    return net.neurons[end]
end

function backprop(net::Network,input::Vector{Float64},label::Vector{Float64})::Vector{Vector{Float64}}
    feedforward!(net,input)

    # δ^L = ∇C_a ⊙ σ'(z^L)
    output_error = net.∇C_a(net.neurons[end],label).*net.act_p(net.z[end])
    errors = [output_error]

    for i=1:(length(net.layers)-1)
        # δ^l=((w^l+1)T δ^l+1) ⊙ σ′(z^l)
        next_error = (net.weights[length(net.layers)-i+1]*errors[1]).*net.act_p(net.z[length(net.layers)-i])
        prepend!(errors,next_error)
    end

    return errors, net.neurons
end

function sgd!(net::Network,train::Vector{Tuple{Vector{Float64},Vector{Float64}}})::Nothing
    for e=1:net.epochs
        train = shuffle(train)
        batches = [train[(1+(i*net.batch)):((i+1)*net.batch)] for i=0:(Int(length(train)/net.batch)-1)]

        for b in batches
            errors = []
            activations = []

            for x in b
                push!(errors,backprop(net,x[1],x[2])[1])
                push!(activations,backprop(net,x[1],x[2])[2])
            end

            # change w & b using errors = []
            for w1=1:length(net.weights)
                net.weights[wl] = net.weights[wl].-(net.η/net.batch).*sum([errors[e][wl]*transpose(activations[e][wl-1]) for e=1:length(errors)])

            end
            for bl=length(net.biases)
                net.biases[bl] = net.biases[bl].-(net.η/net.batch).*sum([errors[e][bl] for e=1:length(errors)])

            end
        end
    end

    return nothing
end

function evaluate(net::Network,test::Vector{Tuple{Vector{Float64},Vector{Float64}}},rule::Function)::Float64
    c = sum([1 for t in test if rule(feedforward(net,t[0])) == rule(t[1])])
    return "$(c) / $(length(test)) = $(100*(c/length(test)))%"
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
taxid = ARGS[1]
outfile = ARGS[2]

regex = "^$taxid\t(.*?)\t.*?\t.*?\t.*?\t(.*?)\t(.*?)\t.*"

open(outfile, "w") do fp
	open("data/gene/gene2accession", "r") do fp2
		for line in eachline(fp2)
		  	m = match(Regex(regex), line)
			if !isnothing(m)
				if m[1] != "-" && m[3] != "-"
					write(fp, "$(m[1])\t$(m[2])\n")
				end
			end
		end
	end
end

max_stars_repo_name: manabuishii/mesh-workflow

max_issues_repo_name: manabuishii/mesh-workflow

max_forks_repo_name: manabuishii/mesh-workflow


text:
"""
$(SIGNATURES)

    Set the local cache.
"""
function setLocalCache(cache::Union{LocalCache, Nothing})::Nothing
    global __localCache
    __localCache = cache
    return nothing
end
"""
$(SIGNATURES)

    Get the keyed local cache.
"""
function getLocalCache()::Union{Nothing, LocalCache}
    global __localCache
    return __localCache
end

"""
$SIGNATURES

    Set to true if GraffSDK should only use local cache.
"""
function forceOnlyLocalCache(isToForce::Bool)::Nothing
    global __forceOnlyLocalCache
    __forceOnlyLocalCache = isToForce
    return nothing
end

"""
$SIGNATURES

    Internal method to get from cache.
"""

function getElement(id::String)::Union{Any, Nothing}
    try
        global __localCache
        if __localCache == nothing
            @debug "Local cache is not set, skipping!"
            return nothing
        end
        numNodes = length(__localCache.collection, Mongoc.BSON("""{ "id": "$(id)" }""") )
        numNodes == 0 && return nothing
        result = Mongoc.find_one(__localCache.collection, Mongoc.BSON("""{ "id": "$(id)" }""") )
        if(typeof(result) == Mongoc.BSON)
            dict = Mongoc.as_dict(result)
            return JSON2.read(dict["data"], BigDataElementResponse)
        else
            @show result
            return JSON2.read(result["data"], BigDataElementResponse)
        end
    catch ex
        # Never throw exception if cache fails - just provide error and carry on.
        @error "[GraffSDK] Get cache failed - $ex"
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        @error "Error! $err"
        return nothing
    end
end

"""
$SIGNATURES

    Internal method to set in cache.
"""
function setElement(id::String, element::BigDataElementResponse)::Nothing
    try
        global __localCache
        if __localCache == nothing
            @debug "Local cache is not set, skipping!"
            return nothing
        end
        selector = Mongoc.BSON("""{ "id": "$(id)" }""")
        numNodes = length(__localCache.collection, selector )
        bsonDoc = Mongoc.BSON()
        bsonDoc["id"] = id
        bsonDoc["data"] = JSON2.write(element)
        if numNodes > 0 # exists
            # Set wrapper...
            setDoc = Mongoc.BSON()
            setDoc["\$set"] = bsonDoc
            Mongoc.update_one(__localCache.collection, selector, setDoc)
        else
            push!(__localCache.collection, bsonDoc)
        end
        return nothing
    catch ex
        # Never throw exception if cache fails - just provide error and carry on.
        @error "[GraffSDK] Set cache failed - $ex"
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        @error "Error! $err"
        return nothing
    end
end

"""
$SIGNATURES

    Internal method to delete from cache.
"""
function deleteElement(id::String)::Nothing
    try
        global __localCache
        if __localCache == nothing
            @debug "Local cache is not set, skipping!"
            return nothing
        end
        selector = Mongoc.BSON("""{ "id": "$(id)" }""")
        numNodes = length(__localCache.collection, selector )
        if numNodes > 0 # exists
            Mongoc.delete_one(__localCache.collection, selector)
        end
        return nothing
    catch ex
        # Never throw exception if cache fails - just provide error and carry on.
        @error "[GraffSDK] Delete cache failed - $ex"
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        @error "Error! $err"
        return nothing
    end
end


text:
using Distributed
using SharedArrays,Plots
addprocs(48)
nprocs()
@everywhere using FRGRealTime,QuadGK,Dierckx,DelimitedFiles

# u0 = [-8.0, 188498.0]
-20.0, 136 * 188498.0
v1 = SharedArray(collect(25:1.0:400.0))
v2 = SharedArray(collect(145:1.0:400.0))

@everywhere function mreturnfun(T)
    sol_temp = FRGRealTime.ZeroSolve(
        T,
        4.0,
        1.0,
        800.0,
        u0=[-8.0, 116585.3],
        # dtmax = 5.0,
        atolint=1e-4,
        rtolint=1e-4,
    )
    println("T=", T)
    return sqrt(-sol_temp[end][2])
end

@everywhere function mreturnfun2(T)
    sol_temp = FRGRealTime.ZeroLPASolve(
        T,
        4.0,
        1.0,
        800.0,
        u0=[-8.0, 116585.3],
    )
    println("T=", T)
    return sqrt(-sol_temp[end][2])
end

outv1 = pmap(mreturnfun, v1)
outv2 = pmap(mreturnfun2, v2)

plot!(v1,outv1)

dir = "/home/wjfu1/tyy/FRG-RealTime-data/Zero/ThermalMass"
writedlm("$dir/realtime_T.dat",v1)
writedlm("$dir/LPA_T.dat",v2)
writedlm("$dir/realtime_m2.dat",outv1)
writedlm("$dir/LPA_m2.dat",outv2)

116585.3
13108
20.43287
temp_T=54.0
sol = FRGRealTime.ZeroSolve(
    temp_T,
    4.0,
    1.0,
    800.0,
    u0=[-8.0, 116585.3],
    # dtmax = 5.0,
    # dtmax=0.2,
    atolint=1e-6,
    rtolint=1e-6,
)

plot(sol.t,-3*sol[1,:])
plot(sol.t,-sol[2,:])

<<<<<<< HEAD
plot(sol.t,-sol[1,:])

plot(sol.t,-sol[2,:])
=======
-sol[2,:][end]

-sol[2,:][end]
>>>>>>> c1b5f5caeab5b4fe6926368d82a33f6b833eefa9

plot(k->sqrt(-sol(k)[2]+k^2),1.0,2.0)

sqrt(-sol[2,:][end])

sol = FRGRealTime.ZeroLPASolve(
        145.0,
        4.0,
        1.0,
        800.0,
        u0=[-8.0, 116585.3],
)

<<<<<<< HEAD
sol=FRGRealTime.ZeroLPASolve(
    145.0,
    4.0,
    1.0,
    800.0,
    u0 = [-8.0, 116585.3],
)

sol=FRGRealTime.ZeroSolve(
    100.0,
    4.0,
    1.0,
    800.0,
    u0 = [-8.0, 116585.3],
    # dtmax = 5.0,
    atolint = 1e-4,
    rtolint=1e-4,
)

sqrt(-sol[end][2])
=======
sol = FRGRealTime.ZeroLPASolve(
        145.0,
        4.0,
        1.0,
        800.0,
        u0=[-14.0, 116585.3],
)

>>>>>>> c1b5f5caeab5b4fe6926368d82a33f6b833eefa9

-sol[2,:][end]

1
sol2 =
    FRGRealTime.ZeroSolve(145.0, 4.0, 1.0, 800.0, u0=sol(800.0), dtmax=0.5)

plot(sol.t,-sol[1,:])

<<<<<<< HEAD
-quadgk(
    k -> FRGRealTime.propReZeroflow(k, 60.0, 4.0, sol, 800.0, atol = 1e-6, rtol = 1e-6),
    1.0,
    800.0,
    atol = 1e-8,
    rtol = 1e-8,
)[1]+188498.0
xaxis
=======
plot!(v1,outv1)
plot!(v1,outv2)
>>>>>>> c1b5f5caeab5b4fe6926368d82a33f6b833eefa9

plot(sol.t,sol[1,:])

<<<<<<< HEAD
writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/Zero/realtime_zero_T=100_m2k.dat",-hcat(sol.u...)[2,:])
writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/Zero/realtime_zero_T=100_lamdak.dat",2*hcat(sol.u...)[1,:])
writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/Zero/realtime_zero_T=100_k.dat",sol.t)
=======
dir3
dir3 = "/home/wjfu1/tyy/FRG-RealTime-data/SpectralFunction/ps=0/LPAT=145"
writedlm("$dir3/m2k.dat",-hcat(sol.u...)[2,:])
writedlm("$dir3/lamdak.dat",2 * hcat(sol.u...)[1,:])
writedlm("$dir3/k.dat",sol.t)
writedlm("$dir3/Imp0_small.dat",v1)
writedlm("$dir3/ImGamma2_small.dat",outvIm4)
writedlm("$dir3/Rep0.dat",p0v1)
writedlm("$dir3/Reps.dat",psv1)
writedlm("$dir3/ReGamma2.dat",outvRe2)

>>>>>>> c1b5f5caeab5b4fe6926368d82a33f6b833eefa9

writedlm("$dir3/Imp0.dat",v1)
writedlm("$dir3/ImI2p0-q0.dat",outvIm4)

plot!(sol2.t,-sol2[2,:])
plot!(sol2.t,sol2[1,:])

sol = FRGRealTime.ZeroLPASolve(
    100.0,
    4.0,
    1.0,
    800.0,
    u0=[-8.0, 0.8423 * 116585.3],
)

sqrt(-sol[end][2])

m2fun(x) = -sol2(x)[2]
lamdafun(x) = sol2(x)[1]

k_temp = collect(1:1.0:800)
qs_temp = collect(collect(1:1.0:800)')

Voutvec = FRGRealTime.lamdabar_Zero.(
    Epi.(k_temp, m2fun.(k_temp)),
    qs_temp,
    k_temp,
    145.0,
    4.0,
    m2fun,
    lamdafun,
    800.0;
    atol=1e-6,
    rtol=1e-6,
)

k_temp = 800.0
qs_temp = 0.00001

using FRGRealTime

FRGRealTime.lamdabar_Zero(
    Epi.(k_temp, m2fun.(k_temp)),
    qs_temp,
    k_temp,
    145.0,
    4.0,
    m2fun,
    lamdafun,
    800.0;
    atol=1e-6,
    rtol=1e-6,
)

writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/realtime_zero_vertex.dat",Voutvec)

@eval @everywhere sol2 = $sol2
@everywhere m2fun(x) = -sol2(x)[2]
@everywhere lamdafun(x) = sol2(x)[1]

@everywhere k_temp = collect(1:20.0:800)
qs_temp = SharedArray(collect(collect(1:20.0:800)'))

Voutvec = map((k, qs) -> FRGRealTime.lamdabar_Zero(
    Epi.(k, m2fun.(k)),
    qs,
    k,
    145.0,
    4.0,
    m2fun,
    lamdafun,
    800.0;
    atol=1e-6,
    rtol=1e-6,
),k_temp,qs_temp)

plot(k_temp,outv1)

k -> lam4piflow(k, m2fun(k), 145.0, 4.0, lamdafun(k))
plot(sol2.t[end - 800:end],sol2[1,:][end - 800:end])
plot(sol2.t,sol2[1,:])


text:
"Bill's Gates list of recommended books" => [Dict("author" => b.author,
                                                  "title" => b.title,
                                                  "cover" => b.cover) for b in @vars(:books)]

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
#: by David Thierry
using JuMP
using Clp
tHorz = 2
m = Model()

@variable(m, 0 <= vCapCombTurb)
@variable(m, 0 <= vPowCombTurb[1:tHorz])

@variable(m, 0 <= ePowSteamTurb[1:tHorz])
@variable(m, 0 <= vCapSteamTurb)

@variable(m, 0 <= ePowOut[1:tHorz])
@variable(m, 0 <= ePowHp[1:tHorz])
@variable(m, 0 <= ePowIp[1:tHorz])
@variable(m, 0 <= vCapTransInter)

@variable(m, 0 <= eFuel[1:tHorz])

@variable(m, 0 <= eSteamIpLp[1:tHorz])
@variable(m, 0 <= eSteamLp1[1:tHorz])

@variable(m, 0 <= ePowLp1[1:tHorz])
@variable(m, 0 <= ePowLp2[1:tHorz])

@variable(m, 0 <= eCo2Fuel[1:tHorz])

# Pcc
@variable(m, 0 <= eCo2CapPcc[1:tHorz])

@variable(m, 0 <= eSteamUsePcc[1:tHorz])
@variable(m, 0 <= ePowUsePcc[1:tHorz])

@variable(m, 0 <= eCo2PccOut[1:tHorz])

@variable(m, 0 <= vCo2PccVent[1:tHorz])
@variable(m, 0 <= vCapPcc)

# Dac
@variable(m, 0 <= eCo2DacFlueIn[1:tHorz])

@variable(m, 0 <= eSteamUseDac[1:tHorz])
@variable(m, 0 <= ePowUseDac[1:tHorz])

@variable(m, 0 <= vCapDac)

@variable(m, 0 <= a0[0:tHorz])  # We could generalize this for any number of hours.
@variable(m, 0 <= a1[0:tHorz])
#@variable(m, 0 <= a2[1:tHorz])

@variable(m, 0 <= aR0[0:tHorz])
#@variable(m, 0 <= aR1[1:tHorz])

@variable(m, 0 <= sF[0:tHorz])
@variable(m, 0 <= sS[0:tHorz])

@variable(m, 0 <= vAbs[0:tHorz])
@variable(m, 0 <= vReg[0:tHorz])

@variable(m, 0 <= eCo2StoreDac[1:tHorz])
@variable(m, 0 <= eCo2CapDac[1:tHorz])
@variable(m, 0 <= eCo2DacVent[1:tHorz])

@variable(m, 0 <= eCo2Comp[1:tHorz])
@variable(m, 0 <= ePowUseComp[1:tHorz])
@variable(m, 0 <= vCapComp)
@variable(m, 0 <= eCo2Vent[1:tHorz])
# Parameters
pHeatRateCombTur = 1.
# pSteamRate = 1.
pEmissFactor = 1.

# pCo2CapRatePcc = 1.

#: Turbine Parameters
pFuelPowHp = 1.
pFuelPowIp = 1.
pFuelSteamIpLp = 1.
pHeatRateLp1 = 1.
pFuelPowLp = 1.
#
#
pSteamUseRatePcc = 1.
pPowUseRatePcc = 1.

pSteamUseRateDacAir = 1.
pSteamUseRateDacFlue = 1.

pPowUseRateDacAir = 1.
pPowUseRateDacFlue = 1.

pHeatRateSteamTur = 1.

pCo2PccCapRate = 1.
pSorbCo2Cap = 1.

pPowUseRateDac = 1.

pPowUseRateComp = 1.
#

pCostInvCombTurb = 1e+02
pCostInvSteaTurb = 1e+02
pCostInvTransInter = 1e+02
pCostInvPcc = 1e+02
pCostInvDac = 1e+03
pCostInvComp = 1e+01

# Cost parameters.
pCostFuel = 100.
pEmissionPrice = 1.
pCo2TranspPrice = 0.1

# Constraints
# 1
@constraint(m, capCombTurIn[i = 1:tHorz], vPowCombTurb[i] <= vCapCombTurb)
# 2
@constraint(m, capSteamTurbIn[i = 1:tHorz], ePowSteamTurb[i] <= vCapSteamTurb)
# 3
@constraint(m, capTransmIntercIn[i = 1:tHorz], ePowOut[i] <= vCapTransInter)
# 4
@constraint(m, fuelEq[i = 1:tHorz], eFuel[i] == pHeatRateCombTur * vPowCombTurb[i])
# 5
@constraint(m, co2FuelEq[i = 1:tHorz], eCo2Fuel[i] == pEmissFactor * vPowCombTurb[i])
# 6
@constraint(m, powOutEq[i = 1:tHorz], ePowOut[i] == vPowCombTurb[i] + ePowSteamTurb[i] - ePowUseDac[i] - ePowUsePcc[i] - ePowUseComp[i])

# Steam
# 7a
@constraint(m, powHpEq[i = 1:tHorz], ePowHp[i] == pFuelPowHp * eFuel[i])
# 7b
@constraint(m, powIpEq[i = 1:tHorz], ePowIp[i] == pFuelPowIp * eFuel[i])

# 8
@constraint(m, stamIpLpEq[i = 1:tHorz], eSteamIpLp[i] == pFuelSteamIpLp * eFuel[i])
# 9
@constraint(m, steamLp1[i = 1:tHorz], eSteamLp1[i] == eSteamIpLp[i] - eSteamUsePcc[i] - eSteamUseDac[i])
# 10
@constraint(m, powLp1Eq[i = 1:tHorz], ePowLp1[i] == eSteamLp1[i] / pHeatRateLp1)
# 11
@constraint(m, powLp2Eq[i = 1:tHorz], ePowLp2[i] == pFuelPowLp * eFuel[i])
# 12
@constraint(m, powerSteamEq[i = 1:tHorz], ePowSteamTurb[i] == ePowHp[i] + ePowIp[i] + ePowLp1[i] + ePowLp2[i])

# PCC
# 13
@constraint(m, co2CapPccEq[i = 1:tHorz], eCo2CapPcc[i] == pCo2PccCapRate * eCo2Fuel[i])
# 14
@constraint(m, co2PccOutEq[i = 1:tHorz], eCo2PccOut[i] == eCo2Fuel[i] - eCo2CapPcc[i])
# 15
@constraint(m, co2DacFlueInEq[i = 1:tHorz], eCo2DacFlueIn[i] == eCo2PccOut[i] - vCo2PccVent[i])
# 16
@constraint(m, co2CapPccIn[i = 1:tHorz], eCo2CapPcc[i] <= vCapPcc)
# 17
@constraint(m, steamUsePccEq[i = 1:tHorz], eSteamUsePcc[i] == pSteamUseRatePcc * eCo2CapPcc[i])
# 18
@constraint(m, powerUsePccEq[i = 1:tHorz], ePowUsePcc[i] == pPowUseRatePcc * eCo2CapPcc[i])

# DAC
# v2 of the model.
@constraint(m, a0Eq[i = 1:tHorz], a0[i] == vAbs[i])
@constraint(m, aR0Eq[i = 1:tHorz], aR0[i] == vReg[i])
@constraint(m, a1dEq[i = 0:tHorz - 1], a1[i + 1] == a0[i])
@constraint(m, storeFeq[i = 0:tHorz - 1], sF[i + 1] == sF[i] - vAbs[i] + aR0[i])
@constraint(m, storeSeq[i = 0:tHorz - 1], sS[i + 1] == sS[i] - vReg[i] + a1[i])
# 21
#@constraint(m, capDacEq, sF[0] + sS[0] == vCapDac / pSorbCo2Cap)
@constraint(m, capDacEq, sF[0] == vCapDac / pSorbCo2Cap)
# 22
@constraint(m, icA1Eq, a1[0] == 0.)
# @constraint(m, icA2Eq, aR1[0] == 0.)
#@constraint(m, icSfEq, sF[0] == 100.)
@constraint(m, icSsEq, sS[0] == 0.)
# 23
@constraint(m, co2StoreDacEq[i = 1:tHorz], eCo2StoreDac[i] == pSorbCo2Cap * sS[i])
# 24
@constraint(m, co2CapDacEq[i = 1:tHorz], eCo2CapDac[i] == pSorbCo2Cap * aR0[i])
# 25
@constraint(m, steamUseDacEq[i = 1:tHorz], eSteamUseDac[i] == pSteamUseRateDacFlue * eCo2CapDac[i])
# 26
@constraint(m, powUseDacEq[i = 1:tHorz], ePowUseDac[i] == pPowUseRateDacFlue * eCo2CapDac[i])
# equal to the amount vented at least in flue mode.
@constraint(m, co2DacFlueVentEq[i = 1:tHorz], eCo2DacVent[i] == eCo2DacFlueIn[i] - eCo2CapDac[i])

# Co2 Compression
# 27
@constraint(m, co2CompEq[i = 1:tHorz], eCo2Comp[i] == eCo2CapPcc[i] + eCo2CapDac[i])
# 28
@constraint(m, powUseCompEq[i = 1:tHorz], ePowUseComp[i] == pPowUseRateComp * eCo2Comp[i])
# 29
@constraint(m, powUseCompIn[i = 1:tHorz], ePowUseComp[i] <= vCapComp)

@constraint(m, co2VentEq[i = 1:tHorz], eCo2Vent[i] == vCo2PccVent[i] + eCo2DacVent[i])

rnum = rand(Float64, tHorz)

@constraint(m, epwIn[i = 1:tHorz], ePowOut[i] >= 1 * rnum[i])
#@expression(m, eObjfExpr, pCostInvCombTurb * vCapCombTurb + pCostInvSteaTurb * vCapSteamTurb
#                           + pCostInvTransInter * vCapTransInter + pCostInvPcc * vCapPcc
#                           + pCostInvDac * vCapDac + pCostInvComp * vCapComp
#                           + sum(pCostFuel * eFuel[i] + pEmissionPrice * eCo2Vent[i] + pCo2TranspPrice * eCo2Comp[i] for i in 1:tHorz))
@expression(m, eObjfExpr, sum(pCostFuel * eFuel[i] + pEmissionPrice * eCo2Vent[i] + pCo2TranspPrice * eCo2Comp[i] for i in 1:tHorz))
#@objective(m, Min, eobjf)
@objective(m, Min, eObjfExpr)

set_optimizer(m, Clp.Optimizer)
set_optimizer_attribute(m, "LogLevel", 3)
set_optimizer_attribute(m, "PresolveType", 1)

# set_start_value(vCapCombTurb, 3.)
# set_start_value(vCapSteamTurb, 2.)
# set_start_value(vCapTransInter, 5.)
# set_start_value(vCapPcc, 5.)
# set_start_value(vCapDac, 4.)
# set_start_value(vCapComp, 8.)

fix(vCapCombTurb, 3., force = true)
fix(vCapSteamTurb, 2., force = true)
fix(vCapTransInter, 5., force = true)
fix(vCapPcc, 5., force = true)
fix(vCapDac, 4., force = true)
fix(vCapComp, 8., force = true)

fix(a0[0], 0., force = true)
fix(aR0[0], 0., force = true)
fix(vAbs[0], 0., force = true)
fix(vReg[0], 0., force = true)

optimize!(m)
println(termination_status(m))

println(value.(ePowOut)[:])
println(value.(sF)[:])
println(value.(sS)[:])
println(value.(a1)[:])
f = open("model.lp", "w")
print(f, m)
close(f)

write_to_file(
    m,
    "lp_mk0.mps";
    #format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC
)

# Design decisions.
# vCapCombTurb
# vCapSteamTurb
# vCapTransInter
# vCapPcc
# vCapDac
# vCapComp
#
# Raw materials.
# eFuel
# sF0
# sS0
#
# penalties
# vCo2Vent

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
function parse_line(s)
    m = match(r"(\d+)-(\d+) (.): (.*)$", s).captures
    return (min = parse(Int, m[1]), max = parse(Int, m[2]), letter = m[3][begin], pwd = m[4])
end

records = parse_line.(readlines("day02.txt"))

using StatsBase

# Part1
# Find number of valid passwords with at least and at most of a letter
count(r -> r.min <= get(countmap(r.pwd), r.letter, -1) <= r.max, records)

# Part2
# Find number of valid passwords with the letter at either min/max but not both
count(r -> xor(r.pwd[r.min] == r.letter, r.pwd[r.max] == r.letter) , records)


text:
module TokyoVariantAnalysis4C2
using Logging
using DataFrames
using Dates
using Plots
using Plots.PlotMeasures

using ..TokyoVariantAnalysis4
using ..TokyoVariantAnalysis4A
using ..TokyoVariantAnalysis4B2

@info "--------"
@info "プロット・バックエンド・設定"
pyplot(
    titlefont=font("Meiryo",10),
    guidefont=font("Meiryo",10),
    tickfont=font("Meiryo",9)
)
const plot_date_start = Date("2021-05-25")
const plot_date_end   = plot_date_start + Day(100)
const date_ticks_dict = Dict(
    Date("2021-03-21") => "緊急事態宣言解除",
    Date("2021-04-01") => "2021\u201004\u201001",
    Date("2021-04-15") => "04\u201015",
    Date("2021-04-25") => "緊急事態宣言発令",
    Date("2021-05-01") => "2021\u201005\u201001",
    Date("2021-05-15") => "05\u201015",
    Date("2021-06-01") => "2021\u201006\u201001",
    Date("2021-06-15") => "06\u201015",
    Date("2021-06-20") => "緊急事態宣言解除",
    Date("2021-07-01") => "2021\u201007\u201001",
    Date("2021-07-04") => "都議会議員選挙",
    Date("2021-07-12") => "緊急事態宣言発令",
    Date("2021-07-15") => "07\u201015",
    Date("2021-07-23") => "オリンピック開会",
    Date("2021-08-01") => "2021\u201008\u201001",
    Date("2021-08-08") => "オリンピック閉会",
    Date("2021-08-15") => "08\u201015",
    Date("2021-08-22") => "宣言解除予定",
    Date("2021-09-01") => "2021\u201009\u201001",
)
const date_ticks_dates = sort(collect(keys(date_ticks_dict)))
const date_ticks_dict_en = Dict(
    Date("2021-03-21") => "End\\ of\\ the\\ SoE",
    Date("2021-04-01") => "2021\u201004\u201001",
    Date("2021-04-15") => "04\u201015",
    Date("2021-04-25") => "Decl.\\ of\\ the\\ SoE",
    Date("2021-05-01") => "2021\u201005\u201001",
    Date("2021-05-15") => "05\u201015",
    Date("2021-06-01") => "2021\u201006\u201001",
    Date("2021-06-15") => "06\u201015",
    Date("2021-06-20") => "End\\ of\\ the\\ SoE",
    Date("2021-07-01") => "2021\u201007\u201001",
    Date("2021-07-04") => "Assembly\\ election",
    Date("2021-07-12") => "Decl.\\ of\\ the\\ SoE",
    Date("2021-07-15") => "07\u201015",
    Date("2021-07-23") => "Olympics\\ starts",
    Date("2021-08-01") => "2021\u201008\u201001",
    Date("2021-08-08") => "Olympics\\ ends",
    Date("2021-08-15") => "08\u201015",
    Date("2021-08-22") => "SoE\\ ends\\ (sched)",
    Date("2021-09-01") => "2021\u201009\u201001",
)
const date_ticks_dates_en = sort(collect(keys(date_ticks_dict_en)))

plots = Dict{String,Plots.Plot}()
@info "--------"
@info "L452R 割合尤度プロット"
for once in 1:1
    pl = plot(
        size=(640,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=-32px,
        title="東京（健安研調査・民間検査機関等調査）"
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 0.0, 1.0)
    ylabel!(pl, "L452R 変異のある変異株（主に δ 株）の割合")
    yticks!(pl,
        [0.0, 0.2, 0.4, 0.6, 0.8, 1.0],
        ["0\\%", "20\\%", "40\\%", "60\\%", "80\\%", "100\\%"]
    )
    # 尤度ヒートマップ描画
    hm_xrange = range(hm_x_min, hm_x_max; length=hm_x_size)
    hm_yrange = range(hm_y_min, hm_y_max; length=hm_y_size)
    heatmap!(pl, hm_xrange, hm_yrange, hm_f;
        c=:Wistia, colorbar_title="相対尤度（最大尤度を 1 とする）", colorbar_titlefontfamily="Meiryo")
    plot!(pl, (t)->1.0/(1.0+1.0/exp(likelihood_max.α+likelihood_max.β*t)), lc=:red3, ls=:dot)
    # データ点
    pname = "L452R_prop"
    lmtl = -100.0
    lmtu = +100.0
    lmt(r) = min(max(r, lmtl), lmtu)
    ds  = variants_df.date_value
    #des = df.date_error
    rs  = variants_df[!,Symbol("$(pname)")]
    ls  = variants_df[!,Symbol("$(pname)_cil")]
    us  = variants_df[!,Symbol("$(pname)_ciu")]
    n_rows = nrow(variants_df)
    pm_ds  = Vector{Float64}(); sizehint!(pm_ds,  n_rows)
    #pm_des = Vector{Float64}(); sizehint!(pm_des, n_rows)
    pm_rs  = Vector{Float64}(); sizehint!(pm_rs,  n_rows)
    pm_els = Vector{Float64}(); sizehint!(pm_els, n_rows)
    pm_eus = Vector{Float64}(); sizehint!(pm_eus, n_rows)
    ps_ds  = Vector{Float64}(); sizehint!(ps_ds,  n_rows)
    #ps_des = Vector{Float64}(); sizehint!(ps_des, n_rows)
    ps_rs  = Vector{Float64}(); sizehint!(ps_rs,  n_rows)
    ps_els = Vector{Float64}(); sizehint!(ps_els, n_rows)
    ps_eus = Vector{Float64}(); sizehint!(ps_eus, n_rows)
    for i in 1:n_rows
        if !ismissing(rs[i])
            r  = lmt(rs[i])
            el = r - lmt(ls[i])
            eu = lmt(us[i]) - r
            if TokyoVariantAnalysis4A.lh_sdate_value ≤ ds[i] ≤ TokyoVariantAnalysis4A.lh_edate_value
                push!(pm_ds,  ds[i])
                #push!(pm_des, des[i])
                push!(pm_rs,  r)
                push!(pm_els, el)
                push!(pm_eus, eu)
            else
                push!(ps_ds,  ds[i])
                #push!(ps_des, des[i])
                push!(ps_rs,  r)
                push!(ps_els, el)
                push!(ps_eus, eu)
            end
        end
    end
    scatter!(pl, pm_ds, pm_rs; yerrors=(pm_els, pm_eus),
        color=RGB(0.0,0.5,1.0), msc=RGB(0.4,0.4,0.4), legend=:none)
    scatter!(pl, ps_ds, ps_rs; yerrors=(ps_els, ps_eus),
        color=RGB(0.0,0.5,1.0), msc=RGB(0.7,0.7,0.7), legend=:none)
    # グリッド
    vline!(pl, vdates; lc=RGBA(0.0,0.0,0.0,0.1), ls=:dot)
    hline!(pl, [0.2, 0.4, 0.6, 0.8]; lc=RGBA(0.0,0.0,0.0,0.1), ls=:dot)
    # アノテーション
    #coef = ratio_fit_coefs[pname][2]
    annotate!(pl, date_to_value(Date("2021-05-28")), 0.925, text("""
        「東京都新型コロナウイルス感染症対策本部報」より
        民間検査機関等における L452R の検査が報告されるようになった 2021-06-01 以降のデータ
        ただし、デルタ株クラスターが含まれているとみられる場合は、クラスターあたり 1 人分のみ計上
        赤点線は指数関数モデルでの最尤線. 背景の色は最大尤度に相対的な尤度を表す""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left)
    )
    plname = "L452R"
    plots[plname] = pl
    @info plname
end
@info "--------"
@info "in English"
for once in 1:1
    pl = plot(
        size=(640,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=-32px,
        title="Tokyo (Inst. of Public Health and Private Labs.)"
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates_en)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict_en[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 0.0, 1.0)
    ylabel!(pl, "Proportion of the L452R mutation (mainly the δ variant)")
    yticks!(pl,
        [0.0, 0.2, 0.4, 0.6, 0.8, 1.0],
        ["0\\%", "20\\%", "40\\%", "60\\%", "80\\%", "100\\%"]
    )
    # 尤度ヒートマップ描画
    hm_xrange = range(hm_x_min, hm_x_max; length=hm_x_size)
    hm_yrange = range(hm_y_min, hm_y_max; length=hm_y_size)
    heatmap!(pl, hm_xrange, hm_yrange, hm_f;
        c=:Wistia, colorbar_title="relative likelihood (the maximum likelihood is 1)", colorbar_titlefontfamily="Meiryo")
    plot!(pl, (t)->1.0/(1.0+1.0/exp(likelihood_max.α+likelihood_max.β*t)), lc=:red3, ls=:dot)
    # データ点
    pname = "L452R_prop"
    lmtl = -100.0
    lmtu = +100.0
    lmt(r) = min(max(r, lmtl), lmtu)
    ds  = variants_df.date_value
    #des = df.date_error
    rs  = variants_df[!,Symbol("$(pname)")]
    ls  = variants_df[!,Symbol("$(pname)_cil")]
    us  = variants_df[!,Symbol("$(pname)_ciu")]
    n_rows = nrow(variants_df)
    pm_ds  = Vector{Float64}(); sizehint!(pm_ds,  n_rows)
    #pm_des = Vector{Float64}(); sizehint!(pm_des, n_rows)
    pm_rs  = Vector{Float64}(); sizehint!(pm_rs,  n_rows)
    pm_els = Vector{Float64}(); sizehint!(pm_els, n_rows)
    pm_eus = Vector{Float64}(); sizehint!(pm_eus, n_rows)
    ps_ds  = Vector{Float64}(); sizehint!(ps_ds,  n_rows)
    #ps_des = Vector{Float64}(); sizehint!(ps_des, n_rows)
    ps_rs  = Vector{Float64}(); sizehint!(ps_rs,  n_rows)
    ps_els = Vector{Float64}(); sizehint!(ps_els, n_rows)
    ps_eus = Vector{Float64}(); sizehint!(ps_eus, n_rows)
    for i in 1:n_rows
        if !ismissing(rs[i])
            r  = lmt(rs[i])
            el = r - lmt(ls[i])
            eu = lmt(us[i]) - r
            if TokyoVariantAnalysis4A.lh_sdate_value ≤ ds[i] ≤ TokyoVariantAnalysis4A.lh_edate_value
                push!(pm_ds,  ds[i])
                #push!(pm_des, des[i])
                push!(pm_rs,  r)
                push!(pm_els, el)
                push!(pm_eus, eu)
            else
                push!(ps_ds,  ds[i])
                #push!(ps_des, des[i])
                push!(ps_rs,  r)
                push!(ps_els, el)
                push!(ps_eus, eu)
            end
        end
    end
    scatter!(pl, pm_ds, pm_rs; yerrors=(pm_els, pm_eus),
        color=RGB(0.0,0.5,1.0), msc=RGB(0.4,0.4,0.4), legend=:none)
    scatter!(pl, ps_ds, ps_rs; yerrors=(ps_els, ps_eus),
        color=RGB(0.0,0.5,1.0), msc=RGB(0.7,0.7,0.7), legend=:none)
    # グリッド
    vline!(pl, vdates; lc=RGBA(0.0,0.0,0.0,0.1), ls=:dot)
    hline!(pl, [0.2, 0.4, 0.6, 0.8]; lc=RGBA(0.0,0.0,0.0,0.1), ls=:dot)
    # アノテーション
    #coef = ratio_fit_coefs[pname][2]
    annotate!(pl, date_to_value(Date("2021-05-28")), 0.925, text("""
        Data source: Daily reports of “Tōkyō-to Shingata Corona Virus Kansenshō Taisaku Hombu”
        Using data since 2021-06-01 when private laboratories began reporting L452R tests.
        Only one case included per cluster if case numbers contain Delta clusters.
        The red dotted line is the maximum likelihood line for the exponential model.
        The background color indicates the likelihood relative to the maximum.""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left)
    )
    plname = "L452R_en"
    plots[plname] = pl
    @info plname
end

end # module

max_stars_repo_name: muimuiz/TokyoVariant

max_issues_repo_name: muimuiz/TokyoVariant

max_forks_repo_name: muimuiz/TokyoVariant


text:
# plan-psf.jl

export PlanPSF, plan_psf
import AbstractFFTs
import FFTW

"""
    PlanPSF{T,Tf,Ti}( ; nx::Int, nz::Int, px::Int, pz::Int, T::DataType)
Struct for storing work arrays and factors for 2D convolution for one thread.
Each PSF is `px × pz`
- `T` datatype of work arrays (subtype of `AbstractFloat`)
- `nx::Int = 128` (`ny` implicitly the same)
- `nz::Int = nx` image size is `[nx,nx,nz]`
- `px::Int = 1`
- `pz::Int = px` (PSF size)
- `padsize::Tuple` : `(padup, paddown, padleft, padright)`
- `workmat [nx+padup+paddown, nz+padleft+padright]` 2D padded image for FFT convolution
- `workvecx [nx+padup+paddown,]`: 1D work vector
- `workvecz [nz+padleft+padright,]`: 1D work vector
- `img_compl [nx+padup+paddown, nz+padleft+padright]`: 2D [complex] padded image for FFT
- `ker_compl [nx+padup+paddown, nz+padleft+padright]`: 2D [complex] padded image for FFT
- `fft_plan::Tf` plan for doing FFT; see `plan_fft!`
- `ifft_plan::Ti` plan for doing IFFT; see `plan_ifft!`
"""
struct PlanPSF{T, Tf, Ti}
    nx::Int
    nz::Int
    px::Int
    pz::Int
    padsize::NTuple{4, Int}
    workmat::Matrix{T}
    workvecx::Vector{T}
    workvecz::Vector{T}
    img_compl::Matrix{Complex{T}}
    ker_compl::Matrix{Complex{T}}
    fft_plan::Tf # Union{AbstractFFTs.ScaledPlan, FFTW.cFFTWPlan}
    ifft_plan::Ti # Union{AbstractFFTs.ScaledPlan, FFTW.cFFTWPlan}

    function PlanPSF( ;
        nx::Int = 128,
        nz::Int = nx,
        px::Int = 1,
        pz::Int = px,
        T::DataType = Float32,
    )

        T <: AbstractFloat || throw("invalid T=$T")
        padup = _padup(nx, px)
        paddown = _paddown(nx, px)
        padleft = _padleft(nz, pz)
        padright = _padright(nz, pz)
        padsize = (padup, paddown, padleft, padright)

        workmat = Matrix{T}(undef, nx+padup+paddown, nz+padleft+padright)
        workvecx = Vector{T}(undef, nx+padup+paddown)
        workvecz = Vector{T}(undef, nz+padleft+padright)

        # complex padimg
        img_compl = Matrix{Complex{T}}(undef, nx+padup+paddown, nz+padleft+padright)
        # complex kernel
        ker_compl = Matrix{Complex{T}}(undef, nx+padup+paddown, nz+padleft+padright)

        fft_plan = plan_fft!(ker_compl)
        ifft_plan = plan_ifft!(ker_compl)
        Tf = typeof(fft_plan)
        Ti = typeof(ifft_plan)

        new{T, Tf, Ti}(
            nx,
            nz,
            px,
            pz,
            padsize,
            workmat,
            workvecx,
            workvecz,
            img_compl,
            ker_compl,
            fft_plan,
            ifft_plan,
        )
    end
end

"""
    plan_psf( ; nx::Int, nz::Int, px::Int, pz::Int, nthread::Int, T::DataType)
Make Vector of structs for storing work arrays and factors
for 2D convolution with SPECT depth-dependent PSF model,
threaded across planes parallel to detector.
Option
- `nx::Int = 128`
- `nz::Int = nx`
- `px::Int = 1`
- `pz::Int = px` PSF size is `px × pz`
- `T` : datatype of work arrays, defaults to `Float32`
- `nthread::Int` # of threads, defaults to `Threads.nthreads()`
"""
function plan_psf( ;
    nx::Int = 128,
    nz::Int = nx,
    px::Int = 1,
    pz::Int = px,
    nthread::Int = Threads.nthreads(),
    T::DataType = Float32,
)
    return [PlanPSF( ; nx, nz, px, pz, T) for id = 1:nthread]
end

"""
    show(io::IO, ::MIME"text/plain", plan::PlanPSF)
"""
function Base.show(io::IO, ::MIME"text/plain", plan::PlanPSF{T}) where {T}
    t = typeof(plan)
    println(io, t)
    for f in (:nx, :nz, :px, :pz, :padsize)
        p = getproperty(plan, f)
        t = typeof(p)
        println(io, " ", f, "::", t, " ", p)
    end
    for f in (:workmat, :workvecx, :workvecz, :img_compl, :ker_compl, :fft_plan, :ifft_plan)
        p = getproperty(plan, f)
        println(io, " ", f, ":", " ", summary(p))
    end
    println(io, " (", sizeof(plan), " bytes)")
end

"""
    show(io::IO, mime::MIME"text/plain", vp::Vector{<:PlanPSF})
"""
function Base.show(io::IO, mime::MIME"text/plain", vp::Vector{<: PlanPSF})
    t = typeof(vp)
    println(io, length(vp), "-element ", t, " with N=", vp[1].nx)
#   show(io, mime, vp[1])
end

"""
    sizeof(::PlanPSF)
Show size in bytes of `PlanPSF` object.
"""
function Base.sizeof(ob::T) where {T <: PlanPSF}
    sum(f -> sizeof(getfield(ob, f)), fieldnames(typeof(ob)))
end