text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames,StatsPlots,FileIO,MAT,RecursiveArrayTools
using Revise
import KenyaCoV
using LinearAlgebra:eigen
using Statistics: median, quantile

"""
Consensus modelling
"""

"""
Load age structured data, define callback control measures, and effect of regional lockdown
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Put in the lockdown effects
T_normal = deepcopy(P.T)
T_regional_lockdown = deepcopy(P.T)
#Outgoing travel
#Nairobi
T_regional_lockdown[:,4] *= 0.1;T_regional_lockdown[4,4] += 1 - sum(T_regional_lockdown[:,4])
#Mombasa + Kwale
T_regional_lockdown[:,12] *= 0.1;T_regional_lockdown[12,12] += 1 - sum(T_regional_lockdown[:,12])
#Kilifi + Kwale
T_regional_lockdown[:,20] *= 0.1;T_regional_lockdown[20,20] += 1 - sum(T_regional_lockdown[:,20])
#Incoming travel
for leaving_area in 1:20,arriving_area in 1:20
    if !(leaving_area in [4,12,20]) && arriving_area in [4,12,20]
        amount_reduced = 0.9*T_regional_lockdown[arriving_area,leaving_area]
        T_regional_lockdown[arriving_area,leaving_area] -= amount_reduced
        T_regional_lockdown[leaving_area,leaving_area] += amount_reduced #All avoided trips to locked down areas lead to staying at home
    end
end

@load "data/detection_rates_for_different_taus.jld2" d_0 d_01 d_025 d_05 d_1
@load "data/susceptibility_rates.jld2" σ

#Initial infecteds
"""
Load age mixing matrices (these are all in to (row) from (col) format)
"""

@load "data/agemixingmatrix_Kenya_all_types.jld2" M_Kenya M_Kenya_ho M_Kenya_other M_Kenya_school M_Kenya_work

#Function for changing contacts so as to have -45% over 14 days
function ramp_down(t)
    if t < 0.
        return 1.
    elseif t>= 0. && t <= 14.
        return (1-t/14.) + 0.55*t/14.
    elseif t > 14.
        return 0.55
    end
end

#Put in the regional lockdown on day 25 (April 7th)
function regional_lockdown_timing(u,t,integrator)
  !integrator.p.lockdown && t > 25.
end
function affect_regional_lockdown!(integrator)
  integrator.p.T = T_regional_lockdown
  integrator.p.lockdown = true
end
cb_regional_lockdown = DiscreteCallback(regional_lockdown_timing,affect_regional_lockdown!)

# both_cbs = CallbackSet(cb_iso_limit,cb_SD_limit)

"""
SCENARIO no intervention baseline
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_1
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 1.
#Set the susceptibility vector --- just to specify the correct R₀
sus_matrix = repeat(P.χ,1,17)
R_A = P.ϵ*((1/P.σ₂) + (1/P.γ) ) #effective duration of asymptomatic
R_M = (P.ϵ/P.σ₂) + (P.ϵ_D/P.γ) #effective duration of mild
R_V = (P.ϵ/P.σ₂) + (P.ϵ_V/P.τ) #effective duration of severe
R_vector = [(1-P.rel_detection_rate[a])*R_A + P.rel_detection_rate[a]*(1-P.hₐ[a])*R_M + P.rel_detection_rate[a]*P.hₐ[a]*R_V for a = 1:17]
inf_matrix = repeat(R_vector',17,1)

eigs, = eigen(sus_matrix.*P.M.*inf_matrix)
max_eigval = Real(eigs[end])
P.χ = ones(KenyaCoV.n_a)/max_eigval #This rescales everything so β is the same as R₀

P.β = rand(KenyaCoV.d_R₀) #Choose R₀ randomly from 2-3 range

u0[4,8,3] = 30 #10 initial pre-symptomatics in Nairobi
u0[12,8,3] = 10 #10 initial pre-symptomatics in Mombasa

prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,1*365.),P)

sims_baseline = KenyaCoV.run_consensus_simulations(P::KenyaCoV.CoVParameters_AS,prob,1000,CallbackSet())
P.χ .*= 1.384
sims_baseline_scaled = KenyaCoV.run_consensus_simulations(P::KenyaCoV.CoVParameters_AS,prob,1000,CallbackSet())

@save joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_baseline.jld2") sims_baseline
@save joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_baseline_scaled.jld2") sims_baseline_scaled

println("Finished baseline sims for consensus modelling ")

"""
SCENARIO 2 --- regional lockdowns
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_1
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 1.
#Set the susceptibility vector --- just to specify the correct R₀
sus_matrix = repeat(P.χ,1,17)
R_A = P.ϵ*((1/P.σ₂) + (1/P.γ) ) #effective duration of asymptomatic
R_M = (P.ϵ/P.σ₂) + (P.ϵ_D/P.γ) #effective duration of mild
R_V = (P.ϵ/P.σ₂) + (P.ϵ_V/P.τ) #effective duration of severe
R_vector = [(1-P.rel_detection_rate[a])*R_A + P.rel_detection_rate[a]*(1-P.hₐ[a])*R_M + P.rel_detection_rate[a]*P.hₐ[a]*R_V for a = 1:17]
inf_matrix = repeat(R_vector',17,1)

eigs, = eigen(sus_matrix.*P.M.*inf_matrix)
max_eigval = Real(eigs[end])
P.χ = ones(KenyaCoV.n_a)/max_eigval

P.β = rand(KenyaCoV.d_R₀) #Choose R₀ randomly from mean 2.5 (2-3) range
P.c_t = ramp_down #This implements the social distancing over 14 days from time 0.

u0[4,8,3] = 30 #30 initial Asymptomatics in Nairobi
u0[12,8,3] = 10 #10 initial pre-symptomatics in Mombasa
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,1*365.),P)

sims_controls = KenyaCoV.run_consensus_simulations(P::KenyaCoV.CoVParameters_AS,prob,1000,cb_regional_lockdown)
P.χ .*= 1.384 #This accounts for the difference in scale between China and Kenya of basic contact rate
sims_controls_scaled = KenyaCoV.run_consensus_simulations(P::KenyaCoV.CoVParameters_AS,prob,1000,cb_regional_lockdown)

@save joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_control.jld2") sims_controls
@save joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_control_scaled.jld2") sims_controls_scaled


text:
module entities

abstract type Company end

struct Supplier <: Company
    name::String
end

struct Manufacturer <: Company
    name::String
end

struct Customer <: Company
    name::String
end

end # Module entities

ext: jl


text:
function main()
    a = max(1, 2)
    println(join([a], " "))
    b = min(1, 2)
    println(join([b], " "))
end

main()


text:
import LinearAlgebra
import NonlinearEquations
import PyPlot
import SparseArrays

function analytical(x, t)
	if x == 0
		return 1.0
	elseif x > ((0.5 + 1 / sqrt(2)) * t)
		return 0.0
	else
		return 0.5 * (sqrt(-2 + t * (sqrt(4 * x / t + 1) - 1) / x + 1) + 1)
	end
end

function getfreenodes(n, dirichletnodes)
	isfreenode = fill(true, n)
	isfreenode[dirichletnodes] .= false
	nodei2freenodei = fill(-1, length(isfreenode))
	freenodei2nodei = Array{Int}(undef, sum(isfreenode))
	j = 1
	for i = 1:length(isfreenode)
		if isfreenode[i]
			nodei2freenodei[i] = j
			freenodei2nodei[j] = i
			j += 1
		end
	end
	return isfreenode, nodei2freenodei, freenodei2nodei
end

f(S) = S ^ 2 / (S ^ 2 + (1 - S) ^ 2)
Ss = 0:1e-2:1
fig, ax = PyPlot.subplots()
ax.plot(Ss, f.(Ss))
ax.set_xlabel("S")
ax.set_ylabel("Fractional Flow")
display(fig)
println()
PyPlot.close(fig)

@NonlinearEquations.equations function bl(Sfree, p, Ks, neighbors, dirichletnodes, dirichletS, area, xs, isfreenode, freenodei2nodei, nodei2freenodei)
	numfreenodes = sum(isfreenode)
	NonlinearEquations.setnumequations(numfreenodes)
	for (i, (node_a, node_b)) in enumerate(neighbors)
		for (node1, node2) in [(node_a, node_b), (node_b, node_a)]
			if isfreenode[node1] && isfreenode[node2]
				j1 = nodei2freenodei[node1]
				j2 = nodei2freenodei[node2]
				if p[node2] > p[node1]
					NonlinearEquations.addterm(j1, Ks[i] * ((p[node2] - p[node1]) / (xs[node2] - xs[node1])) * (Sfree[j2] ^ 2 / (Sfree[j2] ^ 2 + (1 - Sfree[j2]) ^ 2) - Sfree[j1] ^ 2 / (Sfree[j1] ^ 2 + (1 - Sfree[j1]) ^ 2)) * area / (xs[node2] - xs[node1]))
				end
			elseif isfreenode[node1] && !isfreenode[node2]
				j1 = nodei2freenodei[node1]
				if p[node2] > p[node1]
					NonlinearEquations.addterm(j1, Ks[i] * ((p[node2] - p[node1]) / (xs[node2] - xs[node1])) * (dirichletS[node2] ^ 2 / (dirichletS[node2] ^ 2 + (1 - dirichletS[node2]) ^ 2) - Sfree[j1] ^ 2 / (Sfree[j1] ^ 2 + (1 - Sfree[j1]) ^ 2)) * area / (xs[node2] - xs[node1]))
				end
			end
		end
	end
end

function bl_simple(S, f=S->S^2/(S^2+(1-S)^2), area=1)
	n = length(S)
	dx = 1 / n
	dS = zeros(n)
	dS[1] = (f(1) - f(S[1])) * area / dx
	for i = 2:length(S)
		dS[i] = (f(S[i - 1]) - f(S[i])) * area / dx
	end
	return dS
end

n = 1000
neighbors = [i=>i + 1 for i = 1:n - 1]
Ks = ones(length(neighbors))
xs = 0:1/(n - 1):1
S = zeros(n)
dirichletnodes = [1]
dirichletS = zeros(length(S))
dirichletS[1] = 1.0
S[dirichletnodes] = dirichletS[dirichletnodes]
isfreenode, nodei2freenodei, freenodei2nodei = getfreenodes(length(S), dirichletnodes)
Sfree = S[isfreenode]
Sfree_implicit = S[isfreenode]
dt = 0.5 / n
finalt = 0.5
ts = 0:dt:finalt
analytical_solution = map(x->analytical(x, finalt), xs)
p = map(x->1 - x, xs)
area = 1
for t in ts
	global r1, r2
	r1 = bl_simple(Sfree)
	r2 = bl_residuals(Sfree, p, Ks, neighbors, dirichletnodes, dirichletS, area, xs, isfreenode, freenodei2nodei, nodei2freenodei)
	if !(r1 ≈ r2)
		@show t, sum((r1 .- r2) .^ 2)
		error("blah")
	end
	Sfree .+= dt * r2
	if any(isnan.(Sfree))
		@show t
		error("nan")
	end
	Sfree_implicit .+= (LinearAlgebra.I - dt * bl_Sfree(Sfree_implicit, p, Ks, neighbors, dirichletnodes, dirichletS, area, xs, isfreenode, freenodei2nodei, nodei2freenodei)) \ (dt * bl_residuals(Sfree_implicit, p, Ks, neighbors, dirichletnodes, dirichletS, area, xs, isfreenode, freenodei2nodei, nodei2freenodei))
end
S[isfreenode] = Sfree
S_implicit = copy(S)
S_implicit[isfreenode] = Sfree_implicit
fig, ax = PyPlot.subplots()
ax.plot(xs, S, label="explicit")
ax.plot(xs, S_implicit, label="implicit")
ax.plot(xs, analytical_solution, label="analytical")
ax.legend()
ax.set_xlabel("x")
ax.set_xlabel("S")
fig.tight_layout()
display(fig)
println()
PyPlot.close(fig)

n = 1000
neighbors = [i=>i + 1 for i = 1:n - 1]
Ks = ones(length(neighbors))
xs = 0:1/(n - 1):1
S = zeros(n)
dirichletnodes = [length(S)]
dirichletS = zeros(length(S))
dirichletS[end] = 1.0
S[dirichletnodes] = dirichletS[dirichletnodes]
isfreenode, nodei2freenodei, freenodei2nodei = getfreenodes(length(S), dirichletnodes)
Sfree = S[isfreenode]
Sfree_implicit = S[isfreenode]
dt = 0.5 / n
finalt = 0.5
ts = 0:dt:finalt
analytical_solution = map(x->analytical(1 - x, finalt), xs)
p = map(x->x, xs)
area = 1
for t in ts
	r2 = bl_residuals(Sfree, p, Ks, neighbors, dirichletnodes, dirichletS, area, xs, isfreenode, freenodei2nodei, nodei2freenodei)
	Sfree .+= dt * r2
	if any(isnan.(Sfree))
		@show t
		error("nan")
	end
	Sfree_implicit .+= (LinearAlgebra.I - dt * bl_Sfree(Sfree_implicit, p, Ks, neighbors, dirichletnodes, dirichletS, area, xs, isfreenode, freenodei2nodei, nodei2freenodei)) \ (dt * bl_residuals(Sfree_implicit, p, Ks, neighbors, dirichletnodes, dirichletS, area, xs, isfreenode, freenodei2nodei, nodei2freenodei))
end
S[isfreenode] = Sfree
S_implicit = copy(S)
S_implicit[isfreenode] = Sfree_implicit
fig, ax = PyPlot.subplots()
ax.plot(xs, S, label="explicit")
ax.plot(xs, S_implicit, label="implicit")
ax.plot(xs, analytical_solution, label="analytical")
ax.legend()
ax.set_xlabel("x")
ax.set_xlabel("S")
fig.tight_layout()
display(fig)
println()
PyPlot.close(fig)


text:
using ProbNumODE
using ProbNumODE: remake_prob_with_jac, multitime
using DifferentialEquations
using Statistics
using LinearAlgebra
using StatsPlots
using LaTeXStrings

FILEDIR = dirname(@__FILE__)

for (prob, name) in [
    (remake_prob_with_jac(lotka_volterra()), "lotkavolterra"),
    (remake_prob_with_jac(fitzhugh_nagumo_iip()), "fitzhughnagumo"),
    (remake_prob_with_jac(logistic_equation()), "logistic"),
    (remake_prob_with_jac(brusselator()), "brusselator"),
]
    appxsol = solve(remake(prob, u0=big.(prob.u0)), abstol=1e-30, reltol=1e-30)

    densetimes = collect(range(prob.tspan[1], stop=prob.tspan[end], length=1000))
    dense_ref = appxsol.(densetimes)

    function get_res(solve_call)
        sol = solve_call()
        runtime = multitime(solve_call; numruns=10)
        errs = sol.(densetimes) - dense_ref
        L2 = norm(sqrt.(mean.([float.(d) .^ 2 for d in errs])))
        nevals = sol.destats.nf + sol.destats.njacs
        return L2, nevals, runtime
    end

    ekf0_results = []
    ekf0_smooth_results = []
    ekf1_results = []
    ekf1_smooth_results = []
    dp5_results = []
    for (i, (_a, _r)) in enumerate(zip(7:17, 2:14))
        @info "Iteration" name _a _r
        solve_call = () -> solve(prob, EKF0(), sigmarule=:dynamicMV, q=5, abstol=10.0^-_a, reltol=10.0^-_r, smooth=false)
        push!(ekf0_results, get_res(solve_call))

        solve_call = () -> solve(prob, EKF0(), sigmarule=:dynamicMV, q=5, abstol=10.0^-_a, reltol=10.0^-_r, smooth=true)
        push!(ekf0_smooth_results, get_res(solve_call))

        solve_call = () -> solve(prob, EKF1(), sigmarule=:dynamic, q=5, abstol=10.0^-_a, reltol=10.0^-_r, smooth=false)
        push!(ekf1_results, get_res(solve_call))

        solve_call = () -> solve(prob, EKF1(), sigmarule=:dynamic, q=5, abstol=10.0^-_a, reltol=10.0^-_r, smooth=true)
        push!(ekf1_smooth_results, get_res(solve_call))

        solve_call = () -> solve(prob, DP5(), abstol=10.0^-_a, reltol=10.0^-_r)
        push!(dp5_results, get_res(solve_call))
    end

    p = plot([r[2] for r in ekf0_results],
             [r[1] for r in ekf0_results],
             xscale=:log10,
             yscale=:log10,
             marker=:auto,
             label="EKF0",
             );
    plot!(p, [r[2] for r in ekf0_smooth_results],
          [r[1] for r in ekf0_smooth_results],
          marker=:auto,
          label="EKS0",
          );
    plot!(p, [r[2] for r in ekf1_results],
          [r[1] for r in ekf1_results],
          marker=:auto,
          label="EKF1",
          );
    plot!(p, [r[2] for r in ekf1_smooth_results],
          [r[1] for r in ekf1_smooth_results],
          marker=:auto,
          label="EKS1",
          );
    plot!(p, [r[2] for r in dp5_results],
          [r[1] for r in dp5_results],
          marker=:auto,
          label="DP5",
          );

    x0_dp5 = dp5_results[1][2]
    y0_dp5 = dp5_results[1][1]
    xmax_dp5 = dp5_results[end][2]
    xrange = 10 .^ range(log10(x0_dp5), log10(xmax_dp5), length=1000)
    plot!(p, xrange, x -> x^(-6)*(x0_dp5^6 * y0_dp5), color=:black, linestyle=:dot, label=L"$h^6$ conv.")

    x0 = ekf1_results[1][2]
    y0 = ekf1_results[1][1]
    xrange = 10 .^ range(log10(x0), log10(xmax_dp5), length=1000)
    plot!(p, xrange, x -> x^(-6)*(x0^6 * y0), color=:black, linestyle=:dot, label="")

    x0 = ekf1_smooth_results[1][2]
    y0 = ekf1_smooth_results[1][1]
    xrange = 10 .^ range(log10(x0), log10(xmax_dp5), length=1000)
    plot!(p, xrange, x -> x^(-6)*(x0^6 * y0), color=:black, linestyle=:dot, label="")

    plot!(p, xlabel=L"\#evaluations", ylabel="L2 error")
    savefig(p, joinpath(FILEDIR, "dp5_comparison_$(name).png"))
end

max_stars_repo_path: experiments/3 Comparison to DP5/comparison_dp5.jl

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
@params mutable struct TrussStress{T} <: AbstractFunction{T}
    σ::AbstractVector{T} # stress vector, axial stress per cell
    u_fn::Displacement
    transf_matrices::AbstractVector{<:AbstractMatrix{T}}
    fevals::Int
    maxfevals::Int
end

function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::TrussStress)
    return println("TopOpt truss stress function")
end

"""
    TrussStress(solver; maxfevals=10^8)

Construct the TrussStress function struct.
"""
function TrussStress(solver::AbstractFEASolver; maxfevals=10^8)
    T = eltype(solver.u)
    dim = getdim(problem)
    dh = solver.problem.ch.dh
    N = getncells(dh.grid)
    σ = zeros(T, N)
    transf_matrices = Matrix{T}[]
    u_fn = Displacement(solver; maxfevals)
    R = zeros(T, (2, 2 * dim))
    for (cellidx, cell) in enumerate(CellIterator(dh))
        u, v = cell.coords[1], cell.coords[2]
        # R ∈ 2 x (2*dim)
        R_coord = compute_local_axes(u, v)
        fill!(R, 0.0)
        R[1, 1:dim] = R_coord[:, 1]
        R[2, (dim + 1):(2 * dim)] = R_coord[:, 2]
        push!(transf_matrices, R)
    end
    return TrussStress(σ, u_fn, transf_matrices, 0, maxfevals)
end

"""
# Arguments
`x` = design variables

# Returns
displacement vector `σ`, compressive stress < 0, tensile stress > 0
"""
function (ts::TrussStress{T})(x) where {T}
    @unpack σ, transf_matrices, u_fn = ts
    @unpack global_dofs, solver = u_fn
    @unpack penalty, problem, xmin = solver
    dh = getdh(problem)
    ts.fevals += 1
    u = u_fn(x)
    As = getA(problem)
    @unpack Kes = solver.elementinfo
    for e in 1:length(x)
        # Ke = R' * K_local * R
        # F = R * (R' * K_local * R) * u
        celldofs!(global_dofs, dh, e)
        σ[e] = -(transf_matrices[e] * Kes[e] * u[global_dofs])[1] / As[e]
    end
    return copy(σ)
end

# TODO complete
"""
rrule for autodiff.

du/dxe = -K^-1 * dK/dxe * u
d(u)/d(x_e) = - K^-1 * d(K)/d(x_e) * u
            = - K^-1 * (Σ_ei d(ρ_ei)/d(x_e) * K_ei) * u
            = - K^-1 * [d(ρ_e)/d(x_e) * K_e * u]
d(u)/d(x_e)' * Δ = -d(ρ_e)/d(x_e) * u' * K_e * (K^-1 * Δ)
"""
# function ChainRulesCore.rrule(dp::TrussStress, x)
#     @unpack dudx_tmp, solver, global_dofs = dp
#     @unpack penalty, problem, u, xmin = solver
#     dh = getdh(problem)
#     @unpack Kes = solver.elementinfo
#     # Forward-pass
#     # Cholesky factorisation
#     u = dp(x)
#     return u, Δ -> begin # v
#         solver.rhs .= Δ
#         solver(reuse_chol = true, assemble_f = false)
#         dudx_tmp .= 0
#         for e in 1:length(x)
#             _, dρe = get_ρ_dρ(x[e], penalty, xmin)
#             celldofs!(global_dofs, dh, e)
#             Keu = bcmatrix(Kes[e]) * u[global_dofs]
#             dudx_tmp[e] = -dρe * dot(Keu, solver.lhs[global_dofs])
#         end
#         return nothing, dudx_tmp # J1' * v, J2' * v
#     end
# end


text:
module MaxHelpingHandRainbowSpinnerColorFadeTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/RainbowSpinnerColorFadeTrigger" RainbowSpinnerColorFadeTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    colorsA::String="89E5AE,88E0E0,87A9DD,9887DB,D088E2", gradientSizeA::Number=280.0, loopColorsA::Bool=false, centerXA::Number=0.0, centerYA::Number=0.0, gradientSpeedA::Number=50.0,
    colorsB::String="89E5AE,88E0E0,87A9DD,9887DB,D088E2", gradientSizeB::Number=280.0, loopColorsB::Bool=false, centerXB::Number=0.0, centerYB::Number=0.0, gradientSpeedB::Number=50.0,
    direction::String="NoEffect", persistent::Bool=false)

const placements = Ahorn.PlacementDict(
    "Rainbow Spinner Color Fade (max480's Helping Hand)" => Ahorn.EntityPlacement(
        RainbowSpinnerColorFadeTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::RainbowSpinnerColorFadeTrigger)
    return Dict{String, Any}(
        "direction" => Maple.trigger_position_modes
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# A set of inplace broadcasting functions

# broadcast along specific dimension(s)

function _vbroadcast_eachcol!(m::Int, n::Int, f::Functor{2}, r::ContiguousArray, a::ContiguousArray, b::ContiguousArray)
    o = 0
    for j = 1 : n
        for i = 1 : m
            @inbounds r[o+i] = evaluate(f, a[o+i], b[i])
        end
        o += m
    end
end

function _vbroadcast_eachrow!(m::Int, n::Int, f::Functor{2}, r::ContiguousArray, a::ContiguousArray, b::ContiguousArray)
    o = 0
    for j = 1 : n
        @inbounds bj = b[j]
        for i = 1 : m
            @inbounds r[o+i] = evaluate(f, a[o+i], bj)
        end
        o += m
    end
end

function _vbroadcast3!{T}(m::Int, n::Int, k::Int,
    f::Functor{2}, r::ContiguousArray, a::ContiguousArray{T,3}, b::ContiguousArray)

    for l = 1 : k
        _vbroadcast_eachrow!(m, n, f, view(r,:,:,l), view(a,:,:,l), b)
    end
end

function vbroadcast!(f::Functor{2}, r::ContiguousArray, a::ContiguousArray, b::ContiguousArray, dim::Int)
    shp = size(a)
    nd = ndims(a)
    size(r) == shp && size(a, dim) == length(b) || error("Inconsistent argument dimensions.")
    1 <= dim <= nd || error("Invalid value of dim.")

    if dim == 1
        m = shp[1]
        n = succ_length(shp, 1)
        _vbroadcast_eachcol!(m, n, f, r, a, b)
    else
        m = prec_length(shp, dim)
        n = shp[dim]
        k = succ_length(shp, dim)
        if k == 1
            _vbroadcast_eachrow!(m, n, f, r, a, b)
        else
            _vbroadcast3!(m, n, k, f, ContiguousView(r, (m,n,k)), a, b)
        end
    end
    return r
end

vbroadcast1!(f::Functor{2}, a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast!(f, a, a, b, dim)

function vbroadcast(f::Functor{2}, a::ContiguousArray, b::ContiguousArray, dim::Int)
    R = result_type(f, eltype(a), eltype(b))
    vbroadcast!(f, zeros(R, size(a)), a, b, dim)
end

# Specific broadcasting function

badd!(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast1!(Add(), a, b, dim)
bsubtract!(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast1!(Subtract(), a, b, dim)
bmultiply!(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast1!(Multiply(), a, b, dim)
bdivide!(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast1!(Divide(), a, b, dim)

badd(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast(Add(), a, b, dim)
bsubtract(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast(Subtract(), a, b, dim)
bmultiply(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast(Multiply(), a, b, dim)
bdivide(a::ContiguousArray, b::ContiguousArray, dim::Int) = vbroadcast(Divide(), a, b, dim)


text:
function write_log_dat_ampl(lambda,M,N)
	theta = Vector{Float64}(N)
	rand!(theta)
	theta = theta.*7

	f = open(string("log",M,"_",N,".dat"),"w")

	write(f, string("# true theta \n"))
	for i=1:N
		write(f,@sprintf("# \t %0.2f \n", theta[i]))
	end

	write(f, "reset data; \n")
	write(f, "data; \n")
	write(f, string("param m := ",M,"; \n"))
	write(f, string("param n := ",N,"; \n"))
	write(f, string("param lambda := ", lambda ,"; \n"))
	write(f, "param \t x: ")
	for i=1:N
		write(f, string("\t ",i))
	end
	write(f, "\t := \n")

	x = randn(M,N)
	x[:,N] = ones(M)
	for i=1:M
		write(f, string("\t",i))
		xx = x[i,:]
		for j=1:N
			write(f, string("\t",xx[j]))
		end
		write(f, "\n");
	end
	write(f, "\t ; \n")
	write(f, "param \t y \t := \n")

	p = 1./(1+exp(-x*theta))
	y = Vector{Float64}(M)
	for i=1:M
		write(f, string("\t",i))
		y[i] = generate_bernoulli(p[i])
		write(f, string("\t", y[i], "\n"))
	end
	write(f, "\t ; \n")

	write(f, "option solver ipopt; \n");
	write(f, "solve; \n");

	close(f)
	return theta, x, y
end

function generate_bernoulli(p)
	a = rand()
	if(a<=p)
		return 1.0
	else
		return -1.0
	end
end

function write_log_dat_julia(lambda,t,x,y)
	M, N = size(x)
	f = open(string("log",M,"_",N,"_dat.jl"),"w")
	write(f, string("lambda = ",lambda,"\n"))
	write(f, string("M=",M,"\n"))
	write(f, string("N=",N,"\n"))
	write(f,string("x=Array{Float64,2}(M,N) \n"))
	for i = 1:M
		write(f, string("x[",i,",:]=",x[i,:], "\n"))
	end
	write(f,string("y=",y))
	write(f,"\n")

	close(f)
end

function genLogData(M,N)
	lambda = 1.0
	t, x, y = write_log_dat_ampl(lambda, M, N )
	write_log_dat_julia(lambda, t , x , y)
end

# ---------------------------
function genRandData(N,K)
	A = write_random_dat_ampl(N,K)
	writedlm(open(string("randset",N,"_",K,".txt"),"w"),A)
end

function rand_set(n,k)
    ret = []
    for i=1:k
        push!(ret,rand(1:n))
    end
    return ret
end

function genArrowHead(N,K)
	f = open(string("arrowhead",N,"_",K,".dat"),"w")
	write(f, "reset data; \n")
	write(f, string("data; \n"))
	write(f, string("param n:=",N,"; \n"))
	write(f, string("param k:=",K,"; \n"))
	write(f,string("option solver ipopt; \n"))
	write(f,string("solve; \n"))
end

function write_random_dat_ampl(N, K)
	f = open(string("randset",N,"_",K,".dat"),"w")

	write(f, string("data; \n"))
	write(f, string("param n:=",N,"; \n"))
	write(f, string("param k:=",K,"; \n"))
	write(f, string("param rand: "));
	for i =1:K
		write(f, string("\t ",i))
	end
	write(f, string(":= \n"))

	A = Array{Int, 2}(N,K);
	for i =1:N
		A[i,:] = rand_set(N,K)
	end

	for i=1:N
		AA = A[i,:]
		write(f, string("\t",i,"\t"))
		for aa in AA
			write(f, string("\t",aa))
		end
		write(f,"\n")
	end
	write(f,string("; \n"))
	write(f,string("option solver ipopt; \n"))
	write(f,string("solve; \n"))

	write(f, "\n")
	close(f)
	return A
# data;
# param n:=4000;
# param k:=2;

# #display N;
# #display K;

# for {i in 1..n}
# {
#     let {b in 1..k} rand[i, b] := round(Uniform(1,n));
# }

# #display rand;

# option solver ipopt;
# solve;

end


text:
using CImGui
using CImGui.CSyntax
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using CImGuiExtensions
using DataFrames # For this example, don't forget to add this package.
using Printf

struct UnspecifiedSchema <: AbstractSchema end

function create_import_context(keystr::String)
    # The controller is responsible for instantiating the dialog GUI,
    # handling user input and tiggering the file importer once the user
    # confirms the selection of a particular file.
    control = FileDialogControl(true)
    # The model store the path that the user has confirmed, as well as the
    # path that the user has selected but not necessarily confirmed.
    model = FileDialogModel(Path(pwd(),""), Path(pwd(),""))
    # The display properties include the title of the file dialog and what
    # action string is to be displayed on the button.
    properties = FileDialogDisplayProperties(caption = "Open File###"*keystr, action = "Open###"*keystr, width = Cfloat(520), height = Cfloat(210))
    # The importer is triggered once the user has confirmed the selection of a
    # particular path.
    load_csv = CSVImporter(false, UnspecifiedSchema())
    # Construct a context and assign roles to complete the requisite import interactions.
    ImportContext(control, model, properties, load_csv)
end

function open_file_example()

    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1280, 720, "Demo")
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    CImGui.StyleColorsDark()

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)

    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    data::Union{Nothing, DataFrame} = nothing
    # The string will be used to distingush between different open file dialogs
    # in case you ever need to have more than one file dialog active on the
    # screen.
    import_csv_key = "Import Generic CSV File"
    import_csv = create_import_context(import_csv_key)

    while !GLFW.WindowShouldClose(window)

        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        CImGui.Begin("Demo: Open CSV File")
        if CImGui.Button("Open File Dialog")
            enable!(import_csv)
        end
        data₀::Union{Nothing, DataFrame}  = import_csv()
        # isrunning(import_csv) == true while the user is browsing for
        # the file. Once they click "Open" or "Cancel" isrunning(import_csv)
        # will return false.
        if !isrunning(import_csv) && !isnothing(data₀)
            data = data₀
            # You might trigger some other event or set a flag to true now
            # that the data has been loaded and is available for use.
            # Here we just display the data in the REPL.
            display(data)
            # Disable the import use-case once you are finished.
            disable!(import_csv)
        end
        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

    # cleanup
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end

# Run the demo
 open_file_example()


text:
# AbstractVector wrapper around PETSc Vec types
export Vec, comm, NullVec

 """
  Construct a high level Vec object from a low level C.Vec.
  The data field is used to protect things from GC.
  A finalizer is attached to deallocate the memory of the underlying C.Vec, unless
  `first_instance` is set to true.
  `assembled` indicates when values are set via `setindex!` and is reset by
   `AssemblyEnd`
   `verify_assembled` when true, calls to `isassembled` verify all processes
   have `assembled` = true, when false, only the local assembly state is
   checked.  This essentially makes the user responsible for assembling
  the vector before passing it into functions that will use it (like KSP
  solves, etc.).

"""
mutable struct Vec{T} <: AbstractVector{T}
  p::C.Vec{T}
  assembled::Bool # whether are all values have been assembled
  verify_assembled::Bool # check whether all processes are assembled
  insertmode::C.InsertMode # current mode for setindex!
  data::Any # keep a reference to anything needed for the Mat
            # -- needed if the Mat is a wrapper around a Julia object,
            #    to prevent the object from being garbage collected.
  function Vec{T}(p::C.Vec{T}, data=nothing; first_instance::Bool=true,
               verify_assembled::Bool=true) where {T}
    v = new{T}(p, false, verify_assembled, C.INSERT_VALUES, data)
    if first_instance
      if (!deactivate_finalizers)
        finalizer(PetscDestroy,v)
      end
    end
    return v
  end
end

function Base.show(io::IO, x::Vec)
  myrank = MPI.Comm_rank(comm(x))
  if myrank == 0
    println("Petsc Vec of length ", length(x))
  end
  if isassembled(x)
    println(io, "Process ", myrank, " entries:")
    x_arr = LocalVector_readonly(x)
    show(io, x_arr)
    restore(x_arr)
  else
    println(io, "Process ", myrank, " not assembled")
  end
end

"""
  Null vectors, used in place of void pointers in the C
  API
"""
global const NullVec = Dict{DataType, Vec}()

if have_petsc[1]
  global const NullVec1 = Vec{Float64}(C.Vec{Float64}(C_NULL), first_instance=false)
  NullVec[Float64] = NullVec1
end
if have_petsc[2]
  global const NullVec2 = Vec{Float32}(C.Vec{Float32}(C_NULL), first_instance=false)
  NullVec[Float32] = NullVec2
end
if have_petsc[3]
  global const NullVec3 = Vec{Complex128}(C.Vec{Complex128}(C_NULL), first_instance=false)
  NullVec[Complex128] = NullVec3

end
 """
  Gets the MPI communicator of a vector.
"""
function comm(v::Vec{T}) where {T}
  return C.PetscObjectComm(T, v.p.pobj)
end

export gettype

 """
  Get the Symbol that is the format of the vector
"""
function gettype(a::Vec{T}) where {T}
  sym_arr = Array{C.VecType}(undef,1)
  chk(C.VecGetType(a.p, sym_arr))
  return sym_arr[1]
end

 """
  Create an empty, unsized vector.
"""
function Vec(::Type{T}, vtype::C.VecType=C.VECMPI;
                comm::MPI.Comm=MPI.COMM_WORLD) where {T}
  p = Ref{C.Vec{T}}()
  chk(C.VecCreate(comm, p))
  chk(C.VecSetType(p[], vtype))
  v = Vec{T}(p[])
  v
end

 """
  Create a vector, specifying the (global) length len or the local length
  mlocal.  Even if the blocksize is > 1, teh lengths are always number of
  elements in the vector, not number of block elements.  Thus
  len % blocksize must = 0.
"""
function Vec(::Type{T}, len::Integer=C.PETSC_DECIDE;
                         vtype::C.VecType=C.VECMPI,  bs=1,
                         comm::MPI.Comm=MPI.COMM_WORLD,
                         mlocal::Integer=C.PETSC_DECIDE) where {T<:Scalar}
  vec = Vec(T, vtype; comm=comm)
  resize!(vec, len, mlocal=mlocal)
  set_block_size(vec, bs)
  vec
end

 """
  Make a PETSc vector out of an array.  If used in parallel, the array becomes
  the local part of the PETSc vector
"""
# make a Vec that is a wrapper around v, where v stores the local data
function Vec(v::Vector{T}; comm::MPI.Comm=MPI.COMM_WORLD) where {T<:Scalar}
  p = Ref{C.Vec{T}}()
  chk(C.VecCreateMPIWithArray(comm, PetscInt(1), PetscInt(length(v)), PetscInt(C.PETSC_DECIDE), v, p))
  pv = Vec{T}(p[], v)
  return pv
end

function set_block_size(v::Vec{T}, bs::Integer) where {T<:Scalar}
  chk(C.VecSetBlockSize(v.p, PetscInt(bs)))
end

function get_blocksize(v::Vec{T}) where {T<:Scalar}
  bs = Ref{PetscInt}()
  chk(C.VecGetBlockSize(v.p, bs))
  return Int(bs[])
end

export VecGhost, VecLocal, restore

 """
  Make a PETSc vector with space for ghost values.  ghost_idx are the
  global indices that will be copied into the ghost space.
"""
# making mlocal the position and mglobal the keyword argument is inconsistent
# with the other Vec constructors, but it makes more sense here
function VecGhost(::Type{T}, mlocal::Integer,
                  ghost_idx::Array{I,1}; comm=MPI.COMM_WORLD, m=C.PETSC_DECIDE, bs=1, vtype=C.VECMPI) where {T<:Scalar, I <: Integer}

    nghost = length(ghost_idx)
    ghost_idx2 = [ PetscInt(i -1) for i in ghost_idx]

    vref = Ref{C.Vec{T}}()
    if bs == 1
      chk(C.VecCreateGhost(comm, PetscInt(mlocal), PetscInt(m), PetscInt(nghost), ghost_idx2, vref))
    elseif bs > 1
      chk(C.VecCreateGhostBlock(comm, PetscInt(bs), PetscInt(mlocal), PetscInt(mlocal), PetscInt(m), PetscInt(nghost), ghost_idx2, vref))
    else
      println(stderr, "WARNING: unsupported block size requested, bs = ", bs)
    end

    chk(C.VecSetType(vref[], vtype))

    return Vec{T}(vref[])
end

 """
  Create a VECSEQ that contains both the local and the ghost values of the
  original vector.  The underlying memory for the orignal and output vectors
  alias.
"""
function VecLocal( v::Vec{T}) where {T<:Scalar}

  vref = Ref{C.Vec{T}}()
  chk(C.VecGhostGetLocalForm(v.p, vref))
  # store v to use with Get/Restore LocalForm
  # Petsc reference counting solves the gc problem
  return Vec{T}(vref[], v)
end

#TODO: use restore for all types of restoring a local view
 """
  Tell Petsc the VecLocal is no longer needed
"""
function restore(v::Vec{T}) where {T}

  vp = v.data
  vref = Ref(v.p)
  chk(C.VecGhostRestoreLocalForm(vp.p, vref))
end

 """
  The Petsc function to deallocate Vec objects
"""
function PetscDestroy(vec::Vec{T}) where {T}
  if !PetscFinalized(T)  && !isfinalized(vec)
    C.VecDestroy(Ref(vec.p))
    vec.p = C.Vec{T}(C_NULL)  # indicate the vector is finalized
  end
end

 """
  Determine whether a vector has already been finalized
"""
function isfinalized(vec::Vec)
  return isfinalized(vec.p)
end

function isfinalized(vec::C.Vec)
  return vec.pobj == C_NULL
end

global const is_nullvec = isfinalized  # another name for doing the same check

 """
  Use the PETSc routine for printing a vector to stdout
"""
function petscview(vec::Vec{T}) where {T}
  viewer = C.PetscViewer{T}(C_NULL)
  chk(C.VecView(vec.p, viewer))
end

function Base.resize!(x::Vec, m::Integer=C.PETSC_DECIDE; mlocal::Integer=C.PETSC_DECIDE)
  if m == mlocal == C.PETSC_DECIDE
    throw(ArgumentError("either the length (m) or local length (mlocal) must be specified"))
  end
  chk(C.VecSetSizes(x.p, PetscInt(mlocal), PetscInt(m)))
  x
end

###############################################################################
export ghost_begin!, ghost_end!, scatter!, ghost_update!
# ghost vectors: essential methods
 """
  Start communication to update the ghost values (on other processes) from the local
  values
"""
function ghost_begin!(v::Vec{T}; imode=C.INSERT_VALUES,
                               smode=C.SCATTER_FORWARD) where {T<:Scalar}
    chk(C.VecGhostUpdateBegin(v.p, imode, smode))
    return v
end

 """
  Finish communication for updating ghost values
"""
function ghost_end!(v::Vec{T}; imode=C.INSERT_VALUES,
                               smode=C.SCATTER_FORWARD) where {T<:Scalar}
    chk(C.VecGhostUpdateEnd(v.p, imode, smode))
    return v
end

# ghost vectors: helpful methods
 """
  Convenience method for calling both ghost_begin! and ghost_end!
"""
function scatter!(v::Vec{T}; imode=C.INSERT_VALUES, smode=C.SCATTER_FORWARD) where {T<:Scalar}

  ghost_begin!(v, imode=imode, smode=smode)
  ghost_end!(v, imode=imode, smode=smode)
end

# is there a way to specify all varargs must be same type?
# this can't be named scatter! because of ambiguity with the index set scatter!
 """
  Convenience method for calling ghost_begin! and ghost_end! for multiple vectors
"""
function ghost_update!(v...; imode=C.INSERT_VALUES, smode=C.SCATTER_FORWARD)

  for i in v
    ghost_begin!(i, imode=imode, smode=smode)
  end

  for i in v
    ghost_end!(i, imode=imode, smode=smode)
  end

  return v
end

###############################################################################
export lengthlocal, sizelocal, localpart

Base.convert(::Type{C.Vec}, v::Vec) = v.p

import Base.length
 """
  Get the global length of the vector
"""
function length(x::Vec)
  sz = Ref{PetscInt}()
  chk(C.VecGetSize(x.p, sz))
  Int(sz[])
end

 """
  Get the global size of the vector
"""
Base.size(x::Vec) = (length(x),)

 """
  Get the length of the local portion of the vector
"""
function lengthlocal(x::Vec)
  sz = Ref{PetscInt}()
  chk(C.VecGetLocalSize(x.p, sz))
  sz[]
end

"""
  Get the local size of the vector
"""
sizelocal(x::Vec) = (lengthlocal(x),)

"""
  Get local size of the vector
"""
sizelocal(t::AbstractArray{T,n}, d) where {T,n} = (d>n ? 1 : sizelocal(t)[d])

 """
  Get the range of global indices that define the local part of the vector.
  Internally, this calls the Petsc function VecGetOwnershipRange, and has
  the same limitations as that function, namely that some vector formats do
  not have a well defined contiguous range.
"""
function localpart(v::Vec)
  # this function returns a range from the first to the last indicies (1 based)
  # this is different than the Petsc VecGetOwnershipRange function where
  # the max value is one more than the number of entries
  low = Ref{PetscInt}()
  high = Ref{PetscInt}()
  chk(C.VecGetOwnershipRange(v.p, low, high))
  return (low[]+1):(high[])
end

"""
  Similar to localpart, but returns the range of block indices
"""
function localpart_block(v::Vec)
  low = Ref{PetscInt}()
  high = Ref{PetscInt}()
  chk(C.VecGetOwnershipRange(v.p, low, high))
  bs = get_blocksize(v)
  low_b = div(low[], bs); high_b = div(high[]-1, bs)
  ret = (low_b+1):(high_b+1)

  return ret
end

function Base.similar(x::Vec{T}) where {T}
  p = Ref{C.Vec{T}}()
  chk(C.VecDuplicate(x.p, p))
  v=Vec{T}(p[])
  v.assembled        = x.assembled
  v.verify_assembled = x.verify_assembled
  v.data             = x.data
  v.insertmode       = x.insertmode
  v
end

Base.similar(x::Vec{T}, ::Type{T}) where {T} = similar(x)
function Base.similar(x::Vec{T}, T2::Type) where {T}
  VType = gettype(x)
  Vec(T2, length(x), VType; comm=comm(x), mlocal=lengthlocal(x))
end

function Base.similar(x::Vec{T}, T2::Type, len::Union{Int,Dims{1}}) where {T}
  VType = gettype(x)
  len[1]==length(x) && T2==T ? similar(x) : Vec(T2, len[1], vtype=VType; comm=comm(x))
end

function Base.similar(x::Vec{T}, len::Union{Int,Dims{1}}) where {T}
  VType = gettype(x)
  len[1]==length(x) ? similar(x) : Vec(T, len[1], vtype=VType; comm=comm(x))
end

function Base.copy(x::Vec)
  AssemblyBegin(x)
  y = similar(x)
  AssemblyEnd(x)
  chk(C.VecCopy(x.p, y.p))
  y
end

function Base.copy!(dst::Vec,src::Vec)
  chk(C.VecCopy(src.p, dst.p))
end

###############################################################################
export localIS, local_to_global_mapping, set_local_to_global_mapping, has_local_to_global_mapping

"""
  Constructs index set mapping from local indexing to global indexing, based
  on localpart()
"""
function localIS(A::Vec{T}) where {T}

  rows = localpart(A)
  rowis = IS(T, rows, comm=comm(A))
  return rowis
end

"""
  Like localIS, but returns a block index IS
"""
function localIS_block(A::Vec{T}) where {T}
  rows = localpart_block(A)
  bs = get_blocksize(A)
  rowis = ISBlock(T, bs, rows, comm=comm(A))
#  set_blocksize(rowis, get_blocksize(A))
  return rowis
end
"""
  Construct ISLocalToGlobalMappings for the vector.  If a block vector,
  create a block index set
"""
function local_to_global_mapping(A::Vec)

  # localIS creates strided index sets, which require only constant
  # memory
  if get_blocksize(A) == 1
    rowis = localIS(A)
  else
    rowis = localIS_block(A)
  end
  row_ltog = ISLocalToGlobalMapping(rowis)

  return row_ltog
end

# need a better name
"""
  Registers the ISLocalToGlobalMapping with the Vec
"""
function set_local_to_global_mapping(A::Vec{T}, rmap::ISLocalToGlobalMapping{T}) where {T}

  chk(C.VecSetLocalToGlobalMapping(A.p, rmap.p))
end

"""
  Check if the local to global mapping has been registered
"""
function has_local_to_global_mapping(A::Vec{T}) where {T}

  rmap_ref = Ref{C.ISLocalToGlobalMapping{T}}()
  chk(C.VecGetLocalToGlobalMapping(A.p, rmap_re))

  rmap = rmap_ref[]

  return rmap.pobj != C_NULL
end

##########################################################################
import Base: setindex!
export assemble, isassembled, AssemblyBegin, AssemblyEnd

# for efficient vector assembly, put all calls to x[...] = ... inside
# assemble(x) do ... end
 """
  Start communication to assemble stashed values into the vector

  The MatAssemblyType is not needed for vectors, but is provided for
  compatibility with the Mat case.

  Unless vec.verify_assembled == false, users must *never* call the
  C functions VecAssemblyBegin, VecAssemblyEnd and VecSetValues, they must
  call AssemblyBegin, AssemblyEnd, and setindex!.
"""
function AssemblyBegin(x::Vec, t::C.MatAssemblyType=C.MAT_FINAL_ASSEMBLY)
  chk(C.VecAssemblyBegin(x.p))
end

"""
  Generic fallback for AbstractArray, no-op
"""
function AssemblyBegin(x::AbstractArray, t::C.MatAssemblyType=C.MAT_FINAL_ASSEMBLY)

end
 """
  Finish communication for assembling the vector
"""
function AssemblyEnd(x::Vec, t::C.MatAssemblyType=C.MAT_FINAL_ASSEMBLY)
  chk(C.VecAssemblyEnd(x.p))
  x.assembled = true
end

"""
  Check if a vector is assembled (ie. does not have stashed values).  If
  `x.verify_assembled`, the assembly state of all processes is checked,
  otherwise only the local process is checked. `local_only` forces only
  the local process to be checked, regardless of `x.verify_assembled`.
"""
function isassembled(x::Vec, local_only=false)
  myrank = MPI.Comm_rank(comm(x))
  if x.verify_assembled && !local_only
    val = MPI.Allreduce(Int8(x.assembled), MPI.LAND, comm(x))
  else
    val = x.assembled
  end

  return Bool(val)
end

"""
  Generic fallback for AbstractArray, no-op
"""
function AssemblyEnd(x::AbstractArray, t::C.MatAssemblyType=C.MAT_FINAL_ASSEMBLY)

end

isassemble(x::AbstractArray) = true
# assemble(f::Function, x::Vec) is defined in mat.jl

 """
  Like setindex, but requires the indices be 0-base
"""
function setindex0!(x::Vec{T}, v::Array{T}, i::Array{PetscInt}) where {T}
  n = length(v)
  if n != length(i)
    throw(ArgumentError("length(values) != length(indices)"))
  end
  #    println("  in setindex0, passed bounds check")
  chk(C.VecSetValues(x.p, PetscInt(n), i, v, x.insertmode))
  x.assembled = false
  x
end

function setindex!(x::Vec{T}, v::Number, i::Integer) where {T}
  # can't call VecSetValue since that is a static inline function
  setindex0!(x, T[ v ], PetscInt[ i - 1 ])
  v
end

# set multiple entries to a single value
setindex!(x::Vec, v::Number, I::AbstractArray{T}) where {T<:Integer} = assemble(x) do
  for i in I
    x[i] = v
  end
  x
end

function Base.fill!(x::Vec{T}, v::Number) where {T}
  chk(C.VecSet(x.p, T(v)))
  return x
end

function setindex!(x::Vec, v::Number, I::AbstractRange{T}) where {T<:Integer}
  if abs(step(I)) == 1 && minimum(I) == 1 && maximum(I) == length(x)
    fill!(x, v)
    return v
  else
    # use invoke here to avoid a recursion loop
    return invoke(setindex!, Tuple{Vec,typeof(v),AbstractVector{T}}, x,v,I)
  end
end

#TODO: make this a single call to VecSetValues
setindex!(x::Vec, V::AbstractArray, I::AbstractArray{T}) where {T<:Real} =
assemble(x) do
  if length(V) != length(I)
    throw(ArgumentError("length(values) != length(indices)"))
  end
  # possibly faster to make a PetscScalar array from V, and
  # a copy of the I array shifted by 1, to call setindex0! instead?
  c = 1
  for i in I
    x[i] = V[c]
    c += 1
  end
  x
end

# logical indexing
setindex!(A::Vec, x::Number, I::AbstractArray{Bool}) = assemble(A) do
  for i = 1:length(I)
    if I[i]
      A[i] = x
    end
  end
  A
end
for T in (:(Array{T2}),:(AbstractArray{T2})) # avoid method ambiguities
  @eval setindex!(A::Vec, X::$T, I::AbstractArray{Bool}) where {T2<:Scalar} = assemble(A) do
    c = 1
    for i = 1:length(I)
      if I[i]
        A[i] = X[c]
        c += 1
      end
    end
    A
  end
end

##########################################################################
import Base.getindex

# like getindex but for 0-based indices i
function getindex0(x::Vec{T}, i::Vector{PetscInt}) where {T}
  v = similar(i, T)
  chk(C.VecGetValues(x.p, PetscInt(length(i)), i, v))
  v
end

getindex(x::Vec, i::Integer) = getindex0(x, PetscInt[i-1])[1]

getindex(x::Vec, I::AbstractVector{PetscInt}) =
  getindex0(x, PetscInt[ (i-1) for i in I ])

##########################################################################
# more indexing
# 0-based (to avoid temporary copies)
export set_values!, set_values_blocked!, set_values_local!, set_values_blocked_local!

#TODO: in 0.5, use boundscheck macro to verify stride=1

function set_values!(x::Vec{T}, idxs::StridedVecOrMat{PetscInt},
                                 vals::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}

  chk(C.VecSetValues(x.p, PetscInt(length(idxs)), idxs, vals, o))
end

function set_values!(x::Vec{T}, idxs::StridedVecOrMat{I},
                                         vals::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar, I <: Integer}

  # convert idxs to PetscInt
  p_idxs = PetscInt[ i for i in idxs]
  set_values!(x, p_idxs, vals, o)
end

function set_values!(x::AbstractVector, idxs::AbstractArray, vals::AbstractArray,
                     o::C.InsertMode=C.INSERT_VALUES)

  if o == C.INSERT_VALUES
    for i=1:length(idxs)
      x[idxs[i] + 1] = vals[i]
    end
  elseif o == C.ADD_VALUES
    for i=1:length(idxs)
      x[idxs[i] + 1] += vals[i]
    end
  else
    throw(ArgumentError("Unsupported InsertMode"))
  end
end

function set_values_blocked!(x::Vec{T}, idxs::StridedVecOrMat{PetscInt},
                                          vals::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}

  chk(C.VecSetValuesBlocked(x.p, PetscInt(length(idxs)), idxs, vals, o))
end

function set_values_blocked!(x::Vec{T},
                             idxs::StridedVecOrMat{I}, vals::StridedVecOrMat{T},
                             o::C.InsertMode=x.insertmode) where {T <: Scalar, I <: Integer}

  p_idxs = PetscInt[ i for i in idxs]
  set_values_blocked!(x, p_idxs, vals, o)
end

# julia doesn't have blocked vectors, so skip

function set_values_local!(x::Vec{T}, idxs::StridedVecOrMat{PetscInt},
                                       vals::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}

  chk(C.VecSetValuesLocal(x.p, PetscInt(length(idxs)), idxs, vals, o))
end

function set_values_local!(x::Vec{T},
                           idxs::StridedVecOrMat{I}, vals::StridedVecOrMat{T},
                           o::C.InsertMode=x.insertmode) where {T <: Scalar, I <: Integer}

  p_idxs = PetscInt[ i for i in idxs]
  set_values_local!(x, p_idxs, vals, o)
end

# for julia vectors, local = global
function set_values_local!(x::AbstractArray, idxs::AbstractArray,
                           vals::AbstractArray, o::C.InsertMode=C.INSERT_VALUES)

  if o == C.INSERT_VALUES
    for i=1:length(idxs)
      x[idxs[i] + 1] = vals[i]
    end
  elseif o == C.ADD_VALUES
    for i=1:length(idxs)
      x[idxs[i] + 1] += vals[i]
    end
  else
    throw(ArgumentError("Unsupported InsertMode"))
  end

end

function set_values_blocked_local!(x::Vec{T},
                                   idxs::StridedVecOrMat{PetscInt},
                                   vals::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}

  chk(C.VecSetValuesBlockedLocal(x.p, PetscInt(length(idxs)), idxs, vals, o))
end

function set_values_blocked_local!(x::Vec{T},
                           idxs::StridedVecOrMat{I}, vals::StridedVecOrMat{T},
                           o::C.InsertMode=x.insertmode) where {T <: Scalar, I <: Integer}

  p_idxs = PetscInt[ i for i in idxs]
  set_values_blocked_local!(x, p_idxs, vals, o)
end

# julia doesn't have blocked vectors, so skip

###############################################################################
import Base: abs, exp, log, conj, conj!
export abs!, exp!, log!
for (f,pf) in ((:abs,:VecAbs), (:exp,:VecExp), (:log,:VecLog),
  (:conj,:VecConjugate))
  fb = Symbol(string(f, "!"))
  @eval begin
    function $fb(x::Vec)
      chk(C.$pf(x.p))
      x
    end
    $f(x::Vec) = $fb(copy(x))
  end
end

export chop!
function chop!(x::Vec, tol::Real)
  chk(C.VecChop(x.p, tol))
#  chk(ccall((:VecChop, petsc), PetscErrorCode, (pVec, PetscReal), x, tol))
  x
end

for (f, pf, sf) in ((:findmax, :VecMax, :maximum), (:findmin, :VecMin, :minimum))
  @eval begin
    function Base.$f(x::Vec{T}) where {T<:Real}
      i = Ref{PetscInt}()
      v = Ref{T}()
      chk(C.$pf(x.p, i, v))
      (v[], i[]+1)
    end
    Base.$sf(x::Vec{T}) where {T<:Real} = $f(x)[1]
  end
end
# For complex numbers, VecMax and VecMin apparently return the max/min
# real parts, which doesn't match Julia's maximum/minimum semantics.
function LinearAlgebra.norm(x::Union{Vec{T},Vec{Complex{T}}}, p::Real) where {T<:Real}
  v = Ref{T}()
  n = p == 1 ? C.NORM_1 : p == 2 ? C.NORM_2 : p == Inf ? C.NORM_INFINITY :
  throw(ArgumentError("unrecognized Petsc norm $p"))
  chk(C.VecNorm(x.p, n, v))
  v[]
end

#if VERSION >= v"0.5.0-dev+8353" # JuliaLang/julia#13681
#  import Base.normalize!
#else
#  export normalize!
#end

 """
  computes v = norm(x,2), divides x by v, and returns v
"""
#function normalize!(x::Union{Vec{T},Vec{Complex{T}}}) where {T<:Real}
#  v = Ref{T}()
#  chk(C.VecNormalize(x.p, v))
#  v[]
#end

function LinearAlgebra.dot(x::Vec{T}, y::Vec{T}) where {T}
  d = Ref{T}()
  chk(C.VecDot(y.p, x.p, d))
  return d[]
end

# unconjugated dot product (called for x'*y)
#function LinearAlgebra.At_mul_B(x::Vec{T}, y::Vec{T}) where {T<:Complex}
#  d = Array(1)
#  chk(C.VecTDot(x.p, y.p, d))
#  return d
#end

# pointwise operations on pairs of vectors (TODO: support in-place variants?)
import Base: broadcast
#for (f,pf) in ((:max,:VecPointwiseMax), (:min,:VecPointwiseMin),
#  (:.*,:VecPointwiseMult), (:./,:VecPointwiseDivide))
#  @eval function broadcast(::typeof($f), x::Vec, y::Vec)
#    w = similar(x)
#    chk(C.$pf(w.p, x.p, y.p))
#    w
#  end
#end

import Base: +, -, *, /, \
export scale!
function scale!(x::Vec{T}, s::Number) where {T}
  chk(C.VecScale(x.p, T(s)))
  x
end
scale(x::Vec{T},s::Number) where {T} = scale!(copy(x),s)
(*)(x::Vec, a::Number...) = scale(x, prod(a))
(*)(a::Number, x::Vec) = scale(x, a)
(/)(x::Vec, a::Number) = scale(x, inv(a))
(\)(a::Number, x::Vec) = scale(x, inv(a))
function Base.broadcast(::typeof(/), a::Number, x::Vec)
  y = copy(x)
  chk(C.VecReciprocal(y.p))
  if a != 1.0
    scale!(y, a)
  end
  y
end

function (+)(x::Vec{T}, a::Number...) where {T<:Scalar}
  y = copy(x)
  chk(C.VecShift(y.p, T(sum(a))))
  return y
end
(+)(a::Number, x::Vec{T}) where {T<:Scalar} = x + a
(-)(x::Vec{T}, a::Number) where {T<:Scalar} = x + (-a)
(-)(x::Vec) = scale(x, -1)
function (-)(a::Number, x::Vec{T}) where {T<:Scalar}
  y = -x
  chk(C.VecShift(y.p, T(a)))
  return y
end

import Base: ==
function (==)(x::Vec, y::Vec)
  b = Ref{PetscBool}()
  chk(C.VecEqual(x.p, y.p, b))
  b[] != 0
end

function (==)(x::Vec, y::AbstractArray)
  flag = true
  x_arr = LocalVector(x)
  for i=1:length(x_arr)  # do localpart, then MPI reduce
    flag = flag && x_arr[i] == y[i]
  end
  restore(x_arr)

  buf = Int8[flag]
  # process 0 is root
  recbuf = MPI.Reduce(buf, MPI.LAND, 0, comm(x))

  if  MPI.Comm_rank(comm(x)) == 0
    buf[1] = recbuf[1]
  end

  MPI.Bcast!(view(buf,1:1), 0, comm(x))

  return convert(Bool, buf[1])
end

function Base.sum(x::Vec{T}) where {T}
  s = Ref{T}()
  chk(C.VecSum(x.p, s))
  s[]
end

###############################################################################
# map and friends
import Base: map!, map
#map() should be inherited from base

function map!(f, c)
  map!(f, c, c)
end

"""
Applys f element-wise to src to populate dest.  If src is a ghost vector,
then f is applied to the ghost elements as well as the local elements.
"""
function map!(f::F, dest::Vec{T}, src::Vec) where {T,F}
  if length(dest) < length(src)
    throw(ArgumentError("Length of dest must be >= src"))
  end
  if localpart(dest)[1] != localpart(src)[1]
    throw(ArgumentError("start of local part of src and dest must be aligned"))
  end

  dest_arr = LocalVector(dest)
  src_arr = LocalVector_readonly(src)
  try
    for (idx, val) in enumerate(src)
      dest[idx] = f(val)
    end
  finally
    restore(dest_arr)
    restore(src_arr)
  end
end

"""
  Multiple source vector map.  All vectors must have the local and global
  lengths.  If some a ghost vectors and some are not, the map is applied
  only to the local part
"""
function map!(f::F, dest::Vec{T}, src1::Vec{T}, src2::Vec{T2},  src_rest::Vec{T2}...) where {F,T,T2}

  # annoying workaround for #13651
  srcs = (src1, src2, src_rest...)
  # check lengths
  dest_localrange = localpart(dest)
  dest_len = length(dest)
  for src in srcs
    srclen = length(src)
    srcrange_local = localpart(src)
    if dest_len < srclen
      throw(ArgumentError("Length of destination must be greater than source"))
    end

    if dest_localrange[1] != srcrange_local[1]
      throw(ArgumentError("start of local part of src and dest must be aligned"))
    end
  end

  # extract the arrays
  n = length(srcs)
  len = 0
  len_prev = 0
  src_arrs = Array{LocalVectorRead{T2},1}(undef,n)
  use_length_local = false

  dest_arr = LocalVector(dest)
  try
    for (idx, src) in enumerate(srcs)
      src_arrs[idx] = LocalVector_readonly(src)

      # check of length of arrays are same or not
      len = length(src_arrs[idx])
      if len != len_prev && idx != 1 && !use_length_local
        use_length_local = true
      end
      len_prev = len
    end

    # if not all same, do only the local part (which must be the same for all)
    if use_length_local
      min_length = lenth(src1)
    else
      min_length = length(src_arrs[1])
    end
      # do the map
      vals = Array{T,1}(undef,n)
      for i=1:min_length  # TODO: make this the minimum array length
        for j=1:n  # extract values
          vals[j] = src_arrs[j][i]
        end
        dest_arr[i] = f(vals...)
      end
  finally # restore the arrays
    for src_arr in src_arrs
      restore(src_arr)
    end
    restore(dest_arr)
  end
end

##########################################################################
export axpy!, aypx!, axpby!, axpbypcz!
import LinearAlgebra.BLAS.axpy!
import LinearAlgebra.BLAS.axpby!

# y <- alpha*x + y
function axpy!(alpha::Number, x::Vec{T}, y::Vec{T}) where {T}
  chk(C.VecAXPY(y.p, T(alpha), x.p))
  y
end
# w <- alpha*x + y
function axpy!(alpha::Number, x::Vec{T}, y::Vec{T}, w::Vec{T}) where {T}
  chk(C.VecWAXPY(w.p, T(alpha), x.p, y.p))
  y
end
# y <- alpha*y + x
function aypx!(x::Vec{T}, alpha::Number, y::Vec{T}) where {T}
  chk(C.VecAYPX( y.p, T(alpha), x.p))
  y
end
# y <- alpha*x + beta*y
function axpby!(alpha::Number, x::Vec{T}, beta::Number, y::Vec{T}) where {T}
  chk(C.VecAXPBY(y.p, T(alpha), T(beta), x.p))
  y
end
# z <- alpha*x + beta*y + gamma*z
function axpbypcz!(alpha::Number, x::Vec{T}, beta::Number, y::Vec{T},
  gamma::Number, z::Vec{T}) where {T}
  chk(C.VecAXPBYPCZ(z.p, T(alpha), T(beta), T(gamma), x.p, y.p))
  z
end

# y <- y + \sum_i alpha[i] * x[i]
function axpy!(y::V, alpha::AbstractArray, x::AbstractArray{V}) where {V<:Vec}
  n = length(x)
  length(alpha) == n || throw(BoundsError())
  _x = [X.p for X in x]
  _alpha = eltype(y)[a for a in alpha]
  C.VecMAXPY(y.p, PetscInt(n), _alpha, _x)
  y
end

##########################################################################
# element-wise vector operations:
# import Base: .*, ./, .^, +, -

#for (f,pf) in ((:.*,:VecPointwiseMult), (:./,:VecPointwiseDivide), (:.^,:VecPow))
#  @eval function ($f)(x::Vec, y::Vec)
#    z = similar(x)
#    chk(C.$pf(z.p, x.p, y.p))
#    return z
#  end
#end

for (f,s) in ((:+,1), (:-,-1))
  @eval function ($f)(x::Vec{T}, y::Vec{T}) where {T}
    z = similar(x)
    chk(C.VecWAXPY(z.p, T($s), y.p, x.p))
    return z
  end
end

##############################################################################
export LocalVector, LocalVector_readonly, restore

 """
  Object representing the local part of the array, accessing the memory directly.
  Supports all the same indexing as a regular Array
"""
mutable struct LocalVector{T <: Scalar, ReadOnly} <: DenseArray{T, 1}
  a::Array{T, 1}  # the array object constructed around the pointer
  ref::Ref{Ptr{T}}  # reference to the pointer to the data
  pobj::C.Vec{T}
  isfinalized::Bool  # has this been finalized yet
  function LocalVector{T,ReadOnly}(a::Array{T}, ref::Ref, ptr) where {T,ReadOnly}
    varr = new{T,ReadOnly}(a, ref, ptr, false)
    # backup finalizer, shouldn't ever be used because users must call
    # restore before their changes will take effect
    if (!deactivate_finalizers)
      finalizer(restore, varr)
    end
    return varr
  end

end

const LocalVectorRead{T}=LocalVector{T, true}
const LocalVectorWrite{T}=LocalVector{T, false}

"""
  Get the LocalVector of a vector.  Users must call restore when
  finished updating the vector
"""
function LocalVector(vec::Vec{T}) where {T}

  len = lengthlocal(vec)

  ref = Ref{Ptr{T}}()
  chk(C.VecGetArray(vec.p, ref))
  a = unsafe_wrap(Array, ref[], len)
  return LocalVector{T, false}(a, ref, vec.p)
end

"""
  Get the LocalVector of a vector.  Users must call restore when
  finished updating the vector
"""
function LocalVector(vec::Vec{T},len::Int) where {T}
  ref = Ref{Ptr{T}}()
  chk(C.VecGetArray(vec.p, ref))
  a = unsafe_wrap(Array, ref[], len)
  return LocalVector{T, false}(a, ref, vec.p)
end

"""
  Tell Petsc the LocalVector is no longer being used
"""
function restore(varr::LocalVectorWrite{T}) where {T}

  if !varr.isfinalized && !PetscFinalized(T) && !isfinalized(varr.pobj)
    ptr = varr.ref
    chk(C.VecRestoreArray(varr.pobj, ptr))
  end
  varr.isfinalized = true
end

"""
  Get the LocalVector_readonly of a vector.  Users must call restore when
  finished with the object.
"""
function LocalVector_readonly(vec::Vec{T}) where {T}

  len = lengthlocal(vec)

  ref = Ref{Ptr{T}}()
  chk(C.VecGetArrayRead(vec.p, ref))
  a = unsafe_wrap(Array, ref[], len)
  return LocalVector{T, true}(a, ref, vec.p)
end

function restore(varr::LocalVectorRead{T}) where {T}

  if !varr.isfinalized && !PetscFinalized(T) && !isfinalized(varr.pobj)
    ptr = [varr.ref[]]
    chk(C.VecRestoreArrayRead(varr.pobj, ptr))
  end
  varr.isfinalized = true
end

Base.size(varr::LocalVector) = size(varr.a)
# indexing
getindex(varr::LocalVector, i) = getindex(varr.a, i)
setindex!(varr::LocalVectorWrite, v, i) = setindex!(varr.a, v, i)
Base.unsafe_convert(::Type{Ptr{T}}, a::LocalVector{T}) where {T} = Base.unsafe_convert(Ptr{T}, a.a)
Base.stride(a::LocalVector, d::Int64) = stride(a.a, d)
Base.similar(a::LocalVector, T::Type=eltype(a), dims::Dims{1}=size(a)) = similar(a.a, T, dims)

function (==)(x::LocalVector, y::AbstractArray)
  return x.a == y
end


text:
#=
  Lagrange Interpolation
  Newton forward difference interpolation
  Trapezoidal Quadrature
  Simpsons Quadrature
  Composite Trapezoidal Rule Quadrature
  Composite Simpsons Rule Quadrature
=#

function Lagrange(n,x,y)
  #accepts a 2 by n matrix
  #evalates at x
  for i = [2,1]:[2,n]
    a = 1
    k = i
    for j = 1:
      if k == j
        continue
      else
        a *= ((x-x[j]) / (x[k]-x[j]))
      end
    end
    f[i] = a * y[i]
    for i = 1:n-1
      sum += f[i]
    end
  end
  return sum

function Newton(n,x,y)
  #accepts x,y vectors
  for j = 1:n-1
    for i = n-1:-1:j
      y[i] = (y[i]-y[i-1]) / (x[i]-x[i-j-1])
    end
    for i = n-1:-1:0
      a = 1;
      for j= 1:i-1
          a *= a - x[j]
      end
      a *= y[j]
      sum += a
    end
  end
  return sum

function trapezoidalRule(p,a,b)
  #equally spaced nodes
  return h * (y[1]/2 + l1(y[2:end-1]) + y[end]/2)

function simpsonsRule(p,a,b)
  #even number eqaually spaced nodes
  a = 0
  b = 0
  for i = 2:2:end-1
    a += 2 * y[i]
  end
  for j = 3:2:end-1
    b += 4 * y[i]
  end
  return h/3 * (y[1] + a + b + y[end])

function compositeTrapezoidal(p,h,a,b)
  #equally-spaced nodes
  return trapezoidalRule(p,a,b) + (h * l1(y[1:end-1])

function compositeSimpsons(p,n,h,a,b)
  #even number n equally spaced nodes
  a = 0
  for i = 1:n/2
    a += y[2i - 2] + 4y[2i-1] + y[2i])
  end
  return h/3 * a

max_stars_repo_name: K-P-S/Julia-Scientific-Computing

max_issues_repo_name: K-P-S/Julia-Scientific-Computing

max_forks_repo_name: K-P-S/Julia-Scientific-Computing


text:
  #ANNOTAZIONI:
#-Lettere greche: angoli in rad, lettere latine: angoli in gradi
#------------------------------------
using Stripeline
using AstroLib, Dates
using DelimitedFiles
using Printf
using Plots
using Healpix
using FITSIO
using DataFrames
using Rotations
using CSV

output_dir=raw"C:\Users\Leo\OneDrive - studenti.unimi.it\OneDrive - Università degli Studi di Milano\UNI\Tesi"
RPM = 1  #velocità motore zenitale
Ang_Lim = 10
Ang_Zenith = 20
start_day = DateTime( 2022, 1, 1, 0, 0)
telescope_motors(time_s) = (0.0, deg2rad(Ang_Zenith), timetorotang(time_s, RPM))
#Creo le struct dei vari corpi che analizzo
struct SUN end
struct Moon end
#DICTIONARY dei pianeti!
Planets=Dict("Mercury"=>1, "Venus"=>2, "Mars"=>4, "Jupiter"=>5, "Saturn"=>6, "Uranus"=>7, "Neptune"=>8, "Pluto"=>9)

function Horn_dir(H_ID::String, date, time_range_s) #Dove punta l'Horn H_ID a partire da data "date"
    db = InstrumentDB()
    H_o=db.focalplane[H_ID].orientation
    dirs, orientations = genpointings(
        telescope_motors,
        H_o,
        time_range_s,
        date)
    return dirs
end

function Horn_dir(H_o::Array, date, time_range_s)
    dirs, orientations = genpointings(
        telescope_motors,
        H_o,
        time_range_s,
        date)
    return dirs
end

function scalar_prod(v, w)
    aux = sum(v.*w)
end

function mod(v)
    aux = sqrt(sum(v.^2))
end

function ang_dist(α1, δ1, α2, δ2) #Distanza angolare
    θ = acos(sin(δ1)*sin(δ2) + cos(δ1)*cos(δ2)*cos(α1 - α2))
    return θ
end

function ang_dist(v, w) #Distanza angolare tra due vettori
    θ = acos(scalar_prod(v, w)/(mod(v)*mod(w)))
    return θ
end

function Pos(T::Type{SUN}, date)
    α, δ,_,_= AstroLib.sunpos(jdcnv(date), radians=true) #Da valori in rad
    return α, δ
end

function Pos(T::Type{Moon}, date)
    α, δ,_,_,_= AstroLib.moonpos(jdcnv(date), radians=true) #Da valori in rad
    return α, δ
end

function Pos(Pl_ID::Int64, date)
    ra, dec= AstroLib.planet_coords(jdcnv(date), Pl_ID) #Da valori in gradi
    α, δ = (deg2rad(ra), deg2rad(dec))
    return α, δ
end

function Pos(Pl_Name::String, date)
    Pl_ID = Planets[Pl_Name]
    ra, dec= AstroLib.planet_coords(jdcnv(date), Pl_ID) #Da valori in gradi
    α, δ = (deg2rad(ra), deg2rad(dec))
    return α, δ
end

function Dist(Obj, H_dir, date) #Restituisce posizione di un corpo celeste Obj e distanza
    α2, δ2 = Pos(Obj, date)     #dall'Horn con direzione H_dir in data date
    δ1, α1 = H_dir[1], H_dir[2] #dirs sono Dec e AR in rad.
    θ = ang_dist(α1, δ1, α2, δ2)
    return θ, α2, δ2
end

function Map_Plot_Horn(H_ID, start_day, sampling_time_s, lasting_time_s, map, save::Bool)
    filename=joinpath(output_dir,"$H_ID.map")
    time0=datetime2unix(start_day)
    time_range = 0:sampling_time_s:lasting_time_s
    db = InstrumentDB()
    H_o=db.focalplane[H_ID].orientation
    H_dirs=Horn_dir(H_o, start_day, time_range)
    for indx in 1:length(collect(time_range))
        pixel_index = ang2pix(map, π/2 - H_dirs[indx, 1], H_dirs[indx, 2])
        map[pixel_index] += 1 #aumento hit count del telescopio
    end
    if save==true
        saveToFITS(map, "!"*filename)
        return filename
    end
end

function Create_Map(Horns, lasting_time_s, name)
    #filename=joinpath(output_dir,"$name.map")
    map = Healpix.Map{Float64, RingOrder}(128)
    for indx in 1:length(Horns)
        Map_Plot_Horn(Horns[indx], DateTime( 2022, 1, 1, 0, 0), 0.1, lasting_time_s, map, false)
    end
    plot(map, mollweide, title="Map $name")
    saveToFITS(map, "!$name.map")
    return map, "$name.map"
end

function Angular_Ring(ax, θ, map, colour)  #plotta su mappa map anello di apertura θ con asse ax
    R1 = AngleAxis(θ, 1, 1, -(ax[1]+ax[2])/ax[3]) #trovo un qualunque ortogonale ad ax
    v0 = R1*ax   #e ruoto ax di θ attorno ad esso ottenendo v0
    for ϕ in range(0, 2*π, length=100)
        R2 = AngleAxis(ϕ, ax[1], ax[2], ax[3])
        v = R2*v0   #procedo a step ruotando v0 attorno all'asse
        δ, a = vec2ang(v[1], v[2], v[3])
        map[ang2pix(map, δ, a)] = colour #man mano coloro il pixel ottenuto
    end
end

function Map_Plot_Obj(Object, start_day, sampling_time_s, lasting_time_s, map, colour)
    time0=datetime2unix(start_day)
    time_range = 0:sampling_time_s:lasting_time_s
    time_samples = collect(time_range)
    for indx in 1:length(time_samples)
        time = time0 + indx*sampling_time_s
        α, δ= Pos(Object, unix2datetime(time))
        pixel_index_obj = ang2pix(map, π/2 - δ, α) #da declinazione (latitudine) a colat.
        map[pixel_index_obj] = colour #coloro la traccia dell'oggetto
    end
end

function Map_Plot(Object, H_ID, start_day, sampling_time_s, lasting_time_s)
    map = Healpix.Map{Float64, RingOrder}(128)
    Map_Plot_Horn(H_ID, start_day, sampling_time_s, lasting_time_s, map, false)
    max_hit=maximum(map)
    Map_Plot_Obj(Object, start_day, sampling_time_s, lasting_time_s, map, max_hit+max_hit/20)
    plot(map, mollweide, title="$Object - Horn $H_ID for $lasting_time_days days")
    savefig(joinpath(output_dir,"$Object-$H_ID-$lasting_time_days.pdf"))
end

#Overloading per poter passare la mappa dell'Horn già fatta!
function Map_Plot(Object, start_day, sampling_time_s, lasting_time_s, map_name)
    map = readMapFromFITS(map_name, 1, Float64)
    max_hit=maximum(map)
    Map_Plot_Obj(Object, start_day, sampling_time_s, lasting_time_s, map, max_hit+max_hit/20)
    lasting_time_days = lasting_time_s/86400
    plot(map, mollweide, title="$Object on $map_name for $lasting_time_days days")
    savefig(joinpath(output_dir,"$Object-$map_name-$lasting_time_days.pdf"))
end

#Overloading per stampare mappe senza oggetti sopra
function Map_Plot(map_name)
    map = readMapFromFITS(map_name, 1, Float64)
    plot(map, mollweide, title="$map_name")
    savefig("$map_name.pdf")
end

function Animation(Object, H_ID::String, start_day, lasting_time_weeks)
    map = Healpix.Map{Float64, RingOrder}(128)
    Map_Plot_Horn(H_ID, start_day, 11, lasting_time_weeks*7*86400, map, false)
    max_hit=maximum(map)
    sampling_time_s = 3600
    anim = @animate for week in 0:lasting_time_weeks
        start_time_s = datetime2unix(start_day) + week*7*86400
        one_week_s = 86400*7 - sampling_time_s #tolgo l'ultimo sample perchè verrà preso dalla settimana dopo a t=0
        Map_Plot_Obj(Object, unix2datetime(start_time_s), sampling_time_s, one_week_s, map, max_hit+max_hit/20)
        plot(map, mollweide,title="$Object - Horn $H_ID week $week")
    end
    gif(anim, joinpath(output_dir,"$Object-$H_ID-animation.gif"), fps = 10) #velocità
end

#Creo un overloading di Animation, in modo da potergli anche passare la mappa (salvata in FITS)
#con la posizione dell'Horn già plottata. Mi permette di creare una volta sola la
#mappa statica dell'Horn molto più dettagliata e farci sopra animation per
#corpi diversi in modo molto più veloce. Non ha più bisogno di H_ID in imput.
function Animation(Object, start_day::DateTime, lasting_time_weeks, map_name::String)
    map = readMapFromFITS(map_name, 1, Float64)
    max_hit=maximum(map)
    sampling_time_s = 3600
    anim = @animate for week in 0:lasting_time_weeks
        start_time_s = datetime2unix(start_day) + week*7*86400
        one_week_s = 86400*7 - sampling_time_s #tolgo l'ultimo sample perchè verrà preso dalla settimana dopo a t=0
        Map_Plot_Obj(Object, unix2datetime(start_time_s), sampling_time_s, one_week_s, map, max_hit+max_hit/20)
        α, δ = Pos(Object, unix2datetime(start_time_s + one_week_s))
        ax = ang2vec(π/2 - δ, α)
        Angular_Ring(ax, deg2rad(10), map, max_hit+max_hit/20)
        plot(map, mollweide, title="$Object on $map_name week $week")
        Angular_Ring(ax, deg2rad(10), map, 0)
    end
    gif(anim, joinpath(output_dir,"$Object-$map_name-animation.gif"), fps = 10) #velocità
end

function Animation_Sun(Object, start_day::DateTime, lasting_time_weeks, map_name::String)
    map = readMapFromFITS(map_name, 1, Float64)
    max_hit=maximum(map)
    sampling_time_s = 3600
    anim = @animate for week in 0:lasting_time_weeks
        start_time_s = datetime2unix(start_day) + week*7*86400
        one_week_s = 86400*7 - sampling_time_s #tolgo l'ultimo sample perchè verrà preso dalla settimana dopo a t=0
        Map_Plot_Obj(Object, unix2datetime(start_time_s), sampling_time_s, one_week_s, map, max_hit+max_hit/20)
        Map_Plot_Obj(SUN, unix2datetime(start_time_s), sampling_time_s, one_week_s, map, max_hit*2/3)
        α, δ = Pos(SUN, unix2datetime(start_time_s + one_week_s))
        ax = ang2vec(π/2 - δ, α)
        Angular_Ring(ax, deg2rad(10), map, max_hit*2/3)
        plot(map, mollweide, title="$Object + Sun on $map_name week $week")
        map = readMapFromFITS(map_name, 1, Float64)
    end
    gif(anim, joinpath(output_dir,"$Object-Sun-$map_name-animation.gif"), fps = 10) #velocità
end

function Animation_Horns(start_day::DateTime, lasting_time_s)
    map1 = Healpix.Map{Float64, RingOrder}(128)
    start_time_s = datetime2unix(start_day)
    time_range = 0:0.5:lasting_time_s
    time_samples = collect(time_range)
    G6_dirs=Horn_dir("G6", start_day, time_samples)
    I0_dirs=Horn_dir("I0", start_day, time_samples)
    V3_dirs=Horn_dir("V3", start_day, time_samples)
    anim = @animate for indx in 1:length(time_samples)
        pixel_index = ang2pix(map1, π/2 - G6_dirs[indx, 1], G6_dirs[indx, 2])
        map1[pixel_index] += 1
        pixel_index = ang2pix(map1, π/2 - I0_dirs[indx, 1], I0_dirs[indx, 2])
        map1[pixel_index] += 1
        pixel_index = ang2pix(map1, π/2 - V3_dirs[indx, 1], V3_dirs[indx, 2])
        map1[pixel_index] += 1
        plot(map1, mollweide)
    end
    gif(anim, joinpath(output_dir,"Telescope_animation.gif"), fps=20)
end

#time0 è in formato UNIX!!!
function OneH_Simulation(Object, H_o, time0, sampling_time_s, lasting_time_s, f)
    time_range = 0:sampling_time_s:lasting_time_s
    time_samples = collect(time_range)
    H_dirs=Horn_dir(H_o, start_day, time_samples)
    for indx in 1:length(time_samples)
        time = time0 + indx*sampling_time_s
        date=unix2datetime(time)
        θ, α, δ=Dist(Object, H_dirs[indx, :], date)
        if θ < deg2rad(20)
            @printf(f, "%s,%.5e,%.5e,%.5e\n", date, θ, α, δ)
        end
    end
end

function Simulation(Object, H_ID::String, start_day)
    db = InstrumentDB()
    H_o=db.focalplane[H_ID].orientation
    sampling_time_s=0.1 #va messo 0.1
    lasting_time_s=86400
    time0=datetime2unix(start_day)
    filename=joinpath(output_dir,"$Object.$H_ID.txt")
    open(filename, "w") do f
        writedlm(f, ["date" "Ang. Dist. $Object-Horn $H_ID [rad]" "$Object's RA [rad]" "$Object's Dec [rad]"], ',')
        for day in 1:730 #2 anni = 730giorni
            OneH_Simulation(Object, H_o, time0+86400*day, sampling_time_s, lasting_time_s, f)
        end
    end
    filename
end

function Simulation(Object, H_o::Array, start_day)
    sampling_time_s=0.1 #va messo 0.1
    lasting_time_s=86400
    time0=datetime2unix(start_day)
    filename=joinpath(output_dir,"$Object.$H_ID.txt")
    open(filename, "w") do f
        writedlm(f, ["date" "Ang. Dist. $Object-Horn $H_ID [rad]" "$Object's RA [rad]" "$Object's Dec [rad]"], ',')
        for day in 1:730 #2 anni = 730giorni
            OneH_Simulation(Object, H_o, time0+86400*day, sampling_time_s, lasting_time_s, f)
        end
    end
    filename
end

function File2dates(filename, ang_lim)
    dat = readdlm(filename,  ',', skipstart=1)
    dates = dat[dat[:,2] .< deg2rad(ang_lim),1]
    return dates
end

function Dates_Union(dates_matr)
    dates = dates_matr[1]
    for indx in 2:length(dates_matr)
        dates_aux = dates_matr[indx]
        dates = union(dates, dates_aux) #faccio l'unione delle date
    end
    sort!(dates)
    return dates
end

function Dates_Intersection(dates_matr)
    dates = dates_matr[1]
    for indx in 2:length(dates_matr)
        dates_aux = dates_matr[indx]
        dates = intersect(dates, dates_aux) #faccio l'intersezione delle date
    end
    sort!(dates)
    return dates
end

function Time_Count(dates)
    rot_T_s = 60/RPM
    len = length(dates[:,1])
    time_s = 0
    i=1
    while i < len
        t0 = datetime2unix(DateTime(dates[i]))
        while datetime2unix(DateTime(dates[i]))-t0 < rot_T_s && i < len
            i += 1
        end
        time_s += rot_T_s
    end
    return time_s
end

function Planet_Result(files, ang_lims, sampling_time_s, output_name)
    filename = joinpath(output_dir,output_name*".txt")
    start_day_u=datetime2unix(start_day)
    times = []
    dates_matr = []
    for (indx, file) in enumerate(files)  #estraggo le date dai vari file
        dates = File2dates(file, ang_lims[indx])
        push!(dates_matr, dates)
    end
    dates = Dates_Union(dates_matr)    #faccio l'unione
    L_0 = 0
    for week in 1:104   #calcolo i giorni di osservazione in una settimana
        w_end_u = start_day_u + week*7*86400
        L = length(dates[dates .< string(unix2datetime(w_end_u))])
        t = L - L_0
        push!(times, t*sampling_time_s/3600)  #converto in ore
        L_0 = L
    end
    open(filename, "w") do f
        writedlm(f, times)
    end
    return times
end

function Planets_Plot(Planets)
    filename = joinpath(output_dir, "Planets.pdf")
    files = Planets.*"_results.txt"
    plot(1:104, readdlm(files[1]), label = Planets[1], title="Planets Observed ($Ang_Lim ° threshold)", xlabel="Weeks since 01/01/2022", ylabel="Observation time [h/week]")
    for indx in 2:length(files)
        times = readdlm(files[indx])
        plot!(1:104, times, label = Planets[indx])
    end
    savefig(filename)
end

function Get_Results(files, ang_lims, output_name)
    filename = joinpath(output_dir,output_name*".txt")
    times = []
    dates_matr = []
    for (indx, file) in enumerate(files) #faccio il conto dei singoli Horn
        dates = File2dates(file, ang_lims[indx])
        push!(dates_matr, dates)
        t_s = Time_Count(dates)
        push!(times, t_s/86400)
    end
    push!(files, "total")  #faccio il conto totale unendo le info
    push!(ang_lims, Ang_Lim)
    t_s = Time_Count(Dates_Union(dates_matr))
    push!(times, t_s/86400)
    percs = times./730*100
    df = DataFrame(file=files, ang_lim=ang_lims, time_days=times, percentage=percs)
    open(filename, "w") do f
    show(f, MIME("text/latex"), df, eltypes=false)
    end
end

#-------------------------   "MAIN"
#"!" * joinpath(percorso, "$H_ID.map")

files = ["SUN.G6.txt", "SUN.I0.txt", "SUN.V3.txt", "Moon.G6.txt", "Moon.I0.txt", "Moon.V3.txt", "total"]
ang_lims = [Ang_Lim, Ang_Lim, Ang_Lim, Ang_Lim, Ang_Lim, Ang_Lim, Ang_Lim]
times = [58.8361, 50.8889, 42.7972, 58.8965, 51.8431, 44.5215, 135.458]
percs = [8.05974, 6.97108, 5.86263, 8.06802, 7.10179, 6.09884, 18.5559]
df = DataFrame(file=files, ang_lim=ang_lims, time_days=times, percentage=percs)
open("total_latex.txt", "w") do f
    show(f, MIME("text/latex"), df, eltypes=false)
end

#SIMULAZIONE
Object="Neptune"
Simulation(Object, "G6", start_day)
Simulation(Object, "I0", start_day)
Simulation(Object, "V3", start_day)
Simulation(Object, "B5", start_day)

H_files=Array(["Neptune.G6.txt", "Neptune.I0.txt", "Neptune.V3.txt"])
angs=repeat(Ang_Lim:Ang_Lim, outer=length(H_files))
Get_Results(H_files, angs, "Moon_results")

H_files=Array(["SUN.G6.txt"])
angs=Array([Ang_Lim])
Get_Results(H_files, angs, "prova")

Planet_Result(H_files, angs, 0.1, Object*"_results")
Planets_Plot(Array([ "Mercury", "Venus", "Mars", "Jupiter", "Uranus", "Neptune"]))

start_date = DateTime(datas[2, 1])
plot([(DateTime(x) - start_date).value*1e-3/86400 for x in datas[2:length(datas[:,1]), 1]], ones(length(datas[:,1])-1))
#rm("$Object.$Horn.txt")

filename = joinpath(output_dir, "Sun_Moon.pdf")
plot(1:104, Planet_Result(["SUN.G6.txt", "SUN.I0.txt", "SUN.V3.txt"], [10, 10, 10], 0.1, "Sun_res"), label = "Sun", title="Sun/Moon interference with observations ($Ang_Lim ° threshold)", xlabel="Weeks since 01/01/2022", ylabel="Time discarded [h/week]")
plot!(1:104, Planet_Result(["Moon.G6.txt", "Moon.I0.txt", "Moon.V3.txt"], [10, 10, 10], 0.1, "Moon_res"), label = "Moon")
savefig(filename)

#ANIMAZIONI
#Creo la mappa per un Horn, 1 volta sola!! (su 20gg è già abbastanza definita)
lasting_time_days = 1
lasting_time_s=3600
Horns = ["I0"]
m, _ = Create_Map(Horns, lasting_time_s, "I0-1h")
plot(m, mollweide)

#Ci faccio sopra l' animazione del corpo celeste che voglio
Animation_Sun("Mercury", start_day, 104, "G6+I0+V3.map")

#ANIMAZIONE TELESCOPIO
Animation_Horns(start_day, 120)

#PLOT STATICI
Map_Plot("Saturn", start_day, 3600, 86400*365*2, "G6+I0+V3.map")
savefig(joinpath(output_dir,"$Object-$Horn-$lasting_time_days.pdf"))

Map_Plot("I0-1h.map")

map1 = Healpix.Map{Float64, RingOrder}(128)
Angular_Ring(ang2vec(π/2, π/2), deg2rad(10), map1, 4)
plot(map1, mollweide)

#Calcolo horn "fittizio"
db = InstrumentDB()
#i 3 horn che mi interessano
I0 = db.focalplane["I0"].orientation
B0 = db.focalplane["B0"].orientation
O0 = db.focalplane["O0"].orientation
#mi metto nel SR del boresight, lo zenith è 20° indietro nella direzione dell x (asse di simmetria del telesc)
zenith = [sin(deg2rad(-Ang_Zenith)), 0, cos(deg2rad(-Ang_Zenith))]
#calcolo δ: la differenza di distanza tra O0 e lo zenith (con B0 è uguale per simmetria) e I0 e lo zenith
δ = ang_dist(zenith, O0) - ang_dist(zenith, I0)
#creo delle coordinate di un fake_horn nella stessa posizione di I0 ma ruotato di δ
#per avere lo stesso angolo con lo zenith di B0 e O0.
f_h = [sin(δ), 0, cos(δ)]
#posso passare f_h al posto di H_ID grazie ad un overloading di H_dir.

#I1 I5
#O1 B5

max_stars_repo_name: LeeoBianchi/Stripointing

max_issues_repo_name: LeeoBianchi/Stripointing

max_forks_repo_name: LeeoBianchi/Stripointing


text:
#println("module JNI")
#println("import ..JavaCall: JNIEnv, JavaVM, jbyte, jchar, jshort, jint, jlong, jsize, jdouble, jfloat, jboolean")
#for t in ["jobject", "jclass", "jthrowable", "jweak", "jmethodID", "jfieldID",  "jstring", "jarray", "JNINativeMethod"]
#println("$t = Ptr{Nothing}")
#end
#for t in ["object", "boolean", "byte", "short", "int", "long", "float", "double", "char"]
#println("j$(t)Array = Ptr{Nothing}")
#end
#println("jvalue = Int64")
#println()
println("# === Below Generated by ",PROGRAM_FILE," ===")
println("")

function arg_name(m)
  m.captures[3]
end

function arg_value(m)
  #if m.captures[2] == "*" && m.captures[1] == "char" return "String($(m.captures[3]))" end
  if m.captures[3] == "env"
      "penv"
  else
     m.captures[3]
  end
end

decl_arg_type(m) = decl_arg_type(m.captures[1], m.captures[2])
ccall_arg_type(m; r=false) = ccall_arg_type(m.captures[1], m.captures[2], r=r)

function decl_arg_type(t, s)
  if s == "*"
    if t == "char"
      return "AnyString"
    elseif t == "void"
      return "Ptr{Nothing}"
    elseif t == "JNIEnv"
      return "Ptr{JNIEnv}"
    else
      return "Array{$t,1}"
#      return "Ptr{$t}"
    end
  elseif t == "jsize" #|| t == "jint" || t == "jlong" || t == "jshort" || t == "jbyte"
    return Integer
  elseif t == "jobject"
    return "jobject_arg"
  elseif t == "jobjectArray"
    return "jobjectArray_arg"
  end

  if t == "void"
    return "Nothing"
  end

  return s == "" ? t : "Array{$t,1}"
#  return s == "" ? t : "Ptr{$t}"
end

function ccall_arg_type(t, s; r=false)
  if s == "*"
    if t == "char"
      return "Cstring"
    elseif t == "void"
      return "Ptr{Nothing}"
    elseif t == "JNIEnv"
      return "Ptr{JNIEnv}"
    else
      return "Ptr{$t}"
    end
  end

  if t == "void"
    return "Nothing"
  end

  # No asterisk: type
  # If return type, Ptr
  # Else Array
  return s == "" ? t : r ? "Ptr{$t}" : "Array{$t,1}"
end

# julia_arg(m) = string(arg_name(m), "::", decl_arg_type(m))
function julia_arg(m)
    if arg_name(m) == "isCopy"
        "isCopy::PtrIsCopy"
    elseif arg_name(m) == "elems"
        string(arg_name(m), "::", "Ptr{$(m.captures[1])}")
    else
        string(arg_name(m), "::", decl_arg_type(m))
    end
end

for line in open(readlines, "jnienv.jl", "r")
  # m: match comments
  # Example:
  # # jclass ( *GetSuperclass) (JNIEnv *env, jclass sub);
  # Group 1: Return type ((?:void|char|j\w+)) "jclass"
  # Group 2: Asterisk (\**) ""
  # Group 3: Function name (\w+) "GetSuperclass"
  # Group 4: Arguments \((.*)\) "JNIEnv *env, jclass sub"
  m = match(r"\# \s* (?:const\s*)? ((?:void|char|j\w+)) \s* (\**) \s* \( \s* \* (\w+) \s* \) \s* \((.*)\) \s* ;"x, line)
  if m === nothing continue end

  # Ignore functions with variable argument syntax
  # Only process vararg functions that end in A
  if occursin("...", m.captures[4]) continue end
  if occursin("va_list", m.captures[4]) continue end

  # Get return type
  rtype = ccall_arg_type(m.captures[1], m.captures[2], r=true)

  # Function name
  fname = m.captures[3]
  # Split arguments
  args = split(m.captures[4], ",")

  # mm: Analyze arguments
  # Example: [1] "JNIEnv *env", [2] "jclass sub"
  # Group 1: Argument type ((?:void|j\w+|char|JNI\w+|JavaVM))
  # Group 2: Asterisk (\**)
  # Group 3: Argument name (\w+)
  mm = map(x->match(r"^\s* (?:const\s+)? \s* ((?:void|j\w+|char|JNI\w+|JavaVM)) \s*? (\**) \s* (\w+) \s*$"x, x), args)

  # skip the JNIEnv arg for julia since it is passed as a global to ccall
  julia_args = join(push!(map(julia_arg, mm)[2:end],"penv=ppenv[]"), ", ")
  arg_types = join(map(ccall_arg_type, mm), ", ")
  arg_names = join(map(arg_value, mm), ", ")

  # Commented out export command
  # print("#export $fname\n")
  print("$fname($julia_args) =\n  ccall(jniref[].$(fname), $rtype, ($arg_types,), $arg_names)\n\n")
end

# println("end")
println("")
println("# === Above Generated by ",PROGRAM_FILE," ===")


text:
function get_playback_policy(d, s, mdp, backup_policy)
    logpdfs = [logpdf(d, s, i) for i=1:length(first(s)[2])]
    PlaybackPolicy(actions(mdp)[s[:a]], logpdfs, backup_policy)
end

function sample_playback_policy_fn(d, mdp, backup_policy, rng = Random.GLOBAL_RNG)
    () -> get_playback_policy(d, rand(rng, d), mdp, backup_policy)
end

function mdp_loss(mdp, backup_policy)
    function loss(d, s)
        p = get_playback_policy(d, s, mdp, backup_policy)
        -simulate(RolloutSimulator(), mdp, p)
    end
end

function mdp_weight(mdp)
    function weight(d, s)
        exp(logpdf(mdp, actions(mdp)[s[:a]]) - logpdf(d, s))
    end
end

function continous_actions_weight(d_true)
    function weight(d, s)
        exp(logpdf(d_true, s) - logpdf(d, s))
    end
end


text:
module NoisyMNIST

    using Pkg
    using MLDatasets
    using JLD2
    using ..Misc
    using ..MNIST

    const FILEPATH = joinpath(Misc.DATASETS_FOLDER, "noisybitsmnist.jld2")
    const NOISE_RATE = 0.3
    const NOISE_LABEL = -1

    function addnoise(set)
        for line in set
            for i in eachindex(line)
                if rand() < NOISE_RATE
                    line[i] != rand(Bool)
                end
            end
        end
    end

    """
        download()

    Forces downloading the neisy bits mnist dataset
    """
    function download()
        set = mnist()
        addnoise(set["train_x"])
        addnoise(set["test_x"])
        fill!(set["train_y"], NOISE_LABEL)
        fill!(set["test_y"], NOISE_LABEL)
        save(FILEPATH, set)
        set
    end

    """
        noisymnist()

    Downloads and return the noisy bits mnist dataset
    """
    function noisymnist()
        if isfile(FILEPATH)
            return load(FILEPATH)
        else
            return download()
        end
    end

    export noisymnist, NOISE_LABEL

end


text:
"""
    WeightedVariance(T = Float64)

Simple weighted variance, tracked as type `T`.

# Example:
    o = fit!(WeightedVariance(), rand(100), rand(100))
    sum(o)
    mean(o)
    var(o)
    std(o)
"""
mutable struct WeightedVariance{T} <: WeightedOnlineStat{T}
    μ::T
    σ2::T
    W::T
    W2::T
    n::Int
    function WeightedVariance{T}(
            μ = T(0), σ2 = T(0), W = T(0), W2 = T(0), n = 0
        ) where T
        new{T}(T(μ), T(σ2), T(W), T(W2), Int(n))
    end
end

WeightedVariance(μ::T, σ2::T, W::T, W2::T, n::Int) where T =
    WeightedVariance{T}(μ, σ2, W, W2, n)
WeightedVariance(::Type{T}) where T = WeightedVariance(T(0), T(0), T(0), T(0), 0)
WeightedVariance() = WeightedVariance(Float64)

function OnlineStatsBase._fit!(o::WeightedVariance{T}, x, w) where T
    xx = convert(T, x)
    ww = convert(T, w)

    o.n += 1
    γ1 = T(1) / o.n
    o.W = smooth(o.W, ww, γ1)
    o.W2 = smooth(o.W2, ww * ww, γ1)
    γ = ww / (o.W * o.n)
    μ = o.μ

    o.μ = smooth(o.μ, xx, γ)
    o.σ2 = smooth(o.σ2, (xx - o.μ) * (xx - μ), γ)

    return o
end

function OnlineStatsBase._merge!(
        o::WeightedVariance{T},
        o2::WeightedVariance
    ) where T

    o2_μ = convert(T, o2.μ)
    o2_σ2 = convert(T, o2.σ2)
    o2_W = convert(T, o2.W)
    o2_W2 = convert(T, o2.W2)

    n = o.n + o2.n
    W = smooth(o.W, o2_W, o2.n / n)
    γ1 = (o.W * o.n) / (W * n)
    γ2 = (o2_W * o2.n) / (W * n)

    μ = smooth(o.μ, o2_μ, γ2)

    # o.σ2 =
    #     γ1 * ( o.σ2 +  o.μ ^ 2) +
    #     γ2 * (o2_σ2 + o2_μ ^ 2) -
    #     μ ^ 2
    o.σ2 =
        γ1 * ( o.σ2  + (o.μ - μ) ^ 2) +
        γ2 * (o2_σ2 + (o2_μ - μ) ^ 2)

    o.n = n
    o.μ = μ
    o.W = W
    o.W2 = OnlineStats.smooth(o.W2, o2_W2, o2.n / o.n)

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

    # μ = o.μ
    #
    # γ = o2_W / (o2_W + o.W)
    # δ = o2_μ - o.μ
    #
    # o.σ2 = OnlineStats.smooth(o.σ2, o2_σ2, γ) + (δ ^ 2) * γ * (1.0 - γ)
    # o.μ = OnlineStats.smooth(o.μ, o2_μ, γ)
    # # o.σ2 = o.σ2 + (o.W + o2_W) * (μ)
    #
    # o.W += o2_W
    # o.W2 += o2_W2

    return o
end

OnlineStatsBase.value(o::WeightedVariance) = o.σ2
Base.sum(o::WeightedVariance) = mean(o) * meanweight(o) * nobs(o)
Statistics.mean(o::WeightedVariance) = o.μ
function Statistics.var(
        o::WeightedVariance;
        corrected = true,
        weight_type = :analytic
    )
    if corrected
        if weight_type == :analytic
            value(o) / (1 - (o.W2 * nobs(o)) / (weightsum(o) ^ 2))
        elseif weight_type == :frequency
            value(o) / (weightsum(o) - 1) * weightsum(o)
        elseif weight_type == :probability
            error("If you need this, please make a PR or open an issue")
        else
            throw(ArgumentError("weight type $weight_type not implemented"))
        end
    else
        value(o)
    end
end
Statistics.std(o::WeightedVariance; kw...) = sqrt.(var(o; kw...))
Base.copy(o::WeightedVariance) = WeightedVariance(o.μ, o.σ2, o.W, o.W2, o.n)


text:
using Pkg
Pkg.activate(@__DIR__)
CI = get(ENV, "CI", nothing) == "true"
ENV["PYTHON"] = ""
Pkg.build("PyCall")
using Documenter, ClimatePlots

makedocs(sitename = "ClimatePlots.jl",
    doctest = false,
    format = Documenter.HTML(
    prettyurls = CI,
    ),
    pages = [
        "index.md",
        ]

)

if CI
    deploydocs(
    repo   = "github.com/JuliaClimate/ClimatePlots.jl.git",
    target = "build"
    )
end


text:
@testset "check_interval" begin
    @testset "function" begin
        # Unique zero in the interval
        @test ArbExtras.check_interval(sin, Arf(-1), Arf(1)) == (true, true)
        @test ArbExtras.check_interval(sin, Arf(-1e-20), Arf(1e-20)) == (true, true)
        @test ArbExtras.check_interval(cos, Arf(1), Arf(2)) == (true, true)
        @test ArbExtras.check_interval(x -> sin(π * x / 2), Arf(1.5), Arf(2.5)) ==
              (true, true)
        @test ArbExtras.check_interval(
            x -> sin(π * x / 2),
            Arf(2 - 1e-10),
            Arf(2 + 1e-10),
        ) == (true, true)
        @test ArbExtras.check_interval(
            x -> (x - 1) * (x - 2) * (x - 3),
            Arf(1.9),
            Arf(2.1),
        ) == (true, true)

        # No zero in the interval
        @test ArbExtras.check_interval(sin, Arf(1e-10), Arf(2e-10)) == (false, false)
        @test ArbExtras.check_interval(sin, Arf(1), Arf(3)) == (false, false)
        @test ArbExtras.check_interval(cos, Arf(-1), Arf(1)) == (false, false)
        @test ArbExtras.check_interval(x -> sin(π * x / 2), Arf(0.5), Arf(1.5)) ==
              (false, false)
        @test ArbExtras.check_interval(
            x -> (x - 1) * (x - 2) * (x - 3),
            Arf(1.4),
            Arf(1.6),
        ) == (false, false)

        # Multiple simple zeros in the interval
        @test ArbExtras.check_interval(sin, Arf(-1), Arf(4)) == (true, false)
        @test ArbExtras.check_interval(sin, Arf(-4), Arf(1)) == (true, false)
        @test ArbExtras.check_interval(cos, Arf(-2), Arf(2)) == (true, false)
        @test ArbExtras.check_interval(x -> sin(π * x / 2), Arf(1), Arf(2)) == (true, false)
        @test ArbExtras.check_interval(x -> sin(π * x / 2), Arf(2), Arf(3)) == (true, false)
        @test ArbExtras.check_interval(
            x -> (x - 1) * (x - 1 + eps()),
            Arf(1 - eps()),
            Arf(1 + eps()),
        ) == (true, false)

        # Zero of higher order
        @test ArbExtras.check_interval(x -> x^2, Arf(-eps()), Arf(eps())) == (true, false)

        # check_unique = false
        @test ArbExtras.check_interval(sin, Arf(-1), Arf(1), check_unique = false) ==
              (true, false)
        @test ArbExtras.check_interval(cos, Arf(-1), Arf(1), check_unique = false) ==
              (false, false)
        @test ArbExtras.check_interval(sin, Arf(-1), Arf(4), check_unique = false) ==
              (true, false)
        @test ArbExtras.check_interval(
            x -> x^2,
            Arf(-eps()),
            Arf(eps()),
            check_unique = false,
        ) == (true, false)
    end

    @testset "ArbPoly" begin
        x = ArbPoly([1, 0])
        # Unique zero in the interval
        @test ArbExtras.check_interval(Arblib.fromroots(ArbPoly, [0]), Arf(-1), Arf(1)) ==
              (true, true)
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [-2, 0, 2]),
            Arf(-1),
            Arf(1),
        ) == (true, true)
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [0], [2im]),
            Arf(-1),
            Arf(1),
        ) == (true, true)

        # No zero in the interval
        @test ArbExtras.check_interval(ArbPoly(1), Arf(-1), Arf(1)) == (false, false)
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [-2, 2]),
            Arf(-1),
            Arf(1),
        ) == (false, false)

        # Multiple simple zeros in the interval
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [-0.5, 0.5]),
            Arf(-1),
            Arf(1),
        ) == (true, false)
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [-1, 1]),
            Arf(-1),
            Arf(1),
        ) == (true, false)

        # Zero of higher order
        @test ArbExtras.check_interval(ArbPoly(0), Arf(-1), Arf(1)) == (true, false)
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [0, 0]),
            Arf(-1),
            Arf(1),
        ) == (true, false)

        # check_unique = false
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [0]),
            Arf(-1),
            Arf(1),
            check_unique = false,
        ) == (true, false)
        @test ArbExtras.check_interval(ArbPoly(1), Arf(-1), Arf(1), check_unique = false) ==
              (false, false)
        @test ArbExtras.check_interval(
            Arblib.fromroots(ArbPoly, [-0.5, 0.5]),
            Arf(-1),
            Arf(1),
            check_unique = false,
        ) == (true, false)
        @test ArbExtras.check_interval(ArbPoly(0), Arf(-1), Arf(1), check_unique = false) ==
              (true, false)
    end
end

@testset "isolate_roots" begin
    @testset "sin(πx/2)" begin
        f = x -> sin(π * x / 2)
        a = -9.5
        b = 19.5
        roots = -8:2:18
        found, flags = ArbExtras.isolate_roots(f, Arf(a), Arf(b))

        # Check that all roots are accounted for
        found_balls = Arb.(found)
        @test all(any(contains(root), found_balls) for root in roots)

        # Check that all roots are reported as unique and that they
        # indeed are
        @test all(flags)
        @test all(isone(count(contains(root), found_balls)) for root in roots)
    end

    @testset "sin(1/x)" begin
        f = x -> sin(1 / x)
        a = 0
        b = 1
        found1, flags1 = ArbExtras.isolate_roots(f, Arf(a), Arf(b))
        found2, flags2 = ArbExtras.isolate_roots(f, Arf(a), Arf(b), depth = 15)

        # Check that the interval containing zero is included and that
        # it's smaller if we bisect more
        @test !flags1[1] && iszero(found1[1][1])
        @test !flags2[1] && iszero(found2[1][1])
        @test found1[1][2] > found2[1][2]

        # Check that we find isolate more zeros with higher depth
        @test sum(flags1) < sum(flags2)
    end

    @testset "thin interval" begin
        @test ArbExtras.isolate_roots(sin, Arf(0), Arf(0)) == ([(Arf(0), Arf(0))], [true])
        @test ArbExtras.isolate_roots(x -> sin(x) - sin(x), Arf(1), Arf(1)) ==
              ([(Arf(1), Arf(1))], [false])
        @test all(isempty.(ArbExtras.isolate_roots(cos, Arf(0), Arf(0))))
    end

    # Check that the derivative is not computed when check_unique is
    # false and that it splits the expected number of times
    @testset "check_unique = false" begin
        depth = 10
        found, flags = ArbExtras.isolate_roots(
            x -> Arb(0),
            Arf(0),
            Arf(1),
            check_unique = false,
            depth = depth,
        )
        @test length(found) == 2^(depth = 10 - 1)
        @test !any(flags)
        @test all(found[i][2] == found[i+1][1] for i = 1:length(found)-1)
    end

    @testset "ArbPoly" begin
        roots = Arb[-5, 0, ℯ, π, 10]
        p = Arblib.fromroots(ArbPoly, roots)
        a = -10
        b = 20
        found, flags = ArbExtras.isolate_roots(p, Arf(a), Arf(b), depth = 15)

        # Check that all roots are accounted for
        found_balls = Arb.(found)
        @test all(any(contains(root), found_balls) for root in roots)

        # Check that all roots are reported as unique and that they
        # indeed are
        @test all(flags)
        @test all(isone(count(contains(root), found_balls)) for root in roots)
    end

    @testset "ArbPoly - check_unique = false" begin
        depth = 10
        found, flags = ArbExtras.isolate_roots(
            ArbPoly(0),
            Arf(0),
            Arf(1),
            check_unique = false,
            depth = depth,
        )
        @test length(found) == 2^(depth = 10 - 1)
        @test !any(flags)
        @test all(found[i][2] == found[i+1][1] for i = 1:length(found)-1)
    end
end


text:
using YaoBase, YaoArrayRegister

export ReflectGate, reflect

"""
    ReflectGate{N, T, Tr} <: PrimitiveBlock{N}

Reflection operator to target state `psi`.

# Definition

```math
|ψ⟩ → 2 |s⟩⟨s| - 1
```
"""
struct ReflectGate{N,T,Tr<:ArrayReg{1,T}} <: PrimitiveBlock{N}
    psi::Tr
    function ReflectGate{N,T,Tr}(psi) where {N,T,Tr}
        @warn "`ReflectGate` will be moved to `YaoExtensions.jl` in the next release."
        new{N,T,Tr}(psi)
    end
end

"""
    ReflectGate(r::ArrayReg{1})

Create a [`ReflectGate`](@ref) with a quantum register `r`.
"""
ReflectGate(r::ArrayReg{1,T}) where {T} = ReflectGate{nqubits(r),T,typeof(r)}(r)

"""
    ReflectGate(r::AbstractVector)

Create a [`ReflectGate`](@ref) with a quantum state vector `v`.
"""
ReflectGate(v::AbstractVector{<:Complex}) = ReflectGate(ArrayReg(v))

"""
    reflect(r::ArrayReg)

Create a [`ReflectGate`](@ref) with an [`ArrayReg`](@ref).
"""
reflect(r::ArrayReg) = reflect(statevec(r))

"""
    reflect(v::AbstractVector{<:Complex})

Create a [`ReflectGate`](@ref) with an quantum state vector `v`.

# Example

```jldoctest; setup=:(using YaoBlocks; using YaoArrayRegister)
julia> reflect(rand_state(3))
┌ Warning: `ReflectGate` will be moved to `YaoExtensions.jl` in the next release.
└ @ YaoBlocks ~/.julia/dev/YaoBlocks/src/primitive/reflect_gate.jl:19
reflect(ArrayReg{1, Complex{Float64}, Array...})
```
"""
reflect(v::AbstractVector{<:Complex}) = ReflectGate(v)

function apply!(r::ArrayReg, g::ReflectGate{N,T,<:ArrayReg}) where {N,T}
    v = state(g.psi)
    r.state .= 2 .* (v' * r.state) .* v - r.state
    return r
end

# target type is the same with block's
function mat(::Type{T}, r::ReflectGate{N,T}) where {N,T}
    v = statevec(r.psi)
    return 2 * v * v' - IMatrix(length(v))
end

# different
function mat(::Type{T1}, r::ReflectGate{N,T2}) where {N,T1,T2}
    M = mat(T2, r)
    return copyto!(similar(M, T1), M)
end

Base.:(==)(A::ReflectGate, B::ReflectGate) = A.psi == B.psi
Base.copy(r::ReflectGate) = ReflectGate(r.psi)

YaoBase.isreflexive(::ReflectGate) = true
YaoBase.ishermitian(::ReflectGate) = true
YaoBase.isunitary(::ReflectGate) = true