text:
"""
	loglikelihood_contactratemodel_ct(contactrate,θ,model::KenyaSerology.CoVAreaModel)

Calculate the log-posterior density for `contactrate` for parameters fixed as `θ` and data from `model`.
"""
function loglikelihood_contactratemodel_ct(contactrate,θ,model::KenyaSerology.CoVAreaModel,λ)
    @unpack R,E₀,I₀,χ,M_PCR,α,p_test,P_eff = θ #Parameters to fit 1) initial conditions and R₀ 2) parameters for Neg. bin.
    @unpack PCR_cases,sero_cases,N,σ,γ,prob,PCR_array,sero_array,log_priors,PCR_sensitivity,PCR_specificity,sero_specificity,sero_sensitivity,M_BB,relative_testing_rate = model
    T = eltype(contactrate)
    N_eff = N*P_eff
    u0 = convert.(T, [N_eff,E₀,I₀,0.,0.])
    p = convert.(T,vcat([R,σ,γ,N_eff],contactrate))
	smoothness_penalty = T(0.)
	LL = T(0.)

    try
        sol = solve(prob, BS3(); u0=u0, p=p, saveat = 1)
        incidence = KenyaSerology.get_incidence(sol)
        num_PCR = KenyaSerology.simple_conv(incidence,PCR_array)
        p_sero = sero_sensitivity*KenyaSerology.simple_conv(incidence,sero_array)/N
        for t in 59:min(length(incidence),size(PCR_cases,1))
            d_PCR_pos = PCR_cases[t,1]#Detected number of PCR pos that day
            d_PCR_neg = PCR_cases[t,2]#Detected number of PCR neg that day, -1 === no negative tests available
            d_sero_pos = 0
            n_sero = 0
            if t<= size(sero_cases,1)
                d_sero_pos = sero_cases[t,1] #Detected number of sero pos that day
                n_sero = sero_cases[t,1] + sero_cases[t,2] #Number sero tested on that day
            end
            if d_PCR_neg >= 0
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                p_PCR = χ*num_PCR[t]/((χ-1)*num_PCR[t] + N)
                LL += logpdf(BetaBinomial(d_PCR_pos+d_PCR_neg,p_PCR*M_PCR,(1-p_PCR)*M_PCR),d_PCR_pos)#likelihood contribution from PCR testing --- given
            end
            if d_PCR_neg == -1
                #Convert from μ,α parameterisation to p,r parameterisation
                μ = relative_testing_rate[t]*p_test*num_PCR[t] + 0.001
                σ² = μ + α*μ^2
                p_negbin = 1 - (α*μ^2/σ²)
                r_negbin = 1/α
                LL += logpdf(NegativeBinomial(r_negbin,p_negbin),d_PCR_pos)#likelihood contribution from PCR testing
            end
            #Convert from the p_hat,M_BB parameterisation of the Betabinomial distribution to the α_BB,β_BB parameterisation
            if n_sero > 0
                p_hat = p_sero[t] + (1- p_sero[t])*(1-sero_specificity)
                α_BB = M_BB*p_hat
                β_BB = M_BB*(1-p_hat)
                LL += logpdf(BetaBinomial(n_sero,α_BB,β_BB),d_sero_pos)#Likelihood contribution from sero testing
            end
        end
        # return LL + log_priors(θ)
		for t in 2:(length(contactrate)-1)
			smoothness_penalty += (contactrate[t-1] + contactrate[t+1] - 2*contactrate[t])^2
		end
		return LL + log_priors(θ) - (λ*smoothness_penalty/length(contactrate))

    catch
        return T(-Inf)
    end
end

"""
        function loglikelihood_contactratemodel_ct(contactrate,θ,model::KenyaSerology.CoVAreaModel,λ,p_waning)

Calculate the log-posterior density for `contactrate` for parameters fixed as `θ` and data from `model` with serological waning rate `p_waning`.
"""
function loglikelihood_contactratemodel_ct(contactrate,θ,model::KenyaSerology.CoVAreaModel,λ,p_waning)
    @unpack R,E₀,I₀,χ,M_PCR,α,p_test,P_eff = θ #Parameters to fit 1) initial conditions and R₀ 2) parameters for Neg. bin.
    @unpack PCR_cases,sero_cases,N,σ,γ,prob,PCR_array,sero_array,log_priors,PCR_sensitivity,PCR_specificity,sero_specificity,sero_sensitivity,M_BB,relative_testing_rate = model
    T = eltype(contactrate)
    N_eff = N*P_eff
    u0 = convert.(T, [N_eff,E₀,I₀,0.,0.])
    p = convert.(T,vcat([R,σ,γ,N_eff],contactrate))
    smoothness_penalty = T(0.)

    for t = 31:length(sero_array)
        sero_array[t] = (1-p_waning)^(t-30)
    end
    LL = T(0.)

    try
        sol = solve(prob, BS3(); u0=u0, p=p, saveat = 1)
        incidence = KenyaSerology.get_incidence(sol)
        num_PCR = KenyaSerology.simple_conv(incidence,PCR_array)
        p_sero = sero_sensitivity*KenyaSerology.simple_conv(incidence,sero_array)/N
        for t in 59:min(length(incidence),size(PCR_cases,1))
            d_PCR_pos = PCR_cases[t,1]#Detected number of PCR pos that day
            d_PCR_neg = PCR_cases[t,2]#Detected number of PCR neg that day, -1 === no negative tests available
            d_sero_pos = 0
            n_sero = 0
            if t<= size(sero_cases,1)
                d_sero_pos = sero_cases[t,1] #Detected number of sero pos that day
                n_sero = sero_cases[t,1] + sero_cases[t,2] #Number sero tested on that day
            end
            if d_PCR_neg >= 0
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                p_PCR = χ*num_PCR[t]/((χ-1)*num_PCR[t] + N)
                LL += logpdf(BetaBinomial(d_PCR_pos+d_PCR_neg,p_PCR*M_PCR,(1-p_PCR)*M_PCR),d_PCR_pos)#likelihood contribution from PCR testing --- given
            end
            if d_PCR_neg == -1
                #Convert from μ,α parameterisation to p,r parameterisation
                μ = relative_testing_rate[t]*p_test*num_PCR[t] + 0.001
                σ² = μ + α*μ^2
                p_negbin = 1 - (α*μ^2/σ²)
                r_negbin = 1/α
                LL += logpdf(NegativeBinomial(r_negbin,p_negbin),d_PCR_pos)#likelihood contribution from PCR testing
            end
            #Convert from the p_hat,M_BB parameterisation of the Betabinomial distribution to the α_BB,β_BB parameterisation
            if n_sero > 0
                p_hat = p_sero[t] + (1- p_sero[t])*(1-sero_specificity)
                α_BB = M_BB*p_hat
                β_BB = M_BB*(1-p_hat)
                LL += logpdf(BetaBinomial(n_sero,α_BB,β_BB),d_sero_pos)#Likelihood contribution from sero testing
            end
        end
        # return LL + log_priors(θ)
		for t in 2:(length(contactrate)-1)
			smoothness_penalty += (contactrate[t-1] + contactrate[t+1] - 2*contactrate[t])^2
		end

		return LL + log_priors(θ) - (λ*smoothness_penalty/length(contactrate))

    catch
        return T(-Inf)
    end
end

"""
        function loglikelihood_contactratemodel_ct(contactrate,θ,model::KenyaSerology.CoVAreaModel,λ,p_waning,σ²_aboveone)

Calculate the log-posterior density for `contactrate` for parameters fixed as `θ` and data from `model` with serological waning rate `p_waning`.
`σ²_aboveone` scales the penalty for the Ct value going above 1.
"""
function loglikelihood_contactratemodel_ct(contactrate,θ,model::KenyaSerology.CoVAreaModel,λ,p_waning,σ²_aboveone)
    @unpack R,E₀,I₀,χ,M_PCR,α,p_test,P_eff = θ #Parameters to fit 1) initial conditions and R₀ 2) parameters for Neg. bin.
    @unpack PCR_cases,sero_cases,N,σ,γ,prob,PCR_array,sero_array,log_priors,PCR_sensitivity,PCR_specificity,sero_specificity,sero_sensitivity,M_BB,relative_testing_rate = model
    T = eltype(contactrate)
    N_eff = N*P_eff
    u0 = convert.(T, [N_eff,E₀,I₀,0.,0.])
    p = convert.(T,vcat([R,σ,γ,N_eff],contactrate))
    smoothness_penalty = T(0.)
	above_one_penalty = T(0.)

    for t = 31:length(sero_array)
        sero_array[t] = (1-p_waning)^(t-30)
    end
    LL = T(0.)

    try
        sol = solve(prob, BS3(); u0=u0, p=p, saveat = 1)
        incidence = KenyaSerology.get_incidence(sol)
        num_PCR = KenyaSerology.simple_conv(incidence,PCR_array)
        p_sero = sero_sensitivity*KenyaSerology.simple_conv(incidence,sero_array)/N
        for t in 59:min(length(incidence),size(PCR_cases,1))
            d_PCR_pos = PCR_cases[t,1]#Detected number of PCR pos that day
            d_PCR_neg = PCR_cases[t,2]#Detected number of PCR neg that day, -1 === no negative tests available
            d_sero_pos = 0
            n_sero = 0
            if t<= size(sero_cases,1)
                d_sero_pos = sero_cases[t,1] #Detected number of sero pos that day
                n_sero = sero_cases[t,1] + sero_cases[t,2] #Number sero tested on that day
            end
            if d_PCR_neg >= 0
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                p_PCR = χ*num_PCR[t]/((χ-1)*num_PCR[t] + N)
                LL += logpdf(BetaBinomial(d_PCR_pos+d_PCR_neg,p_PCR*M_PCR,(1-p_PCR)*M_PCR),d_PCR_pos)#likelihood contribution from PCR testing --- given
            end
            if d_PCR_neg == -1
                #Convert from μ,α parameterisation to p,r parameterisation
                μ = relative_testing_rate[t]*p_test*num_PCR[t] + 0.001
                σ² = μ + α*μ^2
                p_negbin = 1 - (α*μ^2/σ²)
                r_negbin = 1/α
                LL += logpdf(NegativeBinomial(r_negbin,p_negbin),d_PCR_pos)#likelihood contribution from PCR testing
            end
            #Convert from the p_hat,M_BB parameterisation of the Betabinomial distribution to the α_BB,β_BB parameterisation
            if n_sero > 0
                p_hat = p_sero[t] + (1- p_sero[t])*(1-sero_specificity)
                α_BB = M_BB*p_hat
                β_BB = M_BB*(1-p_hat)
                LL += logpdf(BetaBinomial(n_sero,α_BB,β_BB),d_sero_pos)#Likelihood contribution from sero testing
            end
        end
        # return LL + log_priors(θ)
		for t in 2:(length(contactrate)-1)
			smoothness_penalty += (contactrate[t-1] + contactrate[t+1] - 2*contactrate[t])^2
		end
		for t in 1:(length(contactrate))
			if contactrate[t] > 1
				above_one_penalty += (log(contactrate[t])^2)/(2*σ²_aboveone) #Log-normal assumption above Ct = 1,0.174 chosen so that < 50% greater contacts are 99% probable
			end
		end
		return LL + log_priors(θ) - (λ*smoothness_penalty/length(contactrate)) - above_one_penalty

    catch
        return T(-Inf)
    end
end

"""
	C_t_neg_log_density_grad!(grad,ct,θ,model::KenyaSerology.CoVAreaModel,log_prior_ct::Function,baseline_contact)

In-place calculation of gradient of contact rate negative log-posterior density (first `length(baseline_contact)` days are fixed as Google estimate).
"""
function C_t_neg_log_density_grad!(grad,ct,θ,model::KenyaSerology.CoVAreaModel,baseline_contact,λ)
        l = x -> -loglikelihood_contactratemodel_ct(vcat(baseline_contact,x),θ,model,λ)
        ForwardDiff.gradient!(grad,l,ct)
        return nothing
end

"""
        function C_t_neg_log_density_grad!(grad,ct,θ,model::KenyaSerology.CoVAreaModel,baseline_contact,λ,p_waning)

In-place calculation of gradient of contact rate negative log-posterior density (first `length(baseline_contact)` days are fixed as Google estimate).
With sero-waning `p_waning`.
"""
function C_t_neg_log_density_grad!(grad,ct,θ,model::KenyaSerology.CoVAreaModel,baseline_contact,λ,p_waning)
        l = x -> -loglikelihood_contactratemodel_ct(vcat(baseline_contact,x),θ,model,λ,p_waning)
        ForwardDiff.gradient!(grad,l,ct)
        return nothing
end

"""
	ADAM_step!(X_n::Vector{Float64},
						∇V::Vector{Float64},
						∇V_n::Vector{Float64},
						∇V_sq_n::Vector{Float64},
						n::Int64,
						ADAM_params::NamedTuple)

In-place calculation of parameter, gradient and squared gradient estimates using ADAM SGD method.
"""
function ADAM_step!(X_n::Vector{Float64},
					∇V::Vector{Float64},
					∇V_n::Vector{Float64},
					∇V_sq_n::Vector{Float64},
					hat_∇V_n::Vector{Float64},
					hat_∇V_sq_n::Vector{Float64},
					n::Int64,
					ADAM_params::NamedTuple)
	@unpack β₁,β₂,η,ϵ = ADAM_params
	∇V_n .= β₁.*∇V_n .+ (1-β₁).*∇V
	∇V_sq_n .= β₂.*∇V_sq_n .+ (1-β₂).*(∇V.^2)
	hat_∇V_n .= ∇V_n./(1 - β₁^(n+1))
	hat_∇V_sq_n .= ∇V_sq_n./(1 - β₂^(n+1))
	X_n .= X_n .- η.*hat_∇V_n./(sqrt.(hat_∇V_sq_n) .+ ϵ)
	return nothing
end

"""
	function random_grad_log_ct!(grad,x₀,model_ct,baseline_contact,λ)

In-place calculation for gradient of the log-contact rate.
"""
function random_grad_log_ct!(grad,x₀,model_ct,baseline_contact,λ)
    k = rand(1:size(model_ct.MCMC_results.chain,1))
    θ = (R=model_ct.MCMC_results.chain[k,1,1],
        E₀ =model_ct.MCMC_results.chain[k,2,1],
        I₀ =model_ct.MCMC_results.chain[k,3,1],
        χ = model_ct.MCMC_results.chain[k,4,1],
        M_PCR = model_ct.MCMC_results.chain[k,5,1],
        α = model_ct.MCMC_results.chain[k,6,1],
        p_test = model_ct.MCMC_results.chain[k,7,1],
        P_eff = model_ct.MCMC_results.chain[k,8,1])
    KenyaSerology.C_t_neg_log_density_grad!(grad,exp.(x₀),θ,model_ct,baseline_contact,λ)
    grad .= grad./exp.(x₀)
    return nothing
end

"""
	function random_grad_log_ct!(grad,x₀,model_ct,baseline_contact,λ,trans_name,p_serowaning)

In-place calculation for gradient of the log-contact rate with sero-waning.
"""
function random_grad_log_ct!(grad,x₀,model_ct,baseline_contact,λ,trans_name,p_serowaning)
    k = rand(1:size(model_ct.MCMC_results.chain,1))
	θ = TransformVariables.transform(trans_name,[model_ct.MCMC_results.chain[1,n,1] for n = 1:size(model_ct.MCMC_results.chain,2)])
    KenyaSerology.C_t_neg_log_density_grad!(grad,exp.(x₀),θ,model_ct,baseline_contact,λ,p_serowaning)
    grad .= grad./exp.(x₀)
    return grad
end

"""
	ADAM_optim(X_opt::Vector{Float64},
					grad_func!::Function,
					total_num_steps::Int64,
					averaging_num_steps::Int64,
					ADAM_params::NamedTuple,
					model::KenyaSerology.CoVAreaModel,
					baseline_contact::Vector{Float64},λ)

Optimize the vector initialized at `X_opt` using ADAM method with `grad_func!` gradient.
"""
function ADAM_optim(X_opt::Vector{Float64},
					grad_func!::Function,
					ADAM_params::NamedTuple,
                    model::KenyaSerology.CoVAreaModel,
                    baseline_contact::Vector{Float64})
	@unpack total_num_steps,averaging_num_steps,λ = ADAM_params
	∇V_n = zeros(length(X_opt))
	∇V_sq_n = zeros(length(X_opt))
    hat_∇V_n = zeros(length(X_opt))
    hat_∇V_sq_n = zeros(length(X_opt))
    grad = zeros(length(X_opt))
	X̂ = zeros(length(X_opt))

	for n = 0:(total_num_steps-1)
        grad_func!(grad,X_opt,model,baseline_contact,λ)
        KenyaSerology.ADAM_step!(X_opt,
                                grad,
                                ∇V_n,∇V_sq_n,
                                hat_∇V_n,hat_∇V_sq_n,
                                n,ADAM_params)
		if n > total_num_steps - averaging_num_steps #
			m = n - (total_num_steps - averaging_num_steps)
			X̂ .+= (X_opt - X̂)/m
		end
	end

	return X̂
end

"""
	function ADAM_optim(X_opt::Vector{Float64},
                        grad_func!::Function,
                        ADAM_params::NamedTuple,
                        model::KenyaSerology.CoVAreaModel,
                        baseline_contact::Vector{Float64},
                        trans_name,
                        p_waning)

Optimize the vector initialized at `X_opt` using ADAM method with `grad_func!` gradient with sero-waning.
"""
function ADAM_optim(X_opt::Vector{Float64},
					grad_func!::Function,
					ADAM_params::NamedTuple,
					model::KenyaSerology.CoVAreaModel,
					baseline_contact::Vector{Float64},
					trans_name,
					p_waning)
	println("Starting ADAM SGD optimisation for effective contact rates.")
	@unpack total_num_steps,averaging_num_steps,λ,σ²_aboveone = ADAM_params

	∇V_n = zeros(length(X_opt))
	∇V_sq_n = zeros(length(X_opt))
	hat_∇V_n = zeros(length(X_opt))
	hat_∇V_sq_n = zeros(length(X_opt))
	grad = zeros(length(X_opt))
	X̂ = zeros(length(X_opt))

	for n = 0:(total_num_steps-1)
		grad_func!(grad,X_opt,model,baseline_contact,λ,trans_name,p_waning)
		KenyaSerology.ADAM_step!(X_opt,
								grad,
								∇V_n,∇V_sq_n,
								hat_∇V_n,hat_∇V_sq_n,
								n,ADAM_params)
		if n > total_num_steps - averaging_num_steps #
			m = n - (total_num_steps - averaging_num_steps)
			X̂ .+= (X_opt - X̂)/m
		end

	end
	final_LL = 0
	for k = 1:size(model.MCMC_results.chain,1)
		θ = TransformVariables.transform(trans_name,[model.MCMC_results.chain[k,n,1] for n = 1:size(model.MCMC_results.chain,2)])
		final_LL += KenyaSerology.loglikelihood_contactratemodel_ct(vcat(baseline_contact,exp.(X̂)),θ,model,λ,p_waning,σ²_aboveone) #NB this is without smoothing penalty or priors just predictive ability
	end
	final_LL = final_LL/size(model.MCMC_results.chain,1)

	return X̂,final_LL
end

"""
	function ADAM_optim(X_opt::Vector{Float64},
                        grad_func!::Function,
                        ADAM_params::NamedTuple,
                        model::KenyaSerology.CoVAreaModel,
                        baseline_contact::Vector{Float64},
                        trans_name,
                        p_waning,
                        plt::Plots.Plot)

Optimize the vector initialized at `X_opt` using ADAM method with `grad_func!` gradient.
Every 10 steps plots onto plt so convergence can be visualized.
"""

function ADAM_optim(X_opt::Vector{Float64},
                        grad_func!::Function,
                        ADAM_params::NamedTuple,
                        model::KenyaSerology.CoVAreaModel,
                        baseline_contact::Vector{Float64},
                        trans_name,
                        p_waning,
                        plt::Plots.Plot)

	println("Starting ADAM SGD optimisation for effective contact rates, with plotting output.")

	@unpack total_num_steps,averaging_num_steps,λ = ADAM_params
        ∇V_n = zeros(length(X_opt))
        ∇V_sq_n = zeros(length(X_opt))
        hat_∇V_n = zeros(length(X_opt))
        hat_∇V_sq_n = zeros(length(X_opt))
        grad = zeros(length(X_opt))
        X̂ = zeros(length(X_opt))

	for n = 0:(total_num_steps-1)
                grad_func!(grad,X_opt,model,baseline_contact,λ,trans_name,p_waning)
                KenyaSerology.ADAM_step!(X_opt,
                        grad,
                        ∇V_n,∇V_sq_n,
                        hat_∇V_n,hat_∇V_sq_n,
                        n,ADAM_params)
                if n > total_num_steps - averaging_num_steps #
                        m = n - (total_num_steps - averaging_num_steps)
                        X̂ .+= (X_opt - X̂)/m
                end
                if n%100 == 0
					LL = 0
					for k = 1:size(model.MCMC_results.chain,1)
                        θ = TransformVariables.transform(trans_name,[model.MCMC_results.chain[1,n,1] for n = 1:size(model.MCMC_results.chain,2)])
                        LL -= KenyaSerology.loglikelihood_contactratemodel_ct(vcat(baseline_contact,exp.(X_opt)),θ,model,λ,p_waning)
					end
					LL = LL/size(model.MCMC_results.chain,1)
                    scatter!(plt,[n],[LL],lab = "",color = :black)
					# display(plot(plt))
                end
	end
        plot!(plt,xlabel = "Steps",ylabel = "Mean Neg. LL",yscale = :log10)
	final_LL = 0
	for k = 1:size(model.MCMC_results.chain,1)
		θ = TransformVariables.transform(trans_name,[model.MCMC_results.chain[k,n,1] for n = 1:size(model.MCMC_results.chain,2)])
		final_LL += KenyaSerology.loglikelihood_contactratemodel_ct(vcat(baseline_contact,exp.(X̂)),θ,model,λ,p_waning)
	end
	final_LL = final_LL/size(model.MCMC_results.chain,1)
	return X̂,final_LL,plt
end

"""
	get_model_and_optimise_ADAM(filename,projected_contact_rate;
									varname = "model",
									enddate::Date,
									fixed_days::Int64,
									ADAM_params::NamedTuple)

Load modelfit from `filename`, and use ADAM SGD optimization to improve C(t).
"""
function get_model_and_optimise_ADAM(filename,projected_contact_rate;
									varname = "model",
									enddate::Date,
									fixed_days::Int64,
									ADAM_params::NamedTuple)

	model = load(filename)[varname]
	# model_ct = deepcopy(model)
	model.prob = make_odeproblemforinference_parameter_ct(projected_contact_rate;
									startdate=Date(2020,2,20),
									enddate=enddate)

	n = (enddate - Date(2020,2,20)).value
	baseline_contact = copy(projected_contact_rate.contactrate[1:fixed_days])
	log_ct = log.(copy(projected_contact_rate.contactrate[(fixed_days+1):end]))
	len_ct = length(projected_contact_rate.contactrate)
	if n > len_ct
		log_ct = vcat(log_ct,fill(log_ct[end],n-len_ct))
	end
	println("Got fit using Google mobility and starting first ADAM grad. descent for $(model.areaname)")
	log_ct_fit = ADAM_optim(copy(log_ct),
				random_grad_log_ct!,
				ADAM_params,
	            model,
	            baseline_contact)

    smoothed_ct = vcat(baseline_contact,exp.(log_ct_fit))
    return smoothed_ct,model
end

"""
	optimise_ct_model_ADAM(model::CoVAreaModel,projected_contact_rate;
							enddate::Date,
							fixed_days::Int64,
							ADAM_params::NamedTuple)

Fit a contact rate using ADAM SGD.
"""
function optimise_ct_model_ADAM(model::CoVAreaModel,projected_contact_rate;
									enddate::Date,
									fixed_days::Int64,
									ADAM_params::NamedTuple)

	model.prob = make_odeproblemforinference_parameter_ct(projected_contact_rate;
									startdate=Date(2020,2,20),
									enddate=enddate)

	n = (enddate - Date(2020,2,20)).value
	baseline_contact = copy(projected_contact_rate.contactrate[1:fixed_days])
	log_ct = log.(copy(projected_contact_rate.contactrate[(fixed_days+1):end]))
	len_ct = length(projected_contact_rate.contactrate)
	if n > len_ct
		log_ct = vcat(log_ct,fill(log_ct[end],n-len_ct))
	end
	println("Refit eff. contact rate using ADAM grad. descent for $(model.areaname)")

	log_ct_fit = ADAM_optim(copy(log_ct),
				random_grad_log_ct!,
				ADAM_params,
	            model,
	            baseline_contact)

	smoothed_ct = vcat(baseline_contact,exp.(log_ct_fit))

    return smoothed_ct
end

"""
	function EM_optimise_ADAM(filename,projected_contact_rate,trans;
							num_steps=3,
							varname="model",
							enddate::Date,
							fixed_days::Int64,
							ADAM_params::NamedTuple)

Perform `num_steps` EM algorithm iterations.
"""
function EM_optimise_ADAM(filename,projected_contact_rate,trans;
						num_steps=3,
						varname="model",
						enddate::Date,
						fixed_days::Int64,
						ADAM_params::NamedTuple)
	#Initial fit of C(t)
    smoothed_ct,model = get_model_and_optimise_ADAM(filename,projected_contact_rate;
													varname = varname,
													enddate=enddate,
													fixed_days=fixed_days,
													ADAM_params=ADAM_params)

	model.prob = make_odeproblemforinference_simple(smoothed_ct,#Method for defining the ODE problem underlying the inference
													startdate = Date(2020,2,21),#Don't change from Feb 21st as start date!
													enddate = enddate)
    model.contactrate_data = copy(smoothed_ct)
	if num_steps > 1
	    KenyaSerology.inferparameters!(model,1000,trans)
	end
	if num_steps == 1
		KenyaSerology.inferparameters!(model,10000,trans)
	end

	#Middle fits of C(t)
    for i = 2:(num_steps-1)
        smoothed_ct = optimise_ct_model_ADAM(model,projected_contact_rate;
											enddate=enddate,
											fixed_days=fixed_days,
											ADAM_params=ADAM_params)

        model.prob = make_odeproblemforinference_simple(smoothed_ct,#Method for defining the ODE problem underlying the inference
                                                        startdate = Date(2020,2,21),#Don't change from Feb 21st as start date!
                                                        enddate = enddate)

        model.contactrate_data = copy(smoothed_ct)
        KenyaSerology.inferparameters!(model,1000,trans)
    end

	#Final fit of C(t)
	if num_steps > 1
		smoothed_ct = optimise_ct_model_ADAM(model,projected_contact_rate;
											enddate=enddate,
											fixed_days=fixed_days,
											ADAM_params=ADAM_params)

		model.prob = make_odeproblemforinference_simple(smoothed_ct,#Method for defining the ODE problem underlying the inference
														startdate = Date(2020,2,21),#Don't change from Feb 21st as start date!
														enddate = enddate)

		model.contactrate_data = copy(smoothed_ct)
		KenyaSerology.inferparameters!(model,10000,trans)
	end
    return model
end


text:
using Triangulate

"""
    constrained_triangulation2D(V::Lar.Points, EV::Lar.Cells) -> Lar.Cells
"""
function constrained_triangulation2D(V::Lar.Points, EV::Lar.Cells)
	triin = Triangulate.TriangulateIO()
	triin.pointlist = V
	triin.segmentlist = hcat(EV...)
	(triout, vorout) = triangulate("pQ", triin)
	trias = Array{Int64,1}[c[:] for c in eachcol(triout.trianglelist)]
	return trias
end

"""
 	LinearAlgebraicRepresentation.triangulate2d(V::Lar.Points, EV::Lar.Cells)

@overwrite Lar.triangulate2d.
"""
function LinearAlgebraicRepresentation.triangulate2d(V::Lar.Points, EV::Lar.Cells)
   	 # data for Constrained Delaunay Triangulation (CDT)
   	 points = convert(Array{Float64,2}, V')
	 # points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
   	 # edges_list = convert(Array{Int64,2}, hcat(EV...)')
   	 # edge_boundary = [true for k=1:size(edges_list,1)] ## dead code !!

	trias = constrained_triangulation2D(V::Lar.Points, EV::Lar.Cells)

 	#Triangle.constrained_triangulation(points,points_map,edges_list)
	innertriangles = Array{Int64,1}[]
	for (u,v,w) in trias
		point = (points[u,:]+points[v,:]+points[w,:])./3
		copEV = Lar.lar2cop(EV)
		inner = Lar.point_in_face(point, points::Lar.Points, copEV::Lar.ChainOp)
		if inner
			push!(innertriangles,[u,v,w])
		end
	end
    return innertriangles
end

ext: jl


text:
## -------------------------------------------------------------------
# MSE per method
let

    p = plot(;xlabel = "experiment", ylabel = "MSE")
    for method in ALL_METHODS
        MSEs = []
        for exp in EXPS

            sum = 0.0
            N = 0
            glc_flx = DAT[method, :Hd, :flx, "GLC", exp]
            for ider in FLX_IDERS
                model_val = DAT[method, :ep, :flx, ider, exp]
                exp_val = DAT[method, :Hd, :flx, ider, exp]
                sum += (model_val/glc_flx - exp_val/glc_flx)^2
                N += 1
            end
            push!(MSEs, sum / N)
        end
        scatter!(p, EXPS, MSEs; color = method_colors[method],
            label = string(method), m = 8, alpha = 0.8,
            legend = :topleft
        )
        plot!(p, EXPS, MSEs; color = method_colors[method],
            label = "", ls = :dash, alpha = 0.8
        )
    end
    mysavefig(p, "MSE_per_method")
end

## -------------------------------------------------------------------
# MSE per ider
let
    p = plot(;xlabel = "experiment", ylabel = "MSE")
    for method in ALL_METHODS
        MSEs = []

        for ider in FLX_IDERS
            sum = 0.0
            N = 0
            for exp in EXPS
                glc_flx = DAT[method, :Hd, :flx, "GLC", exp]
                model_val = DAT[method, :ep, :flx, ider, exp]
                exp_val = DAT[method, :Hd, :flx, ider, exp]
                sum += (model_val/glc_flx - exp_val/glc_flx)^2
                N += 1
            end
            push!(MSEs, sum / N)
        end

        scatter!(p, FLX_IDERS, MSEs; color = method_colors[method],
            label = string(method), m = 8, alpha = 0.8,
            legend = :topleft
        )
        plot!(p, FLX_IDERS, MSEs; color = method_colors[method],
            label = "", ls = :dash, alpha = 0.8
        )
    end
    mysavefig(p, "MSE_per_ider")
end

## -------------------------------------------------------------------
# MSE per beta
let
    method = ME_Z_EXPECTED_G_MOVING

    ps = Plots.Plot[]
    for exp in EXPS
        p = plot(;title = string("exp: ", exp),
            xlabel = "beta", ylabel = "MSE"
        )

        datfile = ME_INDEX[method, :DFILE, exp]
        dat = deserialize(datfile)
        epouts = dat[:epouts]
        betas = epouts |> keys |> collect |> sort
        exp_beta = dat[:exp_beta]
        model = dat[:model]

        MSEs = []
        for beta in betas
            epout = epouts[beta]

            sum = 0.0
            N = 0

            glc_flx = Hd.uval(:GLC, exp)
            for ider in FLX_IDERS

                model_met = iJR.Hd_mets_map[ider]
                model_exch = iJR.exch_met_map[model_met]
                model_exchi = ChU.rxnindex(model, model_exch)

                model_flx = ChU.av(model, epout, model_exchi)
                exp_flx = Hd.uval(ider, exp)

                sum += (model_flx/glc_flx - exp_flx/glc_flx)^2
                N += 1
            end

            push!(MSEs, sum / N)
        end

        scatter!(p, betas, MSEs; color = :black,
            label = "", m = 8, alpha = 0.8
        )
        plot!(p, betas, MSEs; color = :black,
            label = "", ls = :dash, alpha = 0.8
        )
        vline!(p, [exp_beta]; color = :black,
            label = "", ls = :dot, lw = 3, alpha = 0.9
        )
        push!(ps, p)
    end
    mysavefig(ps, "MSE_vs_beta"; method)
end


text:
using Test, ImageFeatures, Images, TestImages, Distributions

@testset "Testing ORB params" begin
    orb_params = ORB(num_keypoints = 1000, threshold = 0.2)
    @test orb_params.num_keypoints == 1000
    @test orb_params.n_fast == 12
    @test orb_params.threshold == 0.2
    @test orb_params.harris_factor == 0.04
    @test orb_params.downsample == 1.3
    @test orb_params.levels == 8
    @test orb_params.sigma == 1.2
end

@testset "Testing with Standard Images - Lighthouse (Rotation 45)" begin
    img = testimage("lighthouse")
    img_array_1 = Gray.(img)
    img_array_2 = _warp(img_array_1, pi / 4)

    #Test Number of Keypoints returned
    orb_params = ORB(num_keypoints = 1000, threshold = 0.1)
    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, orb_params)
    @test length(desc_1) == 1000
    orb_params = ORB(num_keypoints = 500, threshold = 0.1)
    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, orb_params)
    @test length(desc_1) == 500
    orb_params = ORB(num_keypoints = 300, threshold = 0.1)
    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, orb_params)
    @test length(desc_1) == 300

    orb_params = ORB(num_keypoints = 1000)

    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, orb_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_array_2, orb_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.2)
    reverse_keypoints_1 = [_reverserotate(m[1], pi / 4, (256, 384)) for m in matches]
    @test (grade_matches(reverse_keypoints_1, [match[2] for match in matches], 10, (i,j)->(abs(i[1]-j[1]) + abs(i[2]-j[2]))) > 0.95 )
end

@testset "Testing with Standard Images - Lighthouse (Rotation 45, Translation (50, 40))" begin
    img = testimage("lighthouse")
    img_array_1 = Gray.(img)
    img_temp_2 = _warp(img_array_1, pi / 4)
    img_array_2 = _warp(img_temp_2, 50, 40)

    orb_params = ORB(num_keypoints = 1000)

    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, orb_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_array_2, orb_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.2)
    reverse_keypoints_1 = [_reverserotate(m[1], pi / 4, (256, 384)) + CartesianIndex(50, 40) for m in matches]
    @test (grade_matches(reverse_keypoints_1, [match[2] for match in matches], 10, (i,j)->(abs(i[1]-j[1]) + abs(i[2]-j[2]))) > 0.95 )
end

@testset "Testing with Standard Images - Lighthouse (Rotation 75, Translation (50, 40))" begin
    img = testimage("lighthouse")
    img_array_1 = Gray.(img)
    img_temp_2 = _warp(img_array_1, 5 * pi / 6)
    img_array_2 = _warp(img_temp_2, 50, 40)

    orb_params = ORB(num_keypoints = 1000)

    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, orb_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_array_2, orb_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.2)
    reverse_keypoints_1 = [_reverserotate(m[1], 5 * pi / 6, (256, 384)) + CartesianIndex(50, 40) for m in matches]
    @test (grade_matches(reverse_keypoints_1, [match[2] for match in matches], 10, (i,j)->(abs(i[1]-j[1]) + abs(i[2]-j[2]))) > 0.95 )
end

@testset "Testing with Standard Images - Lena (Rotation 45, Translation (10, 20))" begin
    img = testimage("lena_gray_512")
    img_array_1 = Gray.(img)
    img_temp_2 = _warp(img_array_1, pi / 4)
    img_array_2 = _warp(img_temp_2, 10, 20)

    orb_params = ORB(num_keypoints = 1000, threshold = 0.18)

    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, orb_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_array_2, orb_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.2)
    reverse_keypoints_1 = [_reverserotate(m[1], pi / 4, (256, 256)) + CartesianIndex(10, 20) for m in matches]
    @test (grade_matches(reverse_keypoints_1, [match[2] for match in matches], 10, (i,j)->(abs(i[1]-j[1]) + abs(i[2]-j[2]))) > 0.95 )
end


text:
__precompile__()

module Homotopies
    import MultivariatePolynomials
    const MP = MultivariatePolynomials
    import FixedPolynomials
    const FP = FixedPolynomials
    import Base: length, ==

    abstract type AbstractHomotopy{T<:Number} end
    abstract type AbstractPolynomialHomotopy{T<:Number} <: AbstractHomotopy{T} end
    abstract type AbstractHomotopyConfig{T} end
    export AbstractHomotopy, AbstractPolynomialHomotopy, AbstractHomotopyConfig

    include("interface.jl")

    include("config.jl")
    include("diffresult.jl")

    include("homotopies/straightline.jl")
    include("homotopies/geodesic_on_the_sphere.jl")

    include("constructors.jl")

    include("polynomial.jl")

    include("totaldegree.jl")
    include("randomhomotopy.jl")

    include("condition.jl")

    include("bigfloat_utilities.jl")
end


text:
"""
    conformgt_by_id(reffile, tgtfile, outfile, chrom, posrange, checkfreq)

Match the VCF records in `tgtfile` to those in `reffile` according to ID.
The function will:
0. Find corresponding VCF records in the target and reference files
0. Exclude target VCF records whose ID cannot be matched to any reference VCF record
0. Exclude target VCF records whose test of equal allele frequency is rejected at significance level `checkfreq`
0. Adjust target VCF records so that chromosome strand and allele order match the VCF reference file
0. The matched VCF records are written into files `outfile.tgt.vcf.gz` and
`outfile.ref.vcf.gz`, both with only "GT" data

# Input
- `reffile`: VCF file with reference genotype (GT) data
- `tgtfile`: VCF file with target genotype (GT) data
- `outfile`: the prefix for output filenames
- `chrom`: chromosome name, must be identical in target and reference files
- `posrange`: position range in the reference file
- `checkfreq`: significance level for testing equal alelle frequencies between mached target and reference records. If the test pvalue is  `≤ checkfreq`, the records are not output. Setting `checkfreq=0` or `checkfreq=false` (default) implies not checking allele frequencies. Setting `checkfreq=1` effectively rejects all tests and no matched records are output

# Output
- `lines`: number of matched VCF records
"""
function conformgt_by_id(
    reffile::AbstractString,
    tgtfile::AbstractString,
    outfile::AbstractString,
    chrom::AbstractString,
    posrange::AbstractRange,
    checkfreq::Number = false
    )
    # open reference and target VCF files
    reader_ref = VCF.Reader(openvcf(reffile, "r"))
    reader_tgt = VCF.Reader(openvcf(tgtfile, "r"))
    records_ref, records_tgt = nrecords(reffile), nrecords(tgtfile)
    # create output files
    writer_ref = VCF.Writer(openvcf(join([outfile, ".ref.vcf.gz"]), "w"),
        VCF.header(reader_ref))
    writer_tgt = VCF.Writer(openvcf(join([outfile, ".tgt.vcf.gz"]), "w"),
        VCF.header(reader_tgt))
    # collect IDs in the reference panel
    @info("Scan IDs in reference panel")
    pbar = ProgressMeter.Progress(records_ref, 1)
    id_list_ref = Vector{String}[]
    record_counter = 0
    for record_ref in reader_ref
        # update progress bar
        record_counter += 1
        ProgressMeter.update!(pbar, record_counter)
        # chromosome not matched
        VCF.chrom(record_ref) == chrom || continue
        # if past search range, break
        pos_ref = VCF.pos(record_ref)
        if pos_ref > last(posrange)
            ProgressMeter.update!(pbar, records_ref)
            break
        end
        # if not in search range, skip this record
        pos_ref < first(posrange) && continue
        # if no "GT" field, skip this record
        VCF.findgenokey(record_ref, "GT") == nothing && continue
        # if no ID, skip this record
        VCF.hasid(record_ref) || continue
        # add ID to list
        id_ref = VCF.id(record_ref)
        if length(id_list_ref) > 0 && id_ref == id_list_ref[end]
            print("Duplicate ID in reference panel: $(id_ref); ")
            print("will use the first occurence\n")
            continue
        end
        push!(id_list_ref, id_ref)
    end
    close(reader_ref)
    # match target IDs to reference IDs
    @info("Match target IDs to reference IDs")
    # re-set reference reader
    reader_ref = VCF.Reader(openvcf(reffile, "r"))
    record_ref = VCF.Record()
    iter_state_ref = iterate(reader_ref)
    # loop over target records
    pbar = ProgressMeter.Progress(records_tgt, 1)
    lines = record_counter = 0
    for record_tgt in reader_tgt
        # update progress bar
        record_counter += 1
        ProgressMeter.update!(pbar, record_counter)
        # if no "GT" field, skip this record
        VCF.findgenokey(record_tgt, "GT") == nothing && continue
        # if no ID, skip
        VCF.hasid(record_tgt) || continue
        # if not in reference panel, skip
        # TODO can be optimized here, shrink id_list_ref
        id_tgt = VCF.id(record_tgt)
        id_tgt ∈ id_list_ref || continue
        # search for matching reference record
        while iter_state_ref !== nothing
            record_ref, state_ref = iter_state_ref
            VCF.hasid(record_ref) && VCF.id(record_ref) == id_tgt && break
            iter_state_ref = iterate(reader_ref, state_ref)
        end
        # check REF/ALT label match
        reflabel_tgt = VCF.ref(record_tgt)
        altlabel_tgt = VCF.alt(record_tgt)
        reflabel_ref = VCF.ref(record_ref)
        reflabel_tgt == reflabel_ref || reflabel_ref ∈ altlabel_tgt || continue
        record_gt_ref, record_gt_tgt = match_gt_allele(record_ref, record_tgt)
        # allele frequency test if requested
        if checkfreq > 0
            _, _, _, n0_ref, n1_ref, = gtstats(record_gt_ref, nothing)
            _, _, _, n0_tgt, n1_tgt, = gtstats(record_gt_tgt, nothing)
            pval = binomial_proportion_test(n0_ref, n1_ref, n0_tgt, n1_tgt)
            pval ≤ checkfreq && continue
        end
        # write to matched target and reference file
        lines += 1
        VCF.write(writer_ref, record_gt_ref)
        VCF.write(writer_tgt, record_gt_tgt)
    end
    close(reader_ref); close(reader_tgt)
    flush(writer_ref); flush(writer_tgt)
    close(VCF.BioCore.IO.stream(writer_ref))
    close(VCF.BioCore.IO.stream(writer_tgt))
    # return
    @info("$(lines) records are matched")
    return lines
end

"""
    conformgt_by_pos(reffile, tgtfile, outfile, chrom, posrange, checkfreq)

Match the VCF GT records in `tgtfile` to those in `reffile` according to
chromosome position. The function will:
0. Find corresponding VCF records in the target and reference files
0. Exclude target VCF records whose position cannot be matched to any reference VCF record
0. Exclude target VCF records whose test of equal allele frequency is rejected at significance level `checkfreq`
0. Adjust target VCF records so that chromosome strand and allele order match the VCF reference file
0. The matched VCF records are written into files `outfile.tgt.vcf.gz` and
`outfile.ref.vcf.gz`, both with only "GT" data

# Input
- `reffile`: VCF file with reference genotype (GT) data
- `tgtfile`: VCF file with target genotype (GT) data
- `outfile`: the prefix for output filenames
- `chrom`: chromosome name, must be identical in target and reference files
- `posrange`: position range in the reference file
- `checkfreq`: significance level for testing equal alelle frequencies between mached target and reference records. If the test pvalue is  `≤ checkfreq`, the records are not output. Setting `checkfreq=0` or `checkfreq=false` (default) implies not checking allele frequencies. Setting `checkfreq=1` effectively rejects all tests and no matched records are output

# Output
- `lines`: number of matched VCF records
"""
function conformgt_by_pos(
    reffile::AbstractString,
    tgtfile::AbstractString,
    outfile::AbstractString,
    chrom::AbstractString,
    posrange::AbstractRange,
    checkfreq::Number = false
    )
    # open reference and target VCF files
    reader_ref = VCF.Reader(openvcf(reffile, "r"))
    reader_tgt = VCF.Reader(openvcf(tgtfile, "r"))
    records_ref, records_tgt = nrecords(reffile), nrecords(tgtfile)
    # create output files (.gz currently throws error)
    writer_ref = VCF.Writer(openvcf(join([outfile, ".ref.vcf.gz"]), "w"),
        VCF.header(reader_ref))
    writer_tgt = VCF.Writer(openvcf(join([outfile, ".tgt.vcf.gz"]), "w"),
        VCF.header(reader_tgt))
    # initialize reference reader
    record_ref, iter_state_ref, pos_ref = VCF.Record(), iterate(reader_ref), 0
    # loop over target records
    @info("Match target POS to reference POS")
    pbar = ProgressMeter.Progress(records_tgt, 1)
    lines = record_counter = 0
    for record_tgt in reader_tgt
        # update progress bar
        record_counter += 1
        ProgressMeter.update!(pbar, record_counter)
        # if no "GT" field, skip this record
        VCF.findgenokey(record_tgt, "GT") == nothing && continue
        # chromosome not matched, skip this record
        VCF.chrom(record_tgt) == chrom || continue
        # if past search range, break
        pos_tgt = VCF.pos(record_tgt)
        if pos_tgt > last(posrange)
            ProgressMeter.update!(pbar, records_tgt)
            break
        end
        # if before search range, skip this record
        pos_tgt < first(posrange) && continue
        # if before the current ref position, skip this record
        pos_tgt < pos_ref && continue
        if pos_tgt == pos_ref
            matched = true
        elseif pos_tgt > pos_ref
            matched = false
            # search for matching reference record
            while iter_state_ref !== nothing
                record_ref, state_ref = iter_state_ref
                # if no "GT" field, skip this record
                VCF.findgenokey(record_tgt, "GT") == nothing && continue
                # if chromosome not matched, skip this record
                VCF.chrom(record_ref) == chrom || continue
                pos_ref = VCF.pos(record_ref)
                if pos_ref < pos_tgt
                    iter_state_ref = iterate(reader_ref, state_ref)
                elseif pos_ref == pos_tgt
                    matched = true; break
                elseif pos_ref > pos_tgt
                    break
                end
            end
        end
        # if no matching reference record, skip this record
        matched || continue
        # check REF/ALT label match
        reflabel_tgt = VCF.ref(record_tgt)
        altlabel_tgt = VCF.alt(record_tgt)
        reflabel_ref = VCF.ref(record_ref)
        reflabel_tgt == reflabel_ref || reflabel_ref ∈ altlabel_tgt || continue
        record_gt_ref, record_gt_tgt = match_gt_allele(record_ref, record_tgt)
        # allele frequency test if requested
        if checkfreq > 0
            _, _, _, n0_ref, n1_ref, = gtstats(record_gt_ref, nothing)
            _, _, _, n0_tgt, n1_tgt, = gtstats(record_gt_tgt, nothing)
            pval = binomial_proportion_test(n0_ref, n1_ref, n0_tgt, n1_tgt)
            pval ≤ checkfreq && continue
        end
        # write to matched target and reference file
        lines += 1
        VCF.write(writer_ref, record_gt_ref)
        VCF.write(writer_tgt, record_gt_tgt)
    end
    close(reader_ref); close(reader_tgt)
    flush(writer_ref); flush(writer_tgt)
    close(VCF.BioCore.IO.stream(writer_ref))
    close(VCF.BioCore.IO.stream(writer_tgt))
    # return
    @info("$(lines) records are matched")
    return lines
end

"""
    match_gt_allele(record1, record2)

Match the REF allele label of record 2 to that of record 1. Return the possibly
REF/ALT flipped record 2 with only "GT" genotype data.
"""
function match_gt_allele(record1::VCF.Record, record2::VCF.Record)
    VCF.findgenokey(record1, "GT") == nothing && error("record1 has no GT field")
    VCF.findgenokey(record2, "GT") == nothing && error("record2 has no GT field")
    ref1 = VCF.ref(record1)
    ref2 = VCF.ref(record2)
    alt2 = VCF.alt(record2)
    record1_out = filter_genotype(record1, ["GT"])
    if ref1 == ref2
        # matching REF labels
        record2_out = filter_genotype(record2, ["GT"])
    elseif ref1 ∈ alt2
        # switch REF/ALT label of record2
        record2_out = flip_gt_allele(record2, ref1)
    else
        throw(ArgumentError("$(ref1) not found in record 2"))
    end
    record1_out, record2_out
end

"""
    filter_genotype(record, [genokey=["GT"]])

Filter a VCF record according to `genokey` and output a VCF record with
genotype formats only in `genokey`.
"""
function filter_genotype(
    record::VCF.Record,
    genokey::Vector{T} = ["GT"]
    ) where T <: AbstractString
    format_out = genokey ∩ VCF.format(record)
    gt_out = [Dict(gk => VCF.genotype(record, i, gk) for gk in format_out)
        for i in 1:length(record.genotype)]

    return VCF.Record(record; genotype = gt_out)
end

"""
    filter_genotype(vcffile, outfile, [genokey=["GT"]])

Filter a VCF file according to `genokey` (default `["GT"]`). Output a VCF
file with genotype formats only in `genokey`. Record (markers) with no fields
in `genokey` are skipped.
"""
function filter_genotype(
    vcffile::T,
    outfile::T = "filtered.vcf.gz",
    genokey::Vector{T} = ["GT"]
    ) where T <: AbstractString
    reader = VCF.Reader(openvcf(vcffile, "r"))
    writer = VCF.Writer(openvcf(outfile, "w"), VCF.header(reader))
    for record in reader
        VCF.write(writer, filter_genotype(record, genokey))
    end
    flush(writer); close(reader); close(writer)
    nothing
end

"""
    flip_gt_allele(record[, altlabel])

Flip the REF and ALT allele lable in a VCF record and output a VCF record with
only "GT" genotype format. By default the first label in ALT field is used as
REF label in the flipped record.
"""
function flip_gt_allele(record::VCF.Record, altlabel::String = VCF.alt(record)[1])
    altlabel ∈ VCF.alt(record) || throw(ArgumentError("ALT label not found"))
    VCF.findgenokey(record, "GT") == nothing && throw(ArgumentError("GT format not found"))
    gt_out = [Dict("GT" => VCF.genotype(record, i, "GT"))
        for i in 1:length(record.genotype)]
    record_out = VCF.Record(record; ref = VCF.alt(record)[1],
        alt = [VCF.ref(record)], genotype = gt_out)
    flip_01!(record_out.data, record_out.genotype[1][1][1]:record_out.genotype[end][end][end])
    return record_out
end

"""
    flip_01!(s[, r])
Flip the digits 0 and 1 in a UInt8 vector `s` in range `r`.
"""
function flip_01!(s::Vector{UInt8}, r::AbstractRange = 1:length(s))
    for i in r
        if s[i] == 0x30
            s[i] = 0x31
        elseif s[i] == 0x31
            s[i] = 0x30
        end
    end
    s
end

"""
    binomial_proportion_test(x1, x2, y1, y2)

Perform hypothesis test for the equality of two binomial proportions. The
contingency table is structured as
        | Group 1 | Group 2
Count 1 |    x1   |   x2
Count 2 |    y1   |   y2
It tests whether `x1 / (x1 + y1)` is equal to `x2 / (x2 + y2)`. When any of the
counts is less than 5, it resorts to the Fisher's exact test.
"""
function binomial_proportion_test(x1::T, y1::T, x2::T, y2::T) where T <: Integer
    if min(x1, y1, x2, y2) < 5
        # Fisher exact test
        pval = pvalue(FisherExactTest(x1, x2, y1, y2))
    else
        # large sample test
        n1, n2 = x1 + y1, x2 + y2
        p̂1, p̂2, p̂ = x1 / n1, x2 / n2, (x1 + x2) / (n1 + n2)
        ẑ = (p̂1 - p̂2) / √(p̂ * (1 - p̂) * (1 / n1 + 1 / n2))
        pval = ẑ > 0 ? 2ccdf(Normal(), ẑ) : 2cdf(Normal(), ẑ)
    end
    pval
end


text:
using Documenter, Colors

makedocs(
    modules = [Colors],
    format = :html,
    sitename = "Colors",
    pages    = Any[
        "Introduction"             => "index.md",
        "Colorspaces"              => "colorspaces.md",
        "Colorscales"              => "colorscales.md",
        "Colormaps"                => "colormaps.md",
        "Named colors"             => "namedcolors.md",
        "References"               => "references.md",
        "Migrating from Color.jl"  => "migratingfromcolor.md",
        "Index"                    => "functionindex.md",
        ]
    )

deploydocs(
    repo = "github.com/JuliaGraphics/Colors.jl.git",
    target = "build",
    julia  = "0.7",
    osname = "linux",
    deps = nothing,
    make = nothing
)


text:
import MLBase
import StatsBase

"""
"""
function roccurve(
        ytrue::AbstractVector{<:Integer},
        yscore::AbstractVector{<:Real},
        )
    allrocnums, allthresholds = getallrocnums(
        ytrue,
        yscore,
        )
    result = roccurve(
        allrocnums,
        allthresholds,
        )
    return result
end

"""
"""
function roccurve(
        allrocnums::AbstractVector{<:MLBase.ROCNums},
        allthresholds::AbstractVector{<:Real},
        )
    allfpr = [false_positive_rate(x) for x in allrocnums]
    alltpr = [true_positive_rate(x) for x in allrocnums]
    #
    permutation = sortperm(allthresholds; rev = false)
    allfpr = allfpr[permutation]
    alltpr = alltpr[permutation]
    allthresholds = allthresholds[permutation]
    return allfpr, alltpr, allthresholds
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module SimplePilotResponseImpl

export
    initialize,
    update,

    updatePilotResponse,

    SimplePilotResponse,
    SimplePRResolutionAdvisory,
    SimplePRCommand

using AbstractPilotResponseImpl
using AbstractPilotResponseInterfaces
using CommonInterfaces

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractPilotResponseInterfaces.updatePilotResponse

type SimplePRResolutionAdvisory
  h_d::Float64
end

type SimplePRCommand

  t::Float64
  v_d::Float64
  h_d::Float64
  psi_d::Float64

  logProb::Float64
end

type SimplePilotResponse <: AbstractPilotResponse

  b_CAS_activated::Bool
  RA::SimplePRResolutionAdvisory

  function SimplePilotResponse()

    obj = new()

    obj.b_CAS_activated = false

    return obj
  end
end

function updatePilotResponse(pr::SimplePilotResponse, command::SimplePRCommand, RA::Union{SimplePRResolutionAdvisory,Void})

  t, v_d, h_d, psi_d = command.t, command.v_d, command.h_d, command.psi_d

  if RA != nothing && pr.b_CAS_activated == false
    pr.RA = RA
    pr.b_CAS_activated = true
  end

  if pr.b_CAS_activated
    h_d = pr.RA.h_d
  end

  return SimplePRCommand(t, v_d, h_d, psi_d, 0.0)
end

update(pr::SimplePilotResponse, command, RA) = update(pr,convert(SimplePRCommand, command), convert(SimplePRResolutionAdvisory,RA))

update(pr::SimplePilotResponse, command::SimplePRCommand, RA::Union{SimplePRResolutionAdvisory,Void}) = updatePilotResponse(pr, command, RA)

function initialize(pr::SimplePilotResponse)

  pr.b_CAS_activated = false
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
#########################################
# Mala for comparison with adaptation   #
#########################################
using Random
include("../../src/faberschauder.jl")
include("../../src/asvar.jl")
include("../../src/pdmp.jl")

function get_info(ϕ, ∇ϕ, ξ, L, α, T)
    U = 0.5*dot(ξ, ξ) #gaussian part
    ∇U = deepcopy(ξ) #gaussian part
    dt = T/(2<<L)
    t = dt*0.5
    tt = 0:T/(2<<L):T
    N = Int(T/dt)
    for i in 1:N
        Xₜ = dotψ(ξ, t, L, T)
        U += ϕ(Xₜ, ξ, t, L, α, T)*dt
        for i in 2:(length(ξ)-1)
            l = lvl(i, L)
            k = (i - 1) ÷ (2 << l)
            δ = T / (1 << (L - l))
            if δ*k < t < δ*(k+1)
                ∇U[i] += ∇ϕ(Xₜ, l, t, ξ,  L, α, T)*dt
            end
        end
        t += dt
    end
    ∇U[1] = ∇U[end] = 0.0
    return ∇U, U
end

function mala_sampler(ϕ, ∇ϕ, ξ₀::Vector{Float64}, niter::Int64, args...)
        ξtrace = Array{Float64}(undef, length(ξ₀), niter)
        ξtrace[:,1] = ξ₀
        ∇U₀, U₀ = get_info(ϕ, ∇ϕ, ξ₀, args...)
        ξ₁ = deepcopy(ξ₀)
        τ = 0.015
        count = 0
        for i in 2:niter
                # initial and final point fixed
                for j in 2:length(ξ₀)-1
                    ξ₁[j]  = ξ₀[j] - τ*∇U₀[j] + sqrt(2*τ)*randn()
                end
                ∇U₁, U₁ = get_info(ϕ, ∇ϕ, ξ₁, args...)
                acc_rej =  U₀ - U₁ + (norm(ξ₁ - ξ₀ + τ*∇U₀)^2 - norm(ξ₀ - ξ₁ + τ*∇U₁)^2)/(4τ)
                # @show acc_rej
                # @assert acc_rej == 0.0
                if acc_rej > log(rand())
                        ∇U₀, U₀, ξ₀ = deepcopy(∇U₁), deepcopy(U₁), deepcopy(ξ₁)
                        count += 1
                end
                if i % 100 == 0
                    if count/100 <= 0.6
                        τ = max(0.0001, τ - 0.0001)
                    else
                        τ = min(1.0, τ + 0.0001)
                    end
                    #println("Adaptive step: ar: $(count/100), new tau: $τ")
                    count = 0
                end
                ξtrace[:,i] = ξ₀
        end
        ξtrace
end

# sin application
function ϕ(Xₜ, ξ, t, L, α, T)
    0.5*α*(α*sin(Xₜ)^2 + cos(Xₜ))
end

function ∇ϕ(Xₜ, l, t, ξ,  L, α, T)
    Λ(t, L - l, T)*0.5*(α^2*sin(2.0*Xₜ) - α*sin(Xₜ))
end

function run_mala(df, niter, α)
    L = 6
    T = 50.0
    n = (2 << L) + 1
    ξ₀ =  randn(n)
    u, v = 0.0, 0.0  # initial and fianl point
    ξ₀[1] = u / sqrt(T)
    ξ₀[end] = v / sqrt(T)
    dt = 1/(2 << L)
    mala_time = @elapsed (out = mala_sampler(ϕ, ∇ϕ, ξ₀, niter, L, α, T))
    burn1 = 1000
    xT_2 = sqrt(T)/2*(out[Int((length(ξ₀)+1)/2), burn1+1:niter])
    nn = length(xT_2)
    ave = cumsum(xT_2)./(Float64.(1:nn))
    for batches in [50]
        ess = ESS(out[:,burn1:end], n_batches = batches)
        push!(df, Dict(:sampler => "MALA", :alpha => α, :T => niter, :nbatches => batches,
            :stat => "ess_mean", :y => sum(ess[2:end-1])/(length(ess)-2), :runtime => mala_time))
        push!(df, Dict(:sampler => "MALA", :alpha => α, :T => niter, :nbatches => batches,
            :stat => "ess_median", :y => median(ess[2:end-1]), :runtime => mala_time))
        push!(df, Dict(:sampler => "MALA", :alpha => α, :T => niter, :nbatches => batches,
            :stat => "ess_min", :y => minimum(ess[2:end-1]), :runtime => mala_time))
        push!(df, Dict(:sampler => "MALA", :alpha => α, :T => niter, :nbatches => batches,
            :stat => "ess_x_T2", :y => ess[Int((length(ξ₀)+1)/2)], :runtime => mala_time))
    end
    # S = T*(0:n)/(n+1)
    # p1 = lines(S, [dotψ(xx[:,end], s, L, T) for s in S], linewidth=0.3)
    # for i in 1:100:size(xx)[2]
    #     lines!(p1, S, [dotψ(xx[:, i], s, L, T) for s in S], linewidth=0.1, alpha = 0.1)
    # end
    # display(p1)
    return df, ave
end

function data_collection_mala(df)
    Random.seed!(0)
    niter = 250000
    mala = []
    for α in [0.1, 0.3, 0.5]
        df, x = run_mala(df, niter, α)
        push!(mala, x)
    end
    return df, mala
end
df, mala = data_collection_mala(df)
# With Automatic Integration (TOO SLOW)
using CSV, JLD2
CSV.write("./scripts/zz_diff_bridges/compare/benchamrk_complete.csv", df)
@save  "./scripts/zz_diff_bridges/compare/covenrgence_mala.jld2" mala
error("Stop here")

cumsum(randn(5,2))

using QuadGK

# sin application
function ϕ(ξ, t, L, α, T)
    Xₜ = dotψ(ξ, t, L, T)
    0.5*α*(α*sin(Xₜ)^2 + cos(Xₜ))
end
function ∇ϕ(l, t, ξ,  L, α, T)
    Xₜ = dotψ(ξ, t, L, T)
    #println(Λ(t, L - l, T)*0.5*(α^2*sin(2.0*Xₜ) - α*sin(Xₜ)))
    0.5*(α^2*sin(2.0*Xₜ) - α*sin(Xₜ))*Λ(t, L - l, T)
end

function get_info(ϕ, ∇ϕ, ξ, L, α, T)
    U = 0.5*dot(ξ, ξ) #gaussian part
    ∇U = deepcopy(ξ) #gaussian part
    U += quadgk(t ->  ϕ(ξ, t, L, α, T), 0.0, T, rtol = 1e-1)[1]
    for i in 2:(length(ξ)-1)
        l = lvl(i, L)
        k = (i - 1) ÷ (2 << l)
        δ = T / (1 << (L - l))
        ∇U[i] += quadgk(t -> ∇ϕ(l, t, ξ, L, α, T), δ*k, δ*(k+1), 1e-1)[1]
    end
    ∇U[1] = ∇U[end] = 0.0
    return ∇U, U
end
runall()

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
using Combinatorics

VecTup = Union{AbstractVector, Tuple}

"""
This implements
```
   (N!)^{-1} ∑_{σ1, σ2} < Y_l1^m1 ∘ σ1, Y_l2^m2 ∘ σ2 >
```
"""
symdotYlm(l1, m1, l2, m2) =
   l1 != l2 ? 0 :  sum( (l1[σ] == l2) * (m1[σ] == m2)
                        for σ in permutations(1:length(l1)) )

function randlm(N, maxl=4)
   l = sort( rand(0:maxl, N) )
   m1 = [ rand(-l[a]:l[a]) for a = 1:N ]
   m2 = [ rand(-l[a]:l[a]) for a = 1:N ]
   return tuple(l...), tuple(m1...), tuple(m2...)
end

l, m1, m2 = randlm(6)

using BenchmarkTools
@btime symdotYlm($l, $m1, $l, $m2)

symdotYlm(l, m1, l, m2)

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 9 Samuel Powell

# examples/gain.jl: example to configure manual gain
using Spinnaker

camlist = CameraList()
cam = camlist[0]
triggersource!(cam, "Software")
triggermode!(cam, "On")
exposure!(cam, 1e5)
acquisitionmode!(cam, "Continuous")
start!(cam)
for gainval in 0:2:48
    gain!(cam, gainval)
    @info "Gain set to $(gainval)dB"
    trigger!(cam)
    saveimage(cam, joinpath(@__DIR__, "gain_$(gainval)dB.png"), spinImageFileFormat(6))
    @info "Image saved"
end
stop!(cam)


text:
facts("Testing simple simulation with fixed agents and no challenge") do
    srand(1234)
    results = TcrCore.simSimple(1000, setupAgentsWithFixedAccuracy(100, 50), [TcrCore.benchmarkRegistryMean])
    print(results)
    @fact results[1] --> 72.20784862274684
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
"""
Hardware threads of device
"""
threads(device) = 0

"""
Blocks that group together hardware threads
"""
blocks(device) = 1
"""
Global memory, e.g. VRAM or RAM of device
"""
global_memory(device) = 0

"""
Free global memory. Isn't supported for AMD cards right now, in which case it returns NaN,
so don't rely on the output of this function.
"""
free_global_memory(device) = NaN

"""
Block local memory
"""
local_memory(device) = 0

"""
Hardware name of a device
"""
name(device) = "Undefined"

"""
Summarizes all features of a device and prints it to `io`
"""
function device_summary(io::IO, device)
    println(io, "Device: ", name(device))
    for (n, f) in (:threads => threads, :blocks => blocks)
        @printf(io, "%19s: %s\n", string(n), string(f(device)))
    end
    for (n, f) in (:global_memory => global_memory, :free_global_memory => free_global_memory, :local_memory => local_memory)
        @printf(io, "%19s: %f mb\n", string(n), f(device) / 10^6)
    end
    return
end

################################
# Device selection functions for e.g. devices(filterfuncs)
"""
Returns true if `device` is a gpu
"""
is_gpu(device) = false

"""
Returns true if `device` is a cpu
"""
is_cpu(device) = false

"""
Checks a device for a certain attribute and returns true if it has at least `value`.
Can be used with e.g. `threads`, `blocks`, `global_memory`, `local_memory`
"""
has_atleast(device, attribute, value) = attribute(ctx_or_device) >= value

#################################
# Context filter functions
# Works for context objects as well but is overloaded in the backends
is_opencl(ctx::Symbol) = ctx == :opencl
is_cudanative(ctx::Symbol) =  ctx == :cudanative
is_julia(ctx::Symbol) =  ctx == :threaded
is_opengl(ctx::Symbol) =  ctx == :opengl

is_opencl(ctx) = false
is_cudanative(ctx) = false
is_julia(ctx) = false
is_opengl(ctx) = false

supports_double(ctx) = false

#=
Functions to select contexts
=#

"""
Creates a new context from `device` without caching the resulting context.
"""
function new_context(device)
    error("Device $device not supported")
end

"""
Destroys context, freeing all it's resources.
"""
function destroy!(context)
    error("Device $context not supported")
end

"""
Resets a context freeing all resources and creating a new context.
"""
function reset!(context)
    error("Context $context not supported")
end

"""
Creates a temporary context for `device` and executes `f(context)` while this context is active.
Context gets destroyed afterwards. Note, that creating a temporary context is expensive.
"""
function on_device(f, device = current_device())
    ctx = new_context(device)
    try
        f(ctx)
    finally
        destroy!(ctx)
    end
    return
end

const filterfuncs = """
Device can be filtered by passing `filter_funcs`, e.g. :
`is_gpu`, `is_cpu`, `(dev)-> has_atleast(dev, threads, 512)`
"""

"""
Returns all devices for the current backend.
$filterfuncs
"""
function available_devices(filter_funcs...)
    result = []
    for device in current_backend().devices()
        if all(f-> f(device), filter_funcs)
            push!(result, device)
        end
    end
    result
end

"""
Returns all devices from `backends = active_backends()`.
$filterfuncs
"""
function all_devices(filter_funcs...; backends = active_backends())
    result = []
    for Module in backends
        for device in Module.devices()
            if all(f-> f(device), filter_funcs)
                push!(result, device)
            end
        end
    end
    result
end

"""
Iterates through all available devices and calls `f(context)` after initializing the standard context for that device.
"""
function forall_devices(func, filterfuncs...)
    for device in all_devices(filterfuncs...)
        ctx = init(device)
        func(ctx)
    end
end


text:
import SpecialFunctions: loggamma

@node MatrixDirichlet Stochastic [ out, a ]

@average_energy MatrixDirichlet (q_out::MatrixDirichlet, q_a::PointMass) = begin
    H = mapreduce(+, zip(eachcol(mean(q_a)), eachcol(logmean(q_out)))) do (q_a_column, logmean_q_out_column)
        return -loggamma(sum(q_a_column)) + sum(loggamma.(q_a_column)) - sum((q_a_column .- 1.0) .* logmean_q_out_column)
    end
    return H
end


text:
@everywhere using DistributedArrays

println(workers())

n = length(workers())
println("#workers: $n")

# pixels = dzeros((3, 4), workers(), [1,n])
# pixels = dzeros(Int32, (3, 4, n))

chunks = (1, 1, n)
pixels = DArray(I -> zeros(Int32, map(length, I)), (3, 4, n), workers(), chunks)

# mask = dfill(false, (3, 4, n), workers())
mask = DArray(I -> fill(false, map(length, I)), (3, 4, n), workers(), chunks)

@everywhere function fill_array(d::DArray)
    pid = myid()

    la = localpart(d)
    la .= pid
    la[1,1] = pid - 1

    println("size(la) = $(size(la))")
end

@sync for w in workers()
    @spawnat w fill_array(pixels)
end

println("pixels:", pixels)
println("mask:", mask)

println("indices:", [@fetchfrom p localindices(pixels) for p in workers()])


text:
using DataStructures
using DelimitedFiles

# returns solution of part 1
function x()
	q = []
	f = open("9.txt")
	for i=1:25
		n = parse(Int, readline(f))
		push!(q, n)
	end
	for i=1:1000-25
		n = parse(Int, readline(f))
		l = false
		for a=1:25
			for b=a:25
				if q[a] + q[b] == n
					l = true
				end
			end
		end
		if !l
			return n
		end
		for j=1:24
			q[j] = q[j+1]
		end
		q[25] = n
	end
end

# returns solution of part 2 (slow solution: 9.275257 seconds)
function y()
	t = readdlm("9.txt", Int64)
	n = x()
	p = (0,0)
	for i=1:length(t), j=i+2:length(t)
		s = 0
		for k=i:j
			s += t[k]
		end
		if s == n
			p = (i,j)
			break
		end
	end
	minimum(t[p[1]:p[2]]) + maximum(t[p[1]:p[2]])
end


text:
using CUDA.CUSPARSE
using SparseArrays
using Test
typeSet = [Float32, Float64, ComplexF32, ComplexF64]
n = 10
@testset "similar" begin
    @testset "similar(A::$h{$elty},Tv::$newElty)" for elty in typeSet,
                                     newElty in typeSet,
                                     (h,dimPtr,dimVal) in ((CuSparseMatrixCSR,:rowPtr,:colVal), (CuSparseMatrixCSC,:colPtr,:rowVal))

        A = sprand(elty, n, n, rand())
        dA = h(A)

        C_simple = similar(dA)
        C_dims = similar(dA,(n,n+1))
        C_eltype = similar(dA,newElty)
        C_full = similar(dA,newElty,(n,n+1))
        @test typeof(C_simple) == typeof(dA)
        @test typeof(C_dims) == typeof(dA)
        @test (typeof(C_eltype) == typeof(dA) && elty==newElty) || ((typeof(C_eltype) != typeof(dA) && elty!=newElty))
        @test (typeof(C_full) == typeof(dA) && elty==newElty) || ((typeof(C_full) != typeof(dA) && elty!=newElty))

        properties = Set(propertynames(dA));

        conserved_simple = Set([dimPtr,dimVal,:dims,:nnz])
        structure_conserved_simple = setdiff(properties,conserved_simple);
        @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname simple conserved" for propertyname in conserved_simple
            @test getproperty(C_simple,propertyname) == getproperty(dA,propertyname)
        end

        @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname simple structure conserved" for propertyname in structure_conserved_simple
            @test length(getproperty(C_simple,propertyname)) == length(getproperty(dA,propertyname))
            @test eltype(getproperty(C_simple,propertyname)) == eltype(getproperty(dA,propertyname))
        end

        conserved_dims = Set([:nnz])
        if h==CuSparseMatrixCSR # Making the array one column longer increases colPtr length but not rowPtr length
            structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims])))
        else #CSC
            structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims,dimPtr])))
            @test length(getproperty(C_dims,dimPtr)) == length(getproperty(dA,dimPtr)) + 1
            @test eltype(getproperty(C_dims,dimPtr)) == eltype(getproperty(dA,dimPtr))
        end
        @test getproperty(C_dims,:dims) == (n,n+1)
        @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname dims conserved" for propertyname in conserved_dims
            @test getproperty(C_dims,propertyname) == getproperty(dA,propertyname)
        end

        @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname dims structure conserved" for propertyname in structure_conserved_dims
            @test length(getproperty(C_dims,propertyname)) == length(getproperty(dA,propertyname))
            @test eltype(getproperty(C_dims,propertyname)) == eltype(getproperty(dA,propertyname))
        end

        conserved_eltype = Set([:nnz,:dims,dimPtr,dimVal])
        structure_conserved_eltype = setdiff(properties,union(conserved_eltype,[:nzVal]))
        @test eltype(getproperty(C_eltype,:nzVal)) == newElty
        @test length(getproperty(C_eltype,:nzVal)) == length(getproperty(dA,:nzVal))
        @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname elty conserved" for propertyname in conserved_eltype
            @test getproperty(C_eltype,propertyname) == getproperty(dA,propertyname)
        end

        @testset "similar(A::$h{$elty},Tv::$newElty) $propertyname elty structure conserved" for propertyname in structure_conserved_eltype
            @test length(getproperty(C_eltype,propertyname)) == length(getproperty(dA,propertyname))
            @test eltype(getproperty(C_eltype,propertyname)) == eltype(getproperty(dA,propertyname))
        end

        @testset "similar(A::$h{$elty},Tv::$newElty) full" begin
            @test eltype(getproperty(C_full,:nzVal)) == newElty
            @test length(getproperty(C_full,:nzVal)) == length(getproperty(dA,:nzVal))
            @test h==CuSparseMatrixCSR ? length(getproperty(C_full,dimPtr)) == length(getproperty(dA,dimPtr)) : length(getproperty(C_full,dimPtr)) == length(getproperty(dA,dimPtr))+1
            @test getproperty(C_dims,:nnz) == getproperty(dA,:nnz)
            @test getproperty(C_full,:dims) == (n,n+1)
        end
    end

    @testset "similar(A::$f($h{$elty}),$newElty)" for elty in typeSet,
        newElty in typeSet,
        f in (transpose, x->reshape(x,n,n)),
        (h,dimPtr,dimVal) in ((CuSparseMatrixCSR,:rowPtr,:colVal), (CuSparseMatrixCSC,:colPtr,:rowVal))

        A = sprand(elty, n, n, rand())
        dA = f(h(A))

        C_simple = similar(dA)
        C_dims = similar(dA,(n,n+1))
        C_eltype = similar(dA,newElty)
        C_full = similar(dA,newElty,(n,n+1))
        @test typeof(C_simple) == typeof(parent(dA))
        @test typeof(C_dims) == typeof(parent(dA))
        @test (typeof(C_eltype) == typeof(parent(dA)) && elty==newElty) || ((typeof(C_eltype) != typeof(parent(dA)) && elty!=newElty))
        @test (typeof(C_full) == typeof(parent(dA)) && elty==newElty) || ((typeof(C_full) != typeof(parent(dA)) && elty!=newElty))

        properties = Set(propertynames(parent(dA)));

        conserved_simple = Set([:nnz])
        structure_conserved_simple = setdiff(properties,conserved_simple);
        @testset "similar(A::$f($h{$elty}),$newElty) $propertyname simple conserved" for propertyname in conserved_simple
            @test getproperty(C_simple,propertyname) == getproperty(parent(dA),propertyname)
        end
        @testset "similar(A::$f($h{$elty}),$newElty) $propertyname simple structure conserved" for propertyname in structure_conserved_simple
            @test length(getproperty(C_simple,propertyname)) == length(getproperty(parent(dA),propertyname))
            @test eltype(getproperty(C_simple,propertyname)) == eltype(getproperty(parent(dA),propertyname))
        end

        conserved_dims = Set([:nnz])
        if h==CuSparseMatrixCSR
            structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims])))
        else #CSC
            structure_conserved_dims = setdiff(properties,union(conserved_dims,Set([dimVal,:dims,dimPtr])))
            @test length(getproperty(C_dims,dimPtr)) == length(getproperty(parent(dA),dimPtr)) + 1
            @test eltype(getproperty(C_dims,dimPtr)) == eltype(getproperty(parent(dA),dimPtr))
        end
        @test getproperty(C_dims,:dims) == (n,n+1)
        @testset "similar(A::$f($h{$elty}),$newElty) $propertyname dims conserved" for propertyname in conserved_dims
            @test getproperty(C_dims,propertyname) == getproperty(parent(dA),propertyname)
        end
        @testset "similar(A::$f($h{$elty}),$newElty) $propertyname dims structure conserved" for propertyname in structure_conserved_dims
            @test length(getproperty(C_dims,propertyname)) == length(getproperty(parent(dA),propertyname))
            @test eltype(getproperty(C_dims,propertyname)) == eltype(getproperty(parent(dA),propertyname))
        end

        conserved_eltype = Set([:nnz,:dims])
        structure_conserved_eltype = setdiff(properties,union(conserved_eltype,[:nzVal]))
        @test eltype(getproperty(C_eltype,:nzVal)) == newElty
        @test length(getproperty(C_eltype,:nzVal)) == length(getproperty(parent(dA),:nzVal))
        @testset "similar(A::$f($h{$elty}),$newElty) $propertyname elty conserved" for propertyname in conserved_eltype
            @test getproperty(C_eltype,propertyname) == getproperty(parent(dA),propertyname)
        end
        @testset "similar(A::$f($h{$elty}),$newElty) $propertyname elty structure conserved" for propertyname in structure_conserved_eltype
            @test length(getproperty(C_eltype,propertyname)) == length(getproperty(parent(dA),propertyname))
            @test eltype(getproperty(C_eltype,propertyname)) == eltype(getproperty(parent(dA),propertyname))
        end
        @testset "similar(A::$f($h{$elty}),$newElty) full" begin
            @test eltype(getproperty(C_full,:nzVal)) == newElty
            @test length(getproperty(C_full,:nzVal)) == length(getproperty(parent(dA),:nzVal))
            @test h==CuSparseMatrixCSR ? length(getproperty(C_full,dimPtr)) == length(getproperty(parent(dA),dimPtr)) : length(getproperty(C_full,dimPtr)) == (length(getproperty(parent(dA),dimPtr))+1)
            @test getproperty(C_dims,:nnz) == getproperty(parent(dA),:nnz)
            @test getproperty(C_full,:dims) == (n,n+1)
        end
    end

end


text:
using DotaBot
using HTTP
using ArgParse
using Sockets
using JSON

s = ArgParseSettings()
@add_arg_table! s begin
    "--breezyIP"
    help = "breezy server IP adress"
    default = "127.0.0.1"
    "--breezyPort"
    help = "breezy server port number"
    default = "8085"
    "--agentIP"
    help = "agent server IP adress"
    default = "127.0.0.1"
    "--agentPort"
    help = "agent server port number"
    default = "8086"
end

args = parse_args(ARGS, s)

@async HTTP.serve(parse(IPAddr, args["breezyIP"]), parse(Int, args["breezyPort"])) do request::HTTP.Request
   @show request
   @show request.method
   @show HTTP.header(request, "Content-Type")
   @show HTTP.payload(request)
   try
       return HTTP.Response(200, "okay")
   catch e
       println(e)
       return HTTP.Response(404, "Error: $e")
   end
end

sleep(1)
global i = 1
while true
    if mod(i, 10) == 0
        HTTP.post(string(DotaBot.agent_ip(args), "/update"), ["Content-Type" => "application/json"],
                  JSON.json(Dict("status"=>"DONE")))
    else
        features = rand(310)
        println(features[1])
        HTTP.post(string(DotaBot.agent_ip(args), "/"),
                  ["Content-Type" => "application/json"], JSON.json(features))
    end
    sleep(0.1)
    global i += 1
end


text:
@testset "LocallyIsometricSublattice" begin

  K, a = rationals_as_number_field()
  Kt, t = PolynomialRing(K,"t")

  #
  # Split case
  #

  E, b = NumberField(t^2+7, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [3*b - 28, b + 14, 0]), map(E, [9, 3, 1]), map(E, [-21//2*b + 49//2, 5//2*b - 35//2, 0])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[3]
  @test is_split(g)
  p = g.p
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Mp, Lp, p)
  @test all(i -> Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[i].p), 1:2)

  #
  # Inert dyadic case
  #

  E, b = NumberField(t^2+3, "b", cached = false)
  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [b - 30, b - 6, 0, 0]), map(E, [-1, 1, 1, 0]), map(E, [0, 0, 0, -1]), map(E, [2*b + 15, 3, 0, 0])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[1]
  @test is_inert(g)
  p = g.p
  @test Hecke.is_dyadic(p)
  Lp = representative(g)

  @test !is_maximal(L,p)[1]
  @test_throws AssertionError Hecke.locally_isometric_sublattice(L, Lp, p)
  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[2].p)

  #
  # Inert non-dyadic case
  #

  E, b = NumberField(t^2+1, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [3*b - 17, 5*b + 11, 0]), map(E, [74*b - 412, 241//2*b + 525//2, 1//2*b - 3//2]), map(E, [-886*b + 3051, -701*b - 2077, -5*b + 11])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[1]
  @test is_inert(g)
  p = g.p
  @test !Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[2].p)
  @test genus(L) != genus(Mp) && genus(L,p) != genus(M,p) && genus(L,p) == genus(Mp,p)

  #
  # Ramified non-dyadic case
  #

  # 2 even scales less than 4

  E, b = NumberField(t^2+3, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-1//2*b - 3//2, 15//2*b + 33//2, 0]), map(E, [0, -1, -1]), map(E, [0, 3//2*b + 3//2, 0])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[1]
  @test is_ramified(g)
  @test count(s in [0,2] for s in scales(g)) == 2
  p = g.p
  @test !Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[2].p)

  # 1 even scale less than 4

  E, b = NumberField(t^2+7, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [3*b - 28, b + 14, 0]), map(E, [9, 3, 1]), map(E, [-21//2*b + 49//2, 5//2*b - 35//2, 0])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[1]
  @test is_ramified(g)
  @test count(s in [0,2] for s in scales(g)) == 1
  p = g.p
  @test !Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[2].p)

  E, b = NumberField(t^2+3, "b", cached = false)

  # 2 odd scales less than 4

  D = matrix(E, 5, 5, [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-3*b - 6, -3//2*b - 9//2, 0, 0, 0]), map(E, [-1, -1, 1, 0, 0]), map(E, [-12*b + 45, -23//2*b + 39//2, 0, -15, 0]), map(E, [0, 0, 0, 0, 1]), map(E, [14*b - 66, 15*b - 30, 0, b + 21, 0])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[1]
  @test is_ramified(g)
  @test count(s in [1,3] for s in scales(g)) == 2
  p = g.p
  @test !Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[2].p)

  # 1 odd scale less than 4

  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-1//2*b + 27//2, 3//2*b + 15//2, 0, 0]), map(E, [-1, -1, 1, 0]), map(E, [0, 0, 0, -1]), map(E, [-2*b - 30, -5*b - 15, 0, 0])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[2]
  @test is_ramified(g)
  @test count(s in [1,3] for s in scales(g)) == 1
  p = g.p
  @test !Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[1].p)

  # 1 scale greater or equal to 4

  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [3//2*b - 27//2, -1//2*b + 3//2, 0, 0]), map(E, [1, 1, -1, 0]), map(E, [-3, 0, 0, 3//2*b - 3//2]), map(E, [0, 0, 0, 1//2*b + 3//2])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[1]
  @test is_ramified(g)
  @test count(s >= 4 for s in scales(g)) != 0
  p = g.p
  @test !Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[2].p)

  #
  # Ramified dyadic case
  #

  E, b = NumberField(t^2+1, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [3*b - 17, 5*b + 11, 0]), map(E, [74*b - 412, 241//2*b + 525//2, 1//2*b - 3//2]), map(E, [-886*b + 3051, -701*b - 2077, -5*b + 11])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[2]
  @test is_ramified(g)
  p = g.p
  @test Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[1].p)

  E, b = NumberField(t^2+17, "b", cached = false)
  D = matrix(E, 2, 2, [102, b, -b, 0])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [1, 0]), map(E, [b, 0]), map(E, [0, 1]), map(E, [0, b])]
  L = hermitian_lattice(E, gens, gram = D)

  G = genus(L)
  V = L.space
  M = Hecke.maximal_integral_lattice(V)
  g = G.LGS[2]
  @test is_ramified(g)
  p = g.p
  @test Hecke.is_dyadic(p)
  Lp = representative(g)

  Mp = @inferred Hecke.locally_isometric_sublattice(M, Lp, p)
  @test is_sublattice(M, Mp)
  @test Hecke._islocally_isometric_kirschmer(Lp, Mp, p)
  @test !Hecke._islocally_isometric_kirschmer(M, Mp, p)
  @test Hecke._islocally_isometric_kirschmer(M, Mp, G.LGS[1].p)

end


text:
function init_pseudo_observable_mappings!(m::Model990)

    pseudo_names = [:y_t, :y_f_t, :NaturalRate, :π_t, :OutputGap, :ExAnteRealRate, :LongRunInflation,
                    :MarginalCost, :Wages, :FlexibleWages, :Hours, :FlexibleHours, :z_t, :Z_t, :ztil_t,
                    :NominalFFR, :RealFFR, :ExpectedNominalNaturalRate]

    # Create PseudoObservable objects
    pseudo = OrderedDict{Symbol,PseudoObservable}()
    for k in pseudo_names
        pseudo[k] = PseudoObservable(k)
    end

    # Fill in names and reverse transforms
    pseudo[:y_t].name = "Output Growth"
    pseudo[:y_t].longname = "Output Growth Per Capita"

    pseudo[:y_f_t].name = "Flexible Output Growth"
    pseudo[:y_f_t].longname = "Output that would obtain in a flexible-price economy."

    pseudo[:NaturalRate].name = "Real Natural Rate"
    pseudo[:NaturalRate].longname = "The real interest rate that would prevail in a flexible-price economy."
    pseudo[:NaturalRate].rev_transform = quartertoannual

    pseudo[:π_t].name = "Inflation"
    pseudo[:π_t].longname = "Inflation"
    pseudo[:π_t].rev_transform = quartertoannual

    pseudo[:OutputGap].name = "Output Gap"
    pseudo[:OutputGap].longname = "Output Gap"

    pseudo[:ExAnteRealRate].name = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].longname = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].rev_transform = quartertoannual

    pseudo[:LongRunInflation].name = "Long Run Inflation"
    pseudo[:LongRunInflation].longname = "Long Run Inflation"
    pseudo[:LongRunInflation].rev_transform = quartertoannual

    pseudo[:MarginalCost].name = "Marginal Cost"
    pseudo[:MarginalCost].longname = "Marginal Cost"

    pseudo[:Wages].name = "Wages"
    pseudo[:Wages].longname = "Wages"

    pseudo[:FlexibleWages].name = "Flexible Wages"
    pseudo[:FlexibleWages].longname = "Wages that would obtain in a flexible-wage economy"

    pseudo[:Hours].name = "Hours"
    pseudo[:Hours].longname = "Hours"

    pseudo[:FlexibleHours].name     = "Flexible Hours"
    pseudo[:FlexibleHours].longname = "Flexible Hours"

    pseudo[:z_t].name     = "z_t"
    pseudo[:z_t].longname = "z_t"

    pseudo[:Z_t].name     = "Z_t"
    pseudo[:Z_t].longname = "Z_t"

    pseudo[:z_t].name     = "ztil_t"
    pseudo[:z_t].longname = "ztil_t"

    pseudo[:NominalFFR].name     = "Nominal FFR"
    pseudo[:NominalFFR].longname = "Nominal FFR at an annual rate"
    pseudo[:NominalFFR].rev_transform = quartertoannual

    pseudo[:RealFFR].name     = "Real FFR"
    pseudo[:RealFFR].longname = "Real FFR at an annual rate"
    pseudo[:RealFFR].rev_transform = quartertoannual

    pseudo[:ExpectedNominalNaturalRate].name  = "Expected Nominal Natural Rate"
    pseudo[:ExpectedNominalNaturalRate].longname = "Natural Rate + Expected Inflation"
    pseudo[:ExpectedNominalNaturalRate].rev_transform = quartertoannual

    # Add to model object
    m.pseudo_observable_mappings = pseudo
end


text:
module TaskTests
include("test_daemontask.jl")
include("test_basictaskinfo.jl")
end # module TaskTests

max_forks_repo_name: seung-lab/SimpleTasks


text:
function modinv{T<:Integer}(a::T, b::T)
    b0 = b
    x0, x1 = zero(T), one(T)
    while a > 1
        q = div(a, b)
        a, b = b, a % b
        x0, x1 = x1 - q * x0, x0
    end
    x1 < 0 ? x1 + b0 : x1
end
modinv(a::Integer, b::Integer) = modinv(promote(a,b)...)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    DynamicAxis(n::Int)

A range spanning `1` to `n` by steps of `1`. The value of `n` can be changed after
construction, but cannot be any less than `0`.
"""
mutable struct DynamicAxis <: AbstractUnitRange{Int}
    stop::Int

    DynamicAxis(stop::Real) = new(Int(max(0, stop)))

    function DynamicAxis(r::AbstractRange)
        first(r) == 1 || (Base.@_noinline_meta; throw(ArgumentError("first element must be 1, got $(first(r))")))
        step(r)  == 1 || (Base.@_noinline_meta; throw(ArgumentError("step must be 1, got $(step(r))")))
        return DynamicAxis(last(r))
    end
end

Base.first(::DynamicAxis) = 1

Base.last(r::DynamicAxis) = getfield(r, :stop)
Base.length(x::DynamicAxis) = last(x)
Base.lastindex(r::DynamicAxis) = last(r)::Int

Base.issubset(r::DynamicAxis, s::OneTo) = last(r) <= last(s)
Base.issubset(r::DynamicAxis, s::DynamicAxis) = last(r) <= last(s)
Base.issubset(r::OneTo, s::DynamicAxis) = last(r) <= last(s)

Base.mod(i::Integer, r::DynamicAxis) = Base.mod1(i, last(r))

function Base.setproperty!(x::DynamicAxis, s::Symbol, val)
    error("cannot use setproperty! on DynamicAxis")
end

Base.intersect(r::DynamicAxis, s::DynamicAxis) = OneTo(min(last(r),last(s)))
Base.intersect(r::DynamicAxis, s::OneTo) = OneTo(min(last(r),last(s)))
Base.intersect(r::OneTo, s::DynamicAxis) = OneTo(min(last(r),last(s)))

Base.show(io::IO, r::DynamicAxis) = print(io, "DynamicAxis($(last(r)))")

Base.:(==)(r::DynamicAxis, s::DynamicAxis) = last(r) === last(s)
Base.:(==)(r::OneTo, s::DynamicAxis) = last(r) === last(s)
Base.:(==)(r::DynamicAxis, s::OneTo) = last(r) == last(s)

Base.empty!(r::DynamicAxis) = (setfield!(r, :stop, 0); r)

Base.in(x::Integer, r::DynamicAxis) = !(1 > x) & !(x > last(r))

function Base.in(x::Real, r::DynamicAxis)
    val = round(Integer, x)
    if in(val, r)
        return @inbounds(getindex(r, val)) == x
    else
        return false
    end
end

function Base.Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::DynamicAxis, x::Number)
    return range(first(r)-x, length=length(r))
end

Base.Broadcast.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::DynamicAxis) = r

###
### ArrayInterface
###
ArrayInterface.can_change_size(::Type{DynamicAxis}) = true
ArrayInterface.known_first(::Type{DynamicAxis}) = 1
ArrayInterface.known_step(::Type{DynamicAxis}) = 1

@propagate_inbounds function Base.getindex(v::DynamicAxis, i::Integer)
    @boundscheck ((i > 0) & (i <= last(v))) || throw(BoundsError(v, i))
    return Int(i)
end

@propagate_inbounds function Base.getindex(r::DynamicAxis, s::OneTo)
    @boundscheck checkbounds(r, s)
    return OneTo(Int(last(s)))
end

Base.AbstractUnitRange{Int}(r::DynamicAxis) = DynamicAxis(r)

Base.promote_rule(::Type{DynamicAxis}, ::Type{OneTo{T}}) where {T} = OneTo{promote_type(T,Int)}
Base.promote_rule(::Type{OneTo{T}}, ::Type{DynamicAxis}) where {T} = OneTo{promote_type(T,Int)}

#=
Base.promote_rule(a::Type{UnitMRange{T1}}, b::Type{DynamicAxis{T2}}) where {T1,T2} = UnitMRange{promote_type(T1,T2)}
Base.promote_rule(a::Type{DynamicAxis}, b::Type{UnitMRange{T2}}) where {T1,T2} = UnitMRange{promote_type(T1,T2)}
=#

function unsafe_grow_end!(x::DynamicAxis, n)
    setfield!(x, :stop, Int(last(x) + n))
    return x
end

function unsafe_shrink_end!(x::DynamicAxis, n)
    setfield!(x, :stop, Int(last(x) - n))
    return x
end


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

# Implementation(s) provided in separate file, loaded via require:
function _h5io_keys end
function _h5io_objtype end
function _h5io_open end

function _h5io_add_path_to_dest(dest_with_subpath::Tuple{Any,AbstractString}, path::AbstractString)
    dest, old_path = dest_with_subpath
    @assert !isempty(old_path)
    @assert !isempty(path)
    (dest, "$old_path/$path")
end

_h5io_add_path_to_dest(dest, path::AbstractString) = (dest, path)

function _h5io_write!(dest_with_subpath::Tuple{Any,AbstractString}, data::AbstractArray{<:Real})
    dest, path = dest_with_subpath
    if occursin('/', path) && _h5_track_order_available()
        group_name = dirname(path)
        dataset_name = basename(path)
        g = if !haskey(dest, group_name)
            HDF5.create_group(dest, group_name; track_order = true)
        else
            HDF5.open_group(dest, group_name)
        end
        g[dataset_name] = data
    else
        dest[path] = data
    end
    nothing
end

function _h5io_write!(dest_with_subpath::Tuple{Any,AbstractString}, data::AbstractArray{<:Nothing})
    nothing
end

function _h5io_write!(dest_with_subpath::Tuple{Any,AbstractString}, data::VectorOfSimilarVectors)
    _h5io_write!(dest_with_subpath, Array(flatview(data)))
end

function _h5io_write!(dest, data::NamedTuple{names}) where {names}
    for name in names
        _h5io_write!(_h5io_add_path_to_dest(dest, String(name)), Base.getproperty(data, name))
    end
    nothing
end

function _h5io_write!(dest, data::FillArrays.Fill{<:Real,1})
    _h5io_write!(dest, convert(Array, data))
end

function _h5io_write!(dest, data::FillArrays.Fill{<:AbstractVector{<:Real},1})
    _h5io_write!(dest, convert(VectorOfSimilarVectors, data))
end

function _h5io_write!(dest, data)
    Tables.istable(data) || throw(ArgumentError("data is not a table."))
    cols = Tables.columns(data)::NamedTuple
    @assert cols isa NamedTuple
    _h5io_write!(dest, cols)
end

function _h5io_read end

function _h5io_read(src_with_subpath::Tuple{Any,AbstractString})
    src, subpath = src_with_subpath
    _h5io_read(src[subpath])
end

_h5io_read(src) = _h5io_read(src, _h5io_objtype(src))

_h5io_read(src, ::Val{:dataset}) = _h5io_read_postprocess(read(src))

function _h5io_read(src, ::Val{:datafile})
    names = (_h5io_keys(src)...,)
    values = map(k -> _h5io_read((src, k)), names)
    name_symbols = map(Symbol, names)
    x = NamedTuple{name_symbols}(values)
    _h5io_read_postprocess(x)
end

_h5io_read_postprocess(V::AbstractVector{<:Real}) = V

_h5io_read_postprocess(VV::AbstractMatrix{<:Real}) = VectorOfSimilarVectors(VV)

_h5io_read_postprocess(nt::NamedTuple) = TypedTables.Table(nt)

function _h5io_read_postprocess(nt::NamedTuple{(:v, :logd, :weight, :info)})
    DensitySampleVector((
        _h5io_read_postprocess_samples(nt.v),
        nt.logd, nt.weight, nt.info, Array{Nothing}(undef, size(nt.info)...))
    )
end
function _h5io_read_postprocess(nt::NamedTuple{(:info, :logd, :v, :weight)})
    # This method is needed for older hdf5 files where `track_order` was not yet used
    DensitySampleVector((
        _h5io_read_postprocess_samples(nt.v),
        nt.logd, nt.weight, nt.info, Array{Nothing}(undef, size(nt.info)...))
    )
end

# Column :info will be missing if `eltype(samples.info)` was `Nothing`:
function _h5io_read_postprocess(nt::NamedTuple{(:v, :logd, :weight)})
    _h5io_read_postprocess(merge(nt, (info = Array{Nothing}(undef, size(nt.weight)...),)))
end
function _h5io_read_postprocess(nt::NamedTuple{(:logd, :v, :weight)})
    # This method is needed for older hdf5 files where `track_order` was not yet used
    _h5io_read_postprocess(merge((info = Array{Nothing}(undef, size(nt.weight)...),), nt))
end

_h5io_read_postprocess(nt::NamedTuple{(:chainid, :chaincycle, :stepno, :sampletype)}) =
    MCMCSampleIDVector((nt.chainid, nt.chaincycle, nt.stepno, nt.sampletype))

# This method is needed for older hdf5 files where `track_order` was not yet used
h5io_read_postprocess(nt::NamedTuple{(:chaincycle, :chainid, :sampletype, :stepno)}) =
    MCMCSampleIDVector((nt.chainid, nt.chaincycle, nt.stepno, nt.sampletype))

function _const_col_value(col::AbstractVector)
    r = first(col)
    all(isequal(r), col) ? r : missing
end

_normalize_vs(shape::AbstractValueShape) = shape

function _normalize_vs(shape::NamedTupleShape{names,AT,VT}) where {names,AT,VT}
    NamedTupleShape(VT, map(_normalize_vs, (;shape...)))
end

_normalize_vs(shape::ScalarShape{<:Real}) = ScalarShape{Real}()
_normalize_vs(shape::ArrayShape{<:Real}) = ArrayShape{Real}(size(shape)...)

function _infer_vs_from_table(v::TypedTables.Table)
    cols = Tables.columns(v)
    raw_vs = map(valshape, map(first, cols))
    const_v = map(_const_col_value, cols)
    NamedTupleShape(map((c, shp) -> ismissing(c) ? _normalize_vs(shp) : ConstValueShape(copy(deepcopy(c))), const_v, (;raw_vs...)))
end

_h5io_read_postprocess_samples(v::AbstractVector) = v

function _h5io_read_postprocess_samples(v::TypedTables.Table)
    shp = _infer_vs_from_table(v)
    unshaped_v = VectorOfSimilarVectors(unshaped.(v, Ref(shp)))
    shp.(unshaped_v)
end


text:
"""
    get_constrained_lowdimset(cpa::CartesianProductArray{N, S}, P::AbstractPolyhedron{N}) where {N<:Real, S<:LazySet{N}}

Preprocess step for intersection between cartesian product array and H-polyhedron.
Returns low-dimensional HPolytope in constrained dimensions of original cpa,
constrained variables and variables in corresponding blocks, original block structure
of low-dimensional set and list of constrained blocks.

### Input

- `cpa` -- Cartesian product array of convex sets
- `P` -- polyhedron

### Output

A tuple of low-dimensional set, list of constrained dimensions, original block
structure of low-dimensional set and corresponding blocks indices.
"""
function get_constrained_lowdimset(cpa::CartesianProductArray{N, S}, P::AbstractPolyhedron{N}) where {N<:Real, S<:LazySet{N}}

    if isbounded(P)
        blocks, non_empty_length = block_to_dimension_indices(cpa)
    else
        constrained_vars = constrained_dimensions(P)
        blocks, non_empty_length = block_to_dimension_indices(cpa, constrained_vars)
    end

    array = Vector{S}()
    sizehint!(array, non_empty_length)
    vars = Vector{Int}()
    block_structure = Vector{UnitRange{Int}}()
    sizehint!(block_structure, non_empty_length)

    last_var = 1
    for i in 1:length(blocks)
        start_index, end_index = blocks[i]
        block_end = last_var + end_index - start_index
        if start_index != -1
            push!(array, cpa.array[i])
            append!(vars, start_index : end_index)
            push!(block_structure, last_var : block_end)
            last_var = block_end + 1
        end
    end

    cpa_low_dim = HPolytope(constraints_list(CartesianProductArray(array)));

    return cpa_low_dim, vars, block_structure, blocks
end


text:
# ---
# title: 1524. Number of Sub-arrays With Odd Sum
# id: problem1524
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Math
# link: <https://leetcode.com/problems/number-of-sub-arrays-with-odd-sum/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`. Return _the number of sub-arrays_ with
# **odd** sum.
#
# As the answer may grow large, the answer  **must be**  computed modulo `10^9 +
# 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,3,5]
#     Output: 4
#     Explanation: All sub-arrays are [[1],[1,3],[1,3,5],[3],[3,5],[5]]
#     All sub-arrays sum are [1,4,9,3,8,5].
#     Odd sums are [1,9,3,5] so the answer is 4.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [2,4,6]
#     Output: 0
#     Explanation: All sub-arrays are [[2],[2,4],[2,4,6],[4],[4,6],[6]]
#     All sub-arrays sum are [2,6,12,4,10,6].
#     All sub-arrays have even sum and the answer is 0.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [1,2,3,4,5,6,7]
#     Output: 16
#
#
# **Example 4:**
#
#
#
#     Input: arr = [100,100,99,99]
#     Output: 4
#
#
# **Example 5:**
#
#
#
#     Input: arr = [7]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `1 <= arr[i] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    LegendreUnitCoeff{T}

Precomputed recursion relation coefficients for the standard unit
normalization. Alias for `LegendreNormCoeff{LegendreUnitNorm,T}`.
"""
LegendreUnitCoeff{T} = LegendreNormCoeff{LegendreUnitNorm,T}

"""
    LegendreFourPiCoeff{T}

Table type of precomputed recursion relation coefficients for the ``4\\pi`` spherical
harmonic normalization. Alias for `LegendreNormCoeff{LegendreFourPiNorm,T}`.
"""
LegendreFourPiCoeff{T} = LegendreNormCoeff{LegendreFourPiNorm,T}

"""
    LegendreOrthoCoeff{T}

Table type of precomputed recursion relation coefficients for the orthonormal
normalization. Alias for `LegendreNormCoeff{LegendreOrthoNorm,T}`.
"""
LegendreOrthoCoeff{T} = LegendreNormCoeff{LegendreOrthoNorm,T}

"""
    LegendreSphereCoeff{T}

Table type of precomputed recursion relation coefficients for the spherical
harmonic normalization. Alias for `LegendreNormCoeff{LegendreSphereNorm,T}`.
"""
LegendreSphereCoeff{T} = LegendreNormCoeff{LegendreSphereNorm,T}

"""
    p = Plm(l, m, x)

Computes the associated Legendre polynomials using unit normalization;
equivalent to `p = legendre(LegendreUnitNorm(), l, m, x)`.
"""
const Plm = LegendreUnitNorm()

"""
    Plm!(P, l, m, x)

Fills the array `P` with the unit-normalized associated Legendre polynomial values
``P_ℓ^m(x)``;
equivalent to `legendre!(LegendreUnitNorm(), P, l, m, x)`.
"""
const Plm! = LegendreUnitNorm()

"""
    λ = λlm(l, m, x)

Computes the associated Legendre polynomials using spherical-harmonic normalization;
equivalent to `λ = legendre(LegendreSphereNorm(), l, m, x)`.
"""
const λlm = LegendreSphereNorm()

"""
    λlm!(Λ, l, m, x)

Fills the array `Λ` with the spherical-harmonic normalized associated Legendre polynomial
values ``λ_ℓ^m(x)``;
equivalent to `legendre!(LegendreSphereNorm(), P, l, m, x)`.
"""
const λlm! = LegendreSphereNorm()