text:
"""
Run Gibbs sampler on masked spike train. Alternate
between imputing data in masked regions and updating
the model through classic gibbs_sample!(...) function.
"""
function masked_gibbs!(
        model::Union{SeqModel,DistributedSeqModel},
        masked_spikes::Vector{Spike},
        unmasked_spikes::Vector{Spike},
        masks::Vector{Mask},
        initial_assignments::Vector{Int64},
        num_spike_resamples::Int64,
        samples_per_resample::Int64,
        extra_split_merge_moves::Int64,
        split_merge_window::Float64,
        save_every::Int64;
        verbose::Bool=true
    )

    sampled_spikes = Spike[]
    sampled_assignments = Int64[]

    # Compute proportion of the data that is masked.
    masked_proportion = 0.0
    for (_, (t0, t1)) in masks
        masked_proportion += (t1 - t0)
    end
    masked_proportion /= (model.max_time * num_neurons(model))
    @show masked_proportion

    # Create inverted masks to compute train log likelihood.
    inv_masks = compute_complementary_masks(
        masks, num_neurons(model), model.max_time)

    # Sanity check.
    assert_spikes_in_mask(masked_spikes, masks)
    assert_spikes_in_mask(unmasked_spikes, inv_masks)
    assert_spikes_not_in_mask(masked_spikes, inv_masks)
    assert_spikes_not_in_mask(unmasked_spikes, masks)

    # Compute log-likelihood of a homogeneous Poisson process in
    # the train and test sets.
    train_baseline = homogeneous_baseline_log_like(unmasked_spikes, inv_masks)
    test_baseline = homogeneous_baseline_log_like(masked_spikes, masks)

    @show train_baseline test_baseline

    n_unmasked = length(unmasked_spikes)
    assignment_hist = zeros(Int64, n_unmasked, 0)
    train_log_p_hist = Float64[]
    test_log_p_hist = Float64[]
    latent_event_hist = Vector{EventSummaryInfo}[]
    globals_hist = SeqGlobals[]

    unmasked_assignments = initial_assignments

    for i = 1:num_spike_resamples

        # Sample new spikes in each masked region.
        sample_masked_spikes!(
            sampled_spikes,
            sampled_assignments,
            model,
            masks
        )

        # assert_spikes_in_mask(sampled_spikes, masks)
        # assert_spikes_not_in_mask(sampled_spikes, inv_masks)

        # Run gibbs sampler. Note that sufficient statistics are
        # recomputed at the beginning of gibb_sample!(...) so all
        # events will have the appropriate spike assignments / initialization
        (
            assignments,
            _assgn_hist,
            _lp_hist,
            _latents,
            _globals
        ) =
        gibbs_sample!(
            model,
            vcat(unmasked_spikes, sampled_spikes),
            vcat(unmasked_assignments, sampled_assignments),
            samples_per_resample,
            extra_split_merge_moves,
            split_merge_window,
            save_every;
            verbose=false
        )

        # Update initial assignments for next Gibbs sampler run.
        unmasked_assignments .= view(assignments, 1:n_unmasked)

        # Save history
        assignment_hist = cat(
            assignment_hist,
            view(_assgn_hist, 1:n_unmasked, :),
            dims=2
        )

        # Evaluate model likelihood on observed spikes.
        push!(
            train_log_p_hist,
            log_like(model, unmasked_spikes, inv_masks) - train_baseline
        )

        # Evaluate model likelihood on heldout spikes.
        push!(
            test_log_p_hist,
            log_like(model, masked_spikes, masks) - test_baseline
        )

        append!(latent_event_hist, _latents)
        append!(globals_hist, _globals)

        verbose && print(i * samples_per_resample, "-")
        flush(stdout)

    end

    verbose && println("Done")
    flush(stdout)

    # Before returning, remove assignments assigned to imputed spikes.
    recompute!(model, unmasked_spikes, unmasked_assignments)

    # Rescale train and test log likelihoods.
    train_log_p_hist ./= ((1 - masked_proportion) * model.max_time * num_neurons(model))
    test_log_p_hist ./= (masked_proportion * model.max_time * num_neurons(model))

    return (
        unmasked_assignments,
        assignment_hist,
        train_log_p_hist,
        test_log_p_hist,
        latent_event_hist,
        globals_hist
    )

end

"""
Annealed Gibbs sampling with masked data.
"""
function annealed_masked_gibbs!(
        model::Union{SeqModel,DistributedSeqModel},
        spikes::Vector{Spike},
        masks::Vector{Mask},
        initial_assignments::Vector{Int64},
        num_anneals::Int64,
        max_temperature::Float64,
        num_spike_resamples_per_anneal::Int64,
        samples_per_resample::Int64,
        extra_split_merge_moves::Int64,
        split_merge_window::Float64,
        save_every::Int64;
        verbose::Bool=true
    )

    masked_spikes, unmasked_spikes = split_spikes_by_mask(spikes, masks)

    target_mean = mean(priors(model).seq_event_amplitude)
    target_var = var(priors(model).seq_event_amplitude)

    temperatures = exp10.(range(log10(max_temperature), 0, length=num_anneals))

    unmasked_assignments = fill(-1, length(unmasked_spikes))
    assignment_hist = zeros(Int64, length(unmasked_spikes), 0)
    train_log_p_hist = Float64[]
    test_log_p_hist = Float64[]
    latent_event_hist = Vector{EventSummaryInfo}[]
    globals_hist = SeqGlobals[]

    for temp in temperatures

        # Print progress.
        verbose && println("TEMP:  ", temp)
        flush(stdout)

        # Anneal prior on sequence amplitude.
        prior = priors(model)
        set_priors!(
            model,
            SeqPriors(
                prior.seq_event_rate,
                prior.seq_type_proportions,
                specify_gamma(target_mean, target_var * temp), # Anneal prior.
                prior.neuron_response_proportions,
                prior.neuron_response_profile,
                prior.bkgd_amplitude,
                prior.bkgd_proportions,
                prior.warp_values,
                prior.warp_log_proportions
            )
        )

        # Recompute probability of introducing a new cluster.
        #  ==> TODO: maybe set_priors! should do this automatically?
        prior = priors(model)
        α = prior.seq_event_amplitude.α
        β = prior.seq_event_amplitude.β
        λ = prior.seq_event_rate
        set_new_cluster_log_prob!(
            model,
            log(α) + log(λ) + log(model.max_time) + α * (log(β) - log(1 + β))
        )
        model.bkgd_log_prob = (
            log(model.globals.bkgd_amplitude)
            + log(model.max_time)
            + log(1 + β)
        )

        # Draw gibbs samples.
        (
            unmasked_assignments,
            _assgn,
            _train_hist,
            _test_hist,
            _latents,
            _globals
        ) =
        masked_gibbs!(
            model,
            masked_spikes,
            unmasked_spikes,
            masks,
            unmasked_assignments,
            num_spike_resamples_per_anneal,
            samples_per_resample,
            extra_split_merge_moves,
            split_merge_window,
            save_every;
            verbose=verbose
        )

        # Save samples.
        assignment_hist = cat(assignment_hist, _assgn, dims=2)
        append!(train_log_p_hist, _train_hist)
        append!(test_log_p_hist, _test_hist)
        append!(latent_event_hist, _latents)
        append!(globals_hist, _globals)

    end

    return (
        unmasked_assignments,
        assignment_hist,
        train_log_p_hist,
        test_log_p_hist,
        latent_event_hist,
        globals_hist
    )

end

function masked_gibbs!(
        model::Union{SeqModel,DistributedSeqModel},
        spikes::Vector{Spike},
        masks::Vector{Mask},
        initial_assignments::Vector{Int64},
        num_spike_resamples::Int64,
        samples_per_resample::Int64,
        extra_split_merge_moves::Int64,
        split_merge_window::Float64,
        save_every::Int64;
        verbose::Bool=true
    )

    masked_spikes, unmasked_spikes = split_spikes_by_mask(spikes, masks)

    return masked_gibbs!(
        model,
        masked_spikes,
        unmasked_spikes,
        masks,
        initial_assignments,
        num_spike_resamples,
        samples_per_resample,
        extra_split_merge_moves,
        split_merge_window,
        save_every;
        verbose=verbose
    )
end

"""
Impute missing data
"""
function sample_masked_spikes!(
        spikes::Vector{Spike},
        assignments::Vector{Int64},
        model::SeqModel,
        masks::Vector{Mask}
    )

    empty!(spikes)
    empty!(assignments)
    globals = model.globals

    # Sample background spikes.
    S_bkgd = rand(Poisson(
        globals.bkgd_amplitude * model.max_time))
    bkgd_dist = Categorical(exp.(globals.bkgd_log_proportions))

    n_bkgd = rand(bkgd_dist, S_bkgd)
    t_bkgd = rand(S_bkgd) * model.max_time

    for (sampled_n, sampled_t) in zip(n_bkgd, t_bkgd)
        for (n, (start, stop)) in masks
            if (start < sampled_t < stop) && (sampled_n == n)
                push!(spikes, Spike(n, sampled_t))
                push!(assignments, -1)
                break
            end
        end
    end

    # Compute neuron probabilities.
    neuron_rel_amps =
        exp.(globals.neuron_response_log_proportions)
    neuron_dists =
        [Categorical(neuron_rel_amps[:, i]) for i = 1:size(neuron_rel_amps, 2)]

    # Sample sequence-evoked spikes.
    for (k, event) in enumerate(model.sequence_events)

        # Assignment id for latent event.
        z = model.sequence_events.indices[k]

        # Num spikes evoked by latent event.
        S = rand(Poisson(event.sampled_amplitude))

        # Sample neuron, then spike time, for each spike.
        for n in rand(neuron_dists[event.sampled_type], S)

            # Sample spike time.
            μ = globals.neuron_response_offsets[n, event.sampled_type]
            σ = sqrt(globals.neuron_response_widths[n, event.sampled_type])
            w = model.priors.warp_values[event.sampled_warp]
            t = w * (σ * randn() + μ) + event.sampled_timestamp

            # Check if spike is inside a masked region.
            for (n_mask, (mask_start, mask_stop)) in masks
                if (mask_start < t < mask_stop) && (n == n_mask)
                    push!(spikes, Spike(n, t))
                    push!(assignments, z)
                    break
                end
            end

        end
    end

    return spikes, assignments

end

sample_masked_spikes!(
    spikes::Vector{Spike},
    assignments::Vector{Int64},
    model::DistributedSeqModel,
    masks::Vector{Mask}
) = sample_masked_spikes!(spikes, assignments, model.primary_model, masks)

# ===========
#
# Helper functions to create masks, split spikes.
#
# ===========

function compute_complementary_masks(
        masks::Vector{Mask},
        num_neurons::Int64,
        max_time::Float64
    )

    inverted_masks = [[(n, (0.0, max_time))] for n in 1:num_neurons]

    for (n, (t0, t1)) in masks
        for i = 1:length(inverted_masks[n])
            n_, (v0, v1) = inverted_masks[n][i]
            @assert n_ == n
            if (t0 >= v0) && (t1 <= v1)
                deleteat!(inverted_masks[n], i)
                push!(
                    inverted_masks[n],
                    (n, (v0, t0))
                )
                push!(
                    inverted_masks[n],
                    (n, (t1, v1))
                )
                break
            end
            @assert (i + 1) != length(inverted_masks)
        end
    end

    return vcat(inverted_masks...)

end

function create_random_mask(
        num_neurons::Integer,
        max_time::Real,
        mask_lengths::Real,
        percent_masked::Real
    )

    @assert num_neurons > 0
    @assert max_time > mask_lengths
    @assert mask_lengths > 0
    @assert 0 <= percent_masked < 100

    T_masked = percent_masked * max_time * num_neurons / 100.0
    n_masks = Int(round(T_masked / mask_lengths))
    intervals = Tuple{Float64, Float64}[]
    for start in range(0, max_time - mask_lengths, step=mask_lengths)
        push!(intervals, (start, start + mask_lengths))
    end

    sample(
        collect(Iterators.product(1:num_neurons, intervals)),
        n_masks, replace=false
    )

end

function create_blocked_mask(
        num_neurons::Integer,
        max_time::Real
    )

    masked_neurons = sample(1:num_neurons, num_neurons ÷ 2, replace=false)

    masked_intervals = vcat(
        repeat([(0.0, max_time / 2)], div(num_neurons ÷ 2, 2)),
        repeat([(max_time / 2, max_time)], cld(num_neurons ÷ 2, 2))
    )

    return collect(zip(masked_neurons, masked_intervals))

end

function split_spikes_by_mask(
        spikes::Vector{Spike},
        masks::Vector{Mask}
    )
    # Save list of spikes that are masked out.
    masked_spikes = Spike[]
    unmasked_spikes = Spike[]

    for x in spikes

        # See if x falls within any masked region.
        is_masked = false
        for (n, (start, stop)) in masks
            if (start < x.timestamp < stop) && (x.neuron == n)
                is_masked = true
                push!(masked_spikes, x)
                break
            end
        end

        # Mark x as unmasked if no match was found.
        !(is_masked) && push!(unmasked_spikes, x)
    end

    return masked_spikes, unmasked_spikes
end

function assert_spikes_in_mask(
        spikes::Vector{Spike},
        masks::Vector{Mask}
    )
    # Check that all spikes are in masked region.
    for x in spikes
        is_in_mask = false
        for (n, (t0, t1)) in masks
            if (t0 < x.timestamp < t1) && (x.neuron == n)
                is_in_mask = true
            end
        end
        if !(is_in_mask)
            @show x
            @assert false "spike is falsely excluded from mask..."
        end
    end
end

function assert_spikes_not_in_mask(
        spikes::Vector{Spike},
        masks::Vector{Mask}
    )
    # Check that all spikes are in masked region.
    for x in spikes
        for (n, (t0, t1)) in masks
            if (t0 < x.timestamp < t1) && (x.neuron == n)
                @show x
                @show (n, (t0, t1))
                @assert false "spike is falsely included in mask..."
            end
        end
    end
end

function clean_masks(
    masks::Vector{Mask},
    num_neurons::Integer
)
    new_masks = Mask[]
    for n = 1:num_neurons
        # Find the masks corresponding to thie neuron
        this_neur_masks = Mask[]
        for mask = 1:length(masks)
            if masks[mask][1] == n
                push!(this_neur_masks, masks[mask])
            end
        end

        # Find the start of mask times for this neuron
        start_times = []
        for (n,(t0,t1)) in this_neur_masks
            push!(start_times, t0)
        end

        # Sort them
        start_times_indices = sortperm(start_times)
        sorted_mask = this_neur_masks[start_times_indices]

        # Find the offending masks
        indices_to_delete = []
        for index = 1:(length(sorted_mask)-1)
            if (sorted_mask[index+1][2][1] - sorted_mask[index][2][2] < 0.0001)
                push!(indices_to_delete, index)
            end
        end

        # Go through and remove them
        counter = 0
        for index = 1:length(sorted_mask)
            if index in indices_to_delete
                counter = counter + 1
            else
                push!(new_masks, (sorted_mask[index][1],(sorted_mask[index-counter][2][1],sorted_mask[index][2][2])))
                counter = 0
            end
        end
    end

    return new_masks
end


text:
# timer.jl

using LazyGrids: btime, @timeo
using BenchmarkTools: @benchmark
using Test: @test, @testset, @test_throws, @inferred

@testset "timer" begin
	t = @benchmark sum(1:7)
	@inferred btime(t)
	@test btime(t) isa String

	@test (@timeo sum(1:7)) isa String
end

ext: jl


text:
#closure to keep data local
function func_gen(f::Function, xs::Vector, ys::Vector)
    return function(x::Vector, g::Vector)
        if length(g) > 0
            g = ForwardDiff.Gradient(f,x)
        end
        for i = 1:length(x)
            print("x[$i]: ",x[i], ", ")
            if i == length(x)
                print("\n")
            end
        end
        return sum(abs2.(f(x).(xs)-ys))
    end
end

##### Positional arguments
# xs::Vector - x values to be passed to f in objective
# ys::Vector - y values to be used in objective
# f::function - function whose parameters should be determined
# alg::Symbol - algorithm
# n::Int - dimension of optimization space
# lbs::Vector - lower bounds
# x0::Vector - initial guess
##### Keyword arguments
# ecs::Vector{Function} - equality constraints
function optim_data(f, xs, ys, alg, n, lbs, ubs, x0; ecs=Vector{}(), ics=Vector{}())
    opt = NLopt.Opt(alg, n)
    NLopt.lower_bounds!(opt, lbs)
    NLopt.upper_bounds!(opt, ubs)

    NLopt.min_objective!(opt, func_gen(f,xs,ys))
    for i=1:length(ecs)
        NLopt.equality_constraint!(opt, ecs[i])
    end
    for i=1:length(ics)
        NLopt.equality_constraint!(opt, ics[i])
    end

    return NLopt.optimize(opt, x0)
end


text:
PlotlyBase.trace_map(p::SyncPlot, axis) = trace_map(p.plot, axis)
JSON.lower(sp::SyncPlot) = sp.plot

PlotlyBase._is3d(p::SyncPlot) = _is3d(p.plot)

# subplot methods on syncplot
Base.hcat(sps::SyncPlot...) = SyncPlot(hcat([sp.plot for sp in sps]...))
Base.vcat(sps::SyncPlot...) = SyncPlot(vcat([sp.plot for sp in sps]...))
Base.vect(sps::SyncPlot...) = vcat(sps...)
Base.hvcat(rows::Tuple{Vararg{Int}}, sps::SyncPlot...) =
    SyncPlot(hvcat(rows, [sp.plot for sp in sps]...))

function PlotlyBase.add_recession_bands!(p::SyncPlot; kwargs...)
    new_shapes = add_recession_bands!(p.plot; kwargs...)
    relayout!(p, shapes=new_shapes)
    new_shapes
end


text:
function updateF!(F::AbstractMatrix, mechanism::Mechanism)
    N = length(mechanism.bodies)*6
    for eqc in mechanism.eqconstraints
        L = length(eqc)
        ids = unique(eqc.childids)
        if eqc.parentid !== nothing
            push!(ids, eqc.parentid)
        end
        for id in ids
            Gᵥ = ConstrainedDynamics.∂g∂ʳvel(mechanism, eqc, id)
            offset = (id-1)*6
            F[N+1:N+L, 1+offset:6+offset] = Gᵥ
            F[1+offset:6+offset, N+1:N+L] = Gᵥ'
        end
        N += L  # Advance along first dimension
    end
    return nothing
end

function updateS!(s::AbstractArray, v::AbstractArray, ω::AbstractArray)
    N = length(v)
    for (iv, is) in zip(1:N, 1:6:N*6)
        s[is:is+2] = v[iv]
        s[is+3:is+5] = ω[iv]
    end
end

function updateMechanism!(mechanism::Mechanism, s::AbstractArray)
    for (id, body) in enumerate(mechanism.bodies)
        state = body.state
        offset = (id-1)*6
        state.vsol[2] = s[1+offset:3+offset]
        state.ωsol[2] = s[4+offset:6+offset]
    end
    Nbodies = length(mechanism.bodies)
    offset = Nbodies
    for eqc in mechanism.eqconstraints
        eqcdim = length(eqc.λsol[1])
        eqc.λsol[2] = s[1+offset:eqcdim+offset]
        offset += eqcdim
    end
end

function getvel(s::AbstractArray, Nbodies::Integer)
    v, ω = Vector{SVector{3, Float64}}(), Vector{SVector{3, Float64}}()
    for sid in 1:6:Nbodies*6
        push!(v, SVector{3, Float64}(s[sid:sid+2]))
        push!(ω, SVector{3, Float64}(s[sid+3:sid+5]))
    end
    return v, ω
end

function d(s::AbstractArray, sᵤ::AbstractArray, Gᵥ::AbstractMatrix, Nbodies::Integer)
    return - sᵤ + s[1:Nbodies*6] + Gᵥ'*s[1+Nbodies*6:end]
end

function g(mechanism::Mechanism)
    Ndims = sum([length(eqc) for eqc in mechanism.eqconstraints])
    g = zeros(MVector{Ndims})
    N = 1
    for eqc in mechanism.eqconstraints
        L = length(eqc)
        g[N:N+L-1] = ConstrainedDynamics.g(mechanism, eqc)
        N += L
    end
    return g
end

function resetMechanism!(mechanism::Mechanism, states::Vector{<:ConstrainedDynamics.State})
    for id in 1:length(mechanism.bodies)
        mechanism.bodies[id].state = deepcopy(states[id])
    end
    ConstrainedDynamics.discretizestate!(mechanism)
    foreach(ConstrainedDynamics.setsolution!, mechanism.bodies)
end

function projectv!(vᵤ::Vector{<:SVector}, ωᵤ::Vector{<:SVector}, mechanism::Mechanism; newtonIter::Integer = 100, ϵ::Real = 1e-10, regularizer::Real = 0.)
    Ndims = sum([length(eqc) for eqc in mechanism.eqconstraints])  # Total dimensionality of constraints
    Nbodies = length(mechanism.bodies)
    F = zeros(Nbodies*6 + Ndims, Nbodies*6 + Ndims)  # 3 vel, 3 ω for each body -> 6
    for i in 1:Nbodies*6 F[i,i] = 1 end
    s = zeros(MVector{6*Nbodies + Ndims})
    sᵤ = zeros(MVector{6*Nbodies})
    updateS!(s, vᵤ, ωᵤ)  # Initial s is [v1, ω1, v2, ω2, ...,  λ1, λ2, ...] with λ = 0
    updateS!(sᵤ, vᵤ, ωᵤ)
    updateMechanism!(mechanism, s)
    updateF!(F, mechanism)
    F += I*regularizer
    Gᵥ = (@view F[1+Nbodies*6:end, 1:Nbodies*6])
    f(s) = vcat(d(s, sᵤ, Gᵥ, Nbodies), g(mechanism))
    α = 1.
    for i in 1:newtonIter
        updateF!(F, mechanism)
        Δs = F\f(s)
        s -= α * Δs
        updateMechanism!(mechanism, s)
        if norm(f(s)) < ϵ && norm(Δs) < ϵ
            break
        end
    end
    return getvel(s, Nbodies)
end


text:
using RFPulse;

println("Reading RF Pulse file");
pulse = RFPulse.read("default.json");

println("Simulating...");
RFPulse.simulate!(pulse);

sim = pulse["SimulationParameters"]["Output"];

println("Test finished");


text:
"""
heqtest(n=100, c=.5, printh="yes"; jprecision=Float64, jmaxit=10)

Generates the examples for the SIREV paper.
"""
function heqtest(n = 100, c = 0.5, printh = "yes"; jprecision = Float64, jmaxit = 10)
    FS = ones(n)
    x0 = ones(n)
    hdata = heqinit(x0, n, c)
    FPS = ones(jprecision, n, n)
    nsoldoutfd = nsold01(
        x0,
        FS,
        FPS,
        heqf!;
        atol = 1.e-18,
        rtol = 1.e-18,
        pdata = hdata,
        maxit = jmaxit,
    )
    nsoldoutex = nsold01(
        x0,
        FS,
        FPS,
        heqf!,
        heqJ!;
        atol = 1.e-18,
        rtol = 1.e-18,
        pdata = hdata,
        maxit = jmaxit,
    )
    if printh == "yes"
        sol = nsoldoutex.solution
        cpout = chandprint(sol, hdata)
        for i = 1:21
            @printf("%1.2e      %1.5e \n", cpout[i, 1], cpout[i, 2])
        end
        @printf("\n \n \n")
        println(nsoldoutex.ithist)
    end
    return (exactout = nsoldoutex, fdout = nsoldoutfd)
end

max_stars_repo_name: ctkelley/MPResults


text:
# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
using HDF5
using Plots
using Printf
using LaTeXStrings
using Statistics
# %%
mode = "spherical" # "plane" or "spherical".
Hmin = 10
Hmax = 200 # maximum height from galactic plane in kpc
Hstep = 2
rlim = 10 # limit of radius (active only plane mode)
kmsinkpcyr = 3.1536e7/3.085677581e16 # 1 km/sec in kpc/yr
# timestep = np.linspace(0.01,0.8,80)
time = 100 # number of snapshots
tave = 50 # minimun snapshot number to use for calculating time average
snapinterval = 0.01 # interval of snapshot in Gyr
binwidthtemp = 0.15
rootdirectory = "/home/oku/SimulationData/isogal"
tempmin = 2
tempmax = 8
hotmin = 5 # Temperature to divide gas to Cold and Hot component
usePGFPlots = true

if !usePGFPlots
    using GR
    gr()
    extension = ".png"
    Plots.scalefontsizes(1.2)
    ENV["GKSwstype"] = "100"
else
    using PGFPlotsX
    pgfplotsx(tex_output_standalone = true)
    extension = ".tex"
    PGFPlotsX.enable_interactive(false)
end

print("initialization done. starting...\n")

H = Hmin:Hstep:Hmax
logT = tempmin:binwidthtemp:tempmax

models = [
    ["ss_model/ver10081/fiducial", "Fiducial", "green", "solid"],
    ["ss_model/ver10081/nonstochastic", "Non-stochastic", "blue", "dash"],
    ["ss_model/ver10081/sniikinetic", "SNII kinetic", "purple", "dot"],
    ["ss_model/ver10081/sniithermal", "SNII thermal", "red", "dashdot"],
    ["ss_model/ver10081/nofb", "No FB", "black", "dashdotdot"],
    # ["ss_model/ver10081/med", "High-reso", "red", "dash"],
     ]
# %%
# function W3(r, h)
#   r = abs(r)/h
#   C = 8/h^3/pi
#   if r > 1
#       return 0
#   elseif r > 1/2
#       return C*2*(1-r)^3
#   else
#       return C*(1 - 6*r^2 + 6*r^3)
#   end
# end
# function W5(r, h)
#     r =abs(r)/h
#     C = 3^9/359/pi/h^3
#     if r > 1
#         return 0
#     elseif r > 2/3
#         return C*(1 - r)^5
#     elseif r > 1/3
#         return C*((1 - r)^5 - 6*(2/3 - r)^5)
#     else
#         return C*((1 - r)^5 - 6*(2/3 - r)^5 + 15*(1/3 - r)^5)
#     end
# end
function crosssection(hsml, z)
    # function integrand(x)
    #     return W5(sqrt(z^2 + x^2),hsml)*2*pi*x
    #     # return W3(sqrt(z^2 + x^2),hsml)*2*pi*x
    # end
    # return quadgk(integrand, 0, sqrt(hsml^2 - z^2))[1]
    y = abs(z)
    C = 3^9/359/pi/hsml^3
    function f(a, b)
        return a*pi*hsml*C*(b - y/hsml)^6 * (y + hsml*(b - y/hsml)/7)
    end
    r = y/hsml
    if r > 1
        return 0
    elseif r > 2/3
        return f(1/3, 1)
    elseif r > 1/3
        return f(1/3, 1) + f(-2, 2/3)
    else
        return f(1/3, 1) + f(-2, 2/3) + f(5, 1/3)
    end
end
function getcenter(sph)
    x = y = z = mtot = 0
    for i in 1:length(sph.mass)
      x += sph.rho[i] * sph.mass[i] * sph.pos[1, i]
      y += sph.rho[i] * sph.mass[i] * sph.pos[2, i]
      z += sph.rho[i] * sph.mass[i] * sph.pos[3, i]
      mtot += sph.rho[i] * sph.mass[i]
    end
    x /= mtot
    y /= mtot
    z /= mtot
    return [x, y, z]
end
function getvelocity(sph)
    x = y = z = mtot = 0
    for i in 1:length(sph.mass)
      x += sph.mass[i] * sph.vel[1, i]
      y += sph.mass[i] * sph.vel[2, i]
      z += sph.mass[i] * sph.vel[3, i]
      mtot += sph.mass[i]
    end
    x /= mtot
    y /= mtot
    z /= mtot
    return [x, y, z]
end
function padding(array)
    ret = []
    for i in 1:length(array)
        value = array[i] > 0 ? array[i] : 1e-10
        push!(ret, value)
    end
    return ret
end
function timeaverage_dMdlogT(Results, h)
    dMdlogT = Dict()
    dMdlogT["mean"] = [[0.0 for j in 1:length(logT)] for i in 1:length(models)]
    dMdlogT["std"] = [[0.0 for j in 1:length(logT)] for i in 1:length(models)]
    for i in 1:length(models)
        for j in 1:length(logT)
            tmp = Array{Float32, 1}(undef, 0)
            for k in tave:time
                val = Results["dM/dlogT"][i][h][k][j]
                if val > 0
                    push!(tmp, log10(val))
                end
            end
            if length(tmp) > 25
                dMdlogT["mean"][i][j] = Statistics.mean(tmp)
                dMdlogT["std"][i][j] = Statistics.std(tmp, corrected=false)
            else
                dMdlogT["mean"][i][j] = -Inf
                dMdlogT["std"][i][j] = 0
            end
        end
    end
    return dMdlogT
end
function timeaverage_height(Results, key)
    ret = Dict()
    ret["mean"] = [Array{Float32, 1}(undef, 0) for i in 1:length(models)]
    ret["std"] = [Array{Float32, 1}(undef, 0) for i in 1:length(models)]
    ret["H"] = [Array{Float32, 1}(undef, 0) for i in 1:length(models)]
    for i in 1:length(models)
        for j in 1:length(H)
            tmp = Array{Float32, 1}(undef, 0)
            for k in tave:time
                val = Results[key][i][j][k]
                if val > 0
                    push!(tmp, log10(val))
                end
            end
            if length(tmp) > 25
                push!(ret["mean"][i], Statistics.mean(tmp))
                push!(ret["std"][i], Statistics.std(tmp, corrected=false))
                push!(ret["H"][i], H[j])
            end
        end
    end
    return ret
end

struct Sph
    pos::Array{Float32, 2}
    vel::Array{Float32, 2}
    hsml::Array{Float32, 1}
    mass::Array{Float32, 1}
    rho::Array{Float32, 1}
    temp::Array{Float32, 1}
    u::Array{Float32, 1}
    Z::Array{Float32, 1}
end
struct Galaxy
    pos::Array{Float32, 1}
    vel::Array{Float32, 1}
    sfr::Float32
    stellarmass::Float32
end

function main(H, sph, galaxy)
    ret = Dict()
    numbintemp = length(logT)
    ret["dM/dlogT"] = zeros(Int(numbintemp))
    ret["dM"] = 0.0
    ret["dE"] = 0.0
    ret["dMz"] = 0.0
    ret["dMCold"] = 0.0
    ret["dECold"] = 0.0
    ret["dMzCold"] = 0.0
    ret["dMHot"] = 0.0
    ret["dEHot"] = 0.0
    ret["dMzHot"] = 0.0
    MS = 0.0
    MSCold = 0.0
    MSHot = 0.0
    for i in 1:length(sph.mass)
        # initialization
        dz = sph.hsml[i] + 1.0
        z = 0

        if mode == "plane"
            z = abs(sph.pos[3,i] - galaxy.pos[3,1])
            dz = abs(z - H)
            r = sqrt((sph.pos[1,i]-galaxy.pos[1,1])^2 + (sph.pos[2,i]-galaxy.pos[2,1])^2)
            if r > rlim
                continue
            end
        elseif mode == "spherical"
            z = sqrt((sph.pos[1,i]-galaxy.pos[1,1])^2 + (sph.pos[2,i]-galaxy.pos[2,1])^2 + (sph.pos[3,i]-galaxy.pos[3,1])^2)
            dz = abs(z - H)
        else
            print("set mode plane or spherical")
            break
        end

        if dz < sph.hsml[i] && sph.Z[i] > 0
            v = -1 #initialization.
            if mode == "plane"
                v = (sph.vel[3,i] - galaxy.vel[3,1])*(sph.pos[3,i] - galaxy.pos[3,1])/z
            elseif mode == "spherical"
                v =
                 ((sph.vel[1,i]-galaxy.vel[1,1])*(sph.pos[1,i]-galaxy.pos[1,1])
                + (sph.vel[2,i]-galaxy.vel[2,1])*(sph.pos[2,i]-galaxy.pos[2,1])
                + (sph.vel[3,i]-galaxy.vel[3,1])*(sph.pos[3,i]-galaxy.pos[3,1]))/z
            end
            if v > 0
                S = crosssection(sph.hsml[i], dz)

                dM = sph.mass[i]*1e10*S*v*kmsinkpcyr
                ret["dM"] += dM

                dMS = sph.mass[i]*1e10*S*kmsinkpcyr
                MS += dMS

                energy = sph.mass[i]*(0.5*v^2 + sph.u[i])
                dE = energy*2e53*S*v*kmsinkpcyr # 2e53 is a conversion factor from 10^10 Msun km^2 s^-2 to erg
                ret["dE"] += dE

                dMz = sph.Z[i]*dM
                ret["dMz"] += dMz

                if log10(sph.temp[i]) < hotmin
                    ret["dMCold"] += dM
                    ret["dECold"] += dE
                    ret["dMzCold"] += dMz
                    MSCold += dMS
                else
                    ret["dMHot"] += dM
                    ret["dEHot"] += dE
                    ret["dMzHot"] += dMz
                    MSHot += dMS
                end

                if log10(sph.temp[i]) > tempmin && log10(sph.temp[i]) < tempmax
                    ibintemp = Int(ceil((log10(sph.temp[i]) - tempmin)/binwidthtemp))
                    ret["dM/dlogT"][ibintemp] += dM/binwidthtemp
                end
            end
        end
    end
    ret["Z"] = ret["dMz"]/ret["dM"]
    ret["ZCold"] = ret["dMzCold"]/ret["dMCold"]
    ret["ZHot"] = ret["dMzHot"]/ret["dMHot"]

    ret["V"] = ret["dM"]/MS
    ret["VCold"] = ret["dMCold"]/MSCold
    ret["VHot"] = ret["dMHot"]/MSHot

    return ret
end

## main part
Results = Dict()
Results["dM"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["dE"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["dMz"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["Z"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["MassLoadingFactor"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["EnergyLoadingFactor"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["MetalLoadingFactor"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["V"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]

Results["dMCold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["dECold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["dMzCold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["ZCold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["MassLoadingFactorCold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["EnergyLoadingFactorCold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["MetalLoadingFactorCold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["VCold"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]

Results["dMHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["dEHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["dMzHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["ZHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["MassLoadingFactorHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["EnergyLoadingFactorHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["MetalLoadingFactorHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]
Results["VHot"] = [[[0.0 for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]

Results["dM/dlogT"] = [[[[] for k in 1:time] for j in 1:length(H)] for i in 1:length(models)]

Results["SFR"] = [[0.0 for k in 1:time] for i in 1:length(models)]
Results["StellarMass"] = [[0.0 for k in 1:time] for i in 1:length(models)]
Results["MetalMass"] = [[0.0 for k in 1:time] for i in 1:length(models)]
for i in 1:length(models)
    @time for j in 1:time
        ## load data
        path2snapshot = @sprintf("%s/%s/snapshot_%03d/snapshot_%03d.hdf5", rootdirectory, models[i][1], j, j)
        sfr = 0.0
        massstars = 0.0
        h5open(path2snapshot, "r") do fp
            pos = read(fp,"PartType0/Coordinates")
            vel = read(fp,"PartType0/Velocities")
            hsml = read(fp, "PartType0/SmoothingLength")
            mass = read(fp, "PartType0/Masses")
            rho = read(fp, "PartType0/Density")
            temp = read(fp, "PartType0/Temperature")
            u = read(fp, "PartType0/InternalEnergy")
            Z = read(fp, "PartType0/Metallicity")
            sfr = read(fp, "PartType0/StarFormationRate")
            massstars = read(fp, "PartType4/Masses")
            global sph = Sph(pos, vel, hsml, mass, rho, temp, u, Z)
        end
        galpos = getcenter(sph)
        galvel = getvelocity(sph)
        galsfr = sum(sfr)
        stellarmass = sum(massstars)
        global galaxy = Galaxy(galpos, galvel, galsfr, stellarmass)

        for k in 1:length(H)
            h = H[k]
            ## main calculation
            ret = main(h, sph, galaxy)

            ## save
            for l in collect(keys(ret))
                Results[l][i][k][j] = ret[l]
            end

            Results["MassLoadingFactor"][i][k][j] = ret["dM"]/galaxy.sfr
            Results["EnergyLoadingFactor"][i][k][j] = ret["dE"]/(7.2e48 * galaxy.sfr)
            Results["MetalLoadingFactor"][i][k][j] = ret["dMz"]/(1e-2 * galaxy.sfr)

            Results["MassLoadingFactorCold"][i][k][j] = ret["dMCold"]/galaxy.sfr
            Results["EnergyLoadingFactorCold"][i][k][j] = ret["dECold"]/(7.2e48 * galaxy.sfr)
            Results["MetalLoadingFactorCold"][i][k][j] = ret["dMzCold"]/(1e-2 * galaxy.sfr)

            Results["MassLoadingFactorHot"][i][k][j] = ret["dMHot"]/galaxy.sfr
            Results["EnergyLoadingFactorHot"][i][k][j] = ret["dEHot"]/(7.2e48 * galaxy.sfr)
            Results["MetalLoadingFactorHot"][i][k][j] = ret["dMzHot"]/(1e-2 * galaxy.sfr)
        end
        Results["SFR"][i][j] = galaxy.sfr
        Results["StellarMass"][i][j] = galaxy.stellarmass
        for k in 1:length(sph.mass)
            Results["MetalMass"][i][j] += sph.mass[k]*sph.Z[k]
        end
    end
end

# %%

X = snapinterval:snapinterval:time*snapinterval

# galaxy properties
Plots.plot(ylabel="SFR [Msun/yr]", xlabel="Time [Gyr]", ylim=(0,16))
for i in 1:length(models)
Plots.plot!(X, Results["SFR"][i], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/SFR"*extension)

Plots.plot()
for i in 1:length(models)
Plots.plot!(X, Results["StellarMass"][i], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/StellarMass"*extension)

Plots.plot()
for i in 1:length(models)
Plots.plot!(X, Results["MetalMass"][i], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/MetalMass"*extension)

h = 1
x = H[h]
print("properties at H = $x\n")
dMdlogT = timeaverage_dMdlogT(Results, h)
Plots.plot(ylim=(-4,2), xlim=(2,8), xlabel=latexstring("\$\\log\\,(T\\,\\,[\\mathrm{K}])\$"), ylabel=latexstring("\$\\log\\,(d\\dot{M}_{\\mathrm{out}}/d\\log T\\,\\,[M_{\\odot} \\mathrm{yr}^{-1} \\mathrm{dex}^{-1}])\$"))
for i in 1:length(models)
Plots.plot!(logT, dMdlogT["mean"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]), ribbon=dMdlogT["std"][i], fillalpha=.3)
end
Plots.savefig("results/dMdlogT"*extension)

# time evolution of outflow properties
Plots.plot(ylabel="Mass outflow rate [Msun/yr]", xlabel="Time [Gyr]",legend=:topright, yscale=:log10, ylim=(1e-3,2e3))
for i in 1:length(models)
Plots.plot!(X, Results["dM"][i][h], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/MassOutFlowRate"*extension)

Plots.plot(xlabel="Time [Gyr]", ylabel="Mass loading factor", legend=:topleft, yscale=:log10, ylim=(5e-4,1e3))
for i in 1:length(models)
Plots.plot!(X, Results["MassLoadingFactor"][i][h], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/MassLoadingFactor"*extension)

Plots.plot(ylabel="Energy outflow rate [erg/yr]", xlabel="Time [Gyr]",legend=:bottomright, yscale=:log10, ylim=(1e44,1e49))
for i in 1:length(models)
Plots.plot!(X, Results["dE"][i][h], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/EnergyOutFlowRate"*extension)

Plots.plot(xlabel="Time [Gyr]", ylabel="Energy loading factor", legend=:topleft, yscale=:log10, ylim=(1e-6,2e0))
for i in 1:length(models)
Plots.plot!(X, Results["EnergyLoadingFactor"][i][h], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/EnergyLoadingFactor"*extension)

Plots.plot(ylabel="Metal mass outflow rate [Msun/yr]", xlabel="Time [Gyr]",legend=:bottomright, yscale=:log10, ylim=(1e-7, 2e1))
for i in 1:length(models)
Plots.plot!(X, Results["dMz"][i][h], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/MetalOutFlowRate"*extension)

Plots.plot(xlabel="Time [Gyr]", ylabel="Metal loading factor", legend=:topleft, yscale=:log10, ylim=(1e-3,1e1))
for i in 1:length(models)
Plots.plot!(X, Results["MetalLoadingFactor"][i][h], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/MetalLoadingFactor"*extension)

Plots.plot(xlabel="Time [Gyr]", ylabel="Metallicity", legend=:topleft, yscale=:log10, ylim=(1e-3,1e-1))
for i in 1:length(models)
Plots.plot!(X, Results["Z"][i][h], label=models[i][2], linecolor=Symbol(models[i][3]), linestyle=Symbol(models[i][4]))
end
Plots.savefig("results/Metallicity"*extension)

# outflow profile against height from the galaxy

XLABEL_H = ""
XSCALE = ""
if mode == "spherical"
    # XSCALE = Symbol("log10")
    XSCALE = Symbol("none")
    XLABEL_H = latexstring("\$R\\,\\mathrm{[kpc]}\$")
elseif mode == "plane"
    XSCALE = Symbol("none")
    XLABEL_H = latexstring("\$H\\,\\mathrm{[kpc]}\$")
end

AllPhase = ["MassLoadingFactor", "EnergyLoadingFactor", "MetalLoadingFactor", "V", "Z"]
ColdPhase = ["MassLoadingFactorCold", "EnergyLoadingFactorCold", "MetalLoadingFactorCold", "VCold", "ZCold"]
HotPhase = ["MassLoadingFactorHot", "EnergyLoadingFactorHot", "MetalLoadingFactorHot", "VHot", "ZHot"]
function addannotate(key, xmin, ymax)
    if key in AllPhase
        Plots.annotate!([(xmin, ymax-0.2, (latexstring("\$\\mathrm{All}\$"), 16, :black, :left))])
    elseif key in ColdPhase
        Plots.annotate!([(xmin, ymax-0.2, (latexstring("\$\\mathrm{Cold}\\,(T < 10^5\\,\\mathrm{K})\$"), 16, :black, :left))])
    elseif key in HotPhase
        Plots.annotate!([(xmin, ymax-0.2, (latexstring("\$\\mathrm{Hot}\\,(T > 10^5\\,\\mathrm{K})\$"), 16, :black, :left))])
    end
end

for key in ["MassLoadingFactor", "MassLoadingFactorCold", "MassLoadingFactorHot"]
    eta = timeaverage_height(Results, key)
    ymin = -4
    ymax = 2
    Plots.plot(xlabel=XLABEL_H, xscale=XSCALE, ylabel=latexstring("\$\\log\\,\\eta_m \$"), legend=:topright, ylim=(ymin,ymax))
    addannotate(key, Hstep, ymax)
    for i in 1:length(models)
        Plots.plot!(eta["H"][i], eta["mean"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]), ribbon=eta["std"][i], fillalpha=.3)
    end
    Plots.savefig("results/"*key*"_H"*extension)
end

for key in ["EnergyLoadingFactor", "EnergyLoadingFactorCold", "EnergyLoadingFactorHot"]
    eta = timeaverage_height(Results, key)
    ymin = -6
    ymax = 0
    Plots.plot(xlabel=XLABEL_H, xscale=XSCALE, ylabel=latexstring("\$\\log\\,\\eta_e \$"), legend=:topright, ylim=(ymin, ymax))
    addannotate(key, Hstep, ymax)
    for i in 1:length(models)
        Plots.plot!(eta["H"][i], eta["mean"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]), ribbon=eta["std"][i], fillalpha=.3)
    end
    Plots.savefig("results/"*key*"_H"*extension)
end

for key in ["MetalLoadingFactor", "MetalLoadingFactorCold", "MetalLoadingFactorHot"]
    eta = timeaverage_height(Results, key)
    ymin = -3
    ymax = 1
    Plots.plot(xlabel=XLABEL_H, xscale=XSCALE, ylabel=latexstring("\$\\log\\,\\eta_Z \$"), legend=:topright, ylim=(ymin, ymax))
    addannotate(key, Hstep, ymax)
    for i in 1:length(models)
        Plots.plot!(eta["H"][i], eta["mean"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]), ribbon=eta["std"][i], fillalpha=.3)
    end
    Plots.savefig("results/"*key*"_H"*extension)
end

for key in ["V", "VCold", "VHot"]
    V = timeaverage_height(Results, key)
    ymin = 0
    ymax = 3
    Plots.plot(xlabel=XLABEL_H, xscale=XSCALE, ylabel=latexstring("\$\\log\\,(V \\mathrm{[km}\\,\\mathrm{s}^{-1}\\mathrm{]}\$"), legend=:bottomright, ylim=(ymin, ymax))
    addannotate(key, Hstep, ymax)
    for i in 1:length(models)
        Plots.plot!(V["H"][i], V["mean"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]), ribbon=V["std"][i], fillalpha=.3)
    end
    Plots.savefig("results/"*key*"_H"*extension)
end

for key in ["Z", "ZCold", "ZHot"]
    Z = timeaverage_height(Results, key)
    ymin = -3
    ymax = -1
    Plots.plot(xlabel=XLABEL_H, xscale=XSCALE, ylabel=latexstring("\$\\log\\,Z\$"), legend=:topright, ylim=(ymin, ymax))
    addannotate(key, Hstep, ymax)
    for i in 1:length(models)
        Plots.plot!(Z["H"][i], Z["mean"][i], label=models[i][2], color=Symbol(models[i][3]), linestyle=Symbol(models[i][4]), ribbon=Z["std"][i], fillalpha=.3)
    end
    Plots.savefig("results/"*key*"_H"*extension)
end


text:
using LIKWID
using Base.Threads
using Test

if length(ARGS) < 1
    @info("test_pin.jl: No cpustr provided to `test_pin.jl`. Assuming first N cores.")
    cores = 0:nthreads()-1
else
    cpustr = ARGS[1]
    if contains(cpustr, "-")
        s = split(cpustr, "-")
        start = parse(Int, s[1])
        stop = parse(Int, s[2])
        cores = start:stop
    else
        # assume explicit list of cores, e.g. "0,4,2,3,12,1"
        cores = parse.(Int, split(cpustr, ","))
    end
end

threadids = zeros(Int, nthreads())
procids = zeros(Int, nthreads())
procids_glibc = zeros(Int, nthreads())
@threads :static for i in 1:nthreads()
    threadids[i] = threadid()
    procids[i] = LIKWID.get_processor_id()
    procids_glibc[i] = LIKWID.get_processor_id_glibc()
end

@test threadids == 1:nthreads()
@test procids == procids_glibc == cores


text:
abstract RegressionFunction

# Linear predictive function

type LinearRegressionFunction <: RegressionFunction
end

values(::LinearRegressionFunction, X::AbstractMatrix, w::Vector) = X*w
gradient(::LinearRegressionFunction, X::AbstractMatrix, ::Vector) = X

# TODO: Dual predictive function


text:
function gradient2d(M::Matrix{Float64}, mask::Matrix{Int})
    di, dj = zeros(size(M)), zeros(size(M))

    nr = size(M, 1)
    nc = size(M, 2)
    v = zeros(3)

    p = zeros(3, 2)
    p[1,:] = [1 0]; v[1] = M[2,1] - M[1,1]
    p[2,:] = [0 1]; v[2] = M[1,2] - M[1,1]
    p[3,:] = [1 1]; v[3] = M[2,2] - M[1,1]

    n = p \ v

    di[1,1] = n[1]
    dj[1,1] = n[2]

    p[1,:] = [0 -1]; v[1] =  M[1,nc-1] - M[1,nc]
    p[2,:] = [1 -1]; v[2] =  M[2,nc-1] - M[1,nc]
    p[3,:] = [1 0];  v[3] =  M[2,nc]   - M[1,nc]

    n = p \ v

    di[1,nc] = n[1]
    dj[1,nc] = n[2]

    p[1,:] = [-1 0]; v[1] = M[nr-1,1] - M[nr,1]
    p[2,:] = [-1 1]; v[2] = M[nr-1,2] - M[nr,1]
    p[3,:] = [0 1];  v[3] = M[nr,2]   - M[nr,1]

    n = p \ v

    di[nr,1] = n[1]
    dj[nr,1] = n[2]

    p[1,:] = [-1 -1]; v[1] = M[nr-1,nc-1] - M[nr,nc]
    p[2,:] = [0 -1];  v[2] = M[nr,nc-1]   - M[nr,nc]
    p[3,:] = [-1 0];  v[3] = M[nr-1,nc]   - M[nr,nc]

    n = p \ v

    di[nr,nc] = n[1]
    dj[nr,nc] = n[2]

    # then, the boundaries

    v = zeros(5)
    p = zeros(5, 2)

    for j=2:nc-1
        p[1,:] = [0 -1]; v[1] = M[1,j-1] - M[1,j]
        p[2,:] = [1 -1]; v[2] = M[2,j-1] - M[1,j]
        p[3,:] = [1 0];  v[3] = M[2,j]   - M[1,j]
        p[4,:] = [0 1];  v[4] = M[1,j+1] - M[1,j]
        p[5,:] = [1 1];  v[5] = M[2,j+1] - M[1,j]

        n = p \ v

        di[1,j] = n[1]
        dj[1,j] = n[2]

        p[1,:] = [-1 -1]; v[1] = M[nr-1,j-1] - M[nr,j]
        p[2,:] = [0 -1];  v[2] = M[nr,j-1]   - M[nr,j]
        p[3,:] = [-1 0];  v[3] = M[nr-1,j]   - M[nr,j]
        p[4,:] = [-1 1];  v[4] = M[nr-1,j+1] - M[nr,j]
        p[5,:] = [0 1];   v[5] = M[nr,j+1]   - M[nr,j]

        n = p \ v

        di[nr,j] = n[1]
        dj[nr,j] = n[2]
    end

    for i=2:nr-1
        p[1,:] = [-1 0]; v[1] = M[i-1,1] - M[i,1]
        p[2,:] = [1 0];  v[2] = M[i+1,1] - M[i,1]
        p[3,:] = [-1 1]; v[3] = M[i-1,2] - M[i,1]
        p[4,:] = [0 1];  v[4] = M[i,2]   - M[i,1]
        p[5,:] = [1 1];  v[5] = M[i+1,2] - M[i,1]

        n = p \ v

        di[i,1] = n[1]
        dj[i,1] = n[2]

        p[1,:] = [-1 -1]; v[1] = M[i-1,nc-1] - M[i,nc]
        p[2,:] = [0 -1];  v[2] = M[i,nc-1]   - M[i,nc]
        p[3,:] = [1 -1];  v[3] = M[i+1,nc-1] - M[i,nc]
        p[4,:] = [-1 0];  v[4] = M[i-1,nc]   - M[i,nc]
        p[5,:] = [1 0];   v[5] = M[i+1,nc]   - M[i,nc]

        n = p \ v

        di[i,nc] = n[1]
        dj[i,nc] = n[2]
    end

    # and finally the remaining function values
    p = zeros(8, 2)
    v = zeros(8)
    for i=2:nr-1
        for j=2:nc-1
            if mask[i,j] != 0
                k = 1

                if mask[i-1,j-1] != 0
                    p[k,:] = [-1 -1]; v[k] = M[i-1,j-1] - M[i,j]
                    k = k + 1
                end

                if mask[i-1,j] != 0
                    p[k,:] = [-1 0];  v[k] = M[i-1,j]   - M[i,j]
                    k = k + 1
                end

                if mask[i-1,j+1] != 0
                    p[k,:] = [-1 1];  v[k] = M[i-1,j+1] - M[i,j]
                    k = k + 1
                end

                if mask[i,j-1] != 0
                    p[k,:] = [0 -1];  v[k] = M[i,j-1]   - M[i,j]
                    k = k + 1
                end

                if mask[i,j+1] != 0
                    p[k,:] = [0 1];   v[k] = M[i,j+1]   - M[i,j]
                    k = k + 1
                end

                if mask[i+1,j-1] != 0
                    p[k,:] = [1 -1];  v[k] = M[i+1,j-1] - M[i,j]
                    k = k + 1
                end

                if mask[i+1,j] != 0
                    p[k,:] = [1 0];   v[k] = M[i+1,j]   - M[i,j]
                    k = k + 1
                end

                if mask[i+1,j+1] != 0
                    p[k,:] = [1 1];   v[k] = M[i+1,j+1] - M[i,j]
                    k = k + 1
                end

                if k > 2
                    n = p[1:k-1,:] \ v[1:k-1]

                    di[i,j] = n[1]
                    dj[i,j] = n[2]
                else
                    di[i,j] = 0
                    dj[i,j] = 0
                end
            else
                di[i,j] = 0
                dj[i,j] = 0
            end
        end
    end
    return Grad2D(di, dj)
end


text:
// wypisz liczby parzyste do 10

int main () {
  int i ;
  i = 0 ;
  while (i < 10){
    if (i % 2 == 0) printInt(i) ;
    i++ ;
    } ;
  printInt(i) ;
  return 0 ;
  } ;

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
using BSplines
using Test

@testset "BSplines.jl" begin
    # Write your tests here.
end


text:
function is_print_time(tm, dt_print)
  if mod(tm, dt_print) < dtm
      is_print = true
      println("--writing files, t = ", tm, "; maxTime = ", maxTime, "; dtm = ", dtm)
  else
      is_print = false
  end
  return is_print
end

function is_quiet_time!(time1, time2, tm, tm_quiet, geom, nz, nq, dt_print, dtm)
    if (time1 && tm < tm_quiet)
        if geom=="slab"
            nq = nquiet
        elseif geom=="spherical"
            nq = nz - nquiet
        end
        println("nz = ", nz, " nquiet = ", nquiet, " tm_quiet=", tm_quiet)
        println("dt_print = ", dt_print, " dtm = ", dtm)
        time1 = false
    elseif (time2 && tm >= tm_quiet)
        println("end of quiet time")
        nq = nz
        # dtm = dtm / dt_multi
        # lm = dtm / dr
        time2 = false
    end
    return nq, time1, time2
end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
module TestAqua

using Aqua
using StopWatches
using StopWatchesBase

function test_stopwatches()
    Aqua.test_all(StopWatches; ambiguities = false)
end

function test_stopwatchesbase()
    Aqua.test_all(StopWatchesBase)
end

end  # module


text:
using PyPlot
using Random

ntrials = 10000            # Number of experiments
x = Int64[]                # Array of outcomes
sumx = 0
for i = 1:ntrials
    nrolls = 0
    found_a_six = false
    while !found_a_six
        nrolls += 1
        roll = rand(1:6)
        if roll == 6
            found_a_six = true
        end
    end
    push!(x, nrolls)
    sumx += nrolls
end

plt.hist(x, 50);
xlabel("Number of rolls")
ylabel("Count")
average = sumx / ntrials

function card_rank(card)
    (card - 1) % 13 + 1
end

function card_suit(card)
    (card - 1) ÷ 13 + 1
end

cards = 1:52      # All cards
[card_rank.(cards) card_suit.(cards)]

cards = randperm(52)
hand = cards[1:5]
[card_rank.(hand) card_suit.(hand)]

ntrials = 100000
nflush = 0
for itrial = 1:ntrials
    cards = randperm(52)
    hand = cards[1:5]

    suits = card_suit.(hand)
    same_suit = true
    for i = 2:5
        if suits[i] ≠ suits[1]
            same_suit = false
            break
        end
    end
    if same_suit
        nflush += 1
    end
end

approx_probability = nflush / ntrials

exact_probability = 4 * binomial(13,5) / binomial(52,5)

n = 10000               # Number of darts
x = 2rand(n) .- 1       # Coordinates
y = 2rand(n) .- 1

plot([-1,1,1,-1,-1], [-1,-1,1,1,-1])  # Draw square
theta = 2π*(0:100)./100               # Draw circle
plot(cos.(theta), sin.(theta))
plot(x, y, linestyle="None", marker=".", markersize=1)   # Plot dart points
axis("equal")

# Determine if points are inside the circle (a "hit")
hits = 0
for i = 1:n
    if x[i]^2 + y[i]^2 ≤ 1
        hits += 1
    end
end

approx_pi = 4hits / n

function random_walk(n)
    x = [0]
    y = [0]
    while abs(x[end]) < n && abs(y[end]) < n
        if rand() < 0.5
            if rand() < 0.5  # Up
                push!(x, x[end])
                push!(y, y[end] + 1)
            else             # Down
                push!(x, x[end])
                push!(y, y[end] - 1)
            end
        else
            if rand() < 0.5  # Right
                push!(x, x[end] + 1)
                push!(y, y[end])
            else             # Left
                push!(x, x[end] - 1)
                push!(y, y[end])
            end
        end
    end
    x,y
end

n = 10
x,y = random_walk(n)

plot(x, y, marker=".", markersize=8)     # Draw dots at each point in x,y
grid(true)
axis([-n,n,-n,n]);


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

# connected component of adjacency matrix containing vertex
function component(adjacency::AbstractMatrix{Int}, vertex::Int)
  frontier = [vertex]
  visited  = Int[]
  # breadth-first search
  while !isempty(frontier)
    u = pop!(frontier)
    push!(visited, u)
    for v in findall(!iszero, adjacency[u,:])
      if v ∉ visited
        push!(frontier, v)
      end
    end
  end
  visited
end

"""
    @metasolver solver solvertype body

A helper macro to create a solver named `solver` of type `solvertype`
with parameters specified in `body`.

## Examples

Create a solver with parameters `mean` and `variogram` for each variable
of the problem, and a global parameter that specifies whether or not
to use the GPU:

```julia
julia> @metasolver MySolver AbstractSimulationSolver begin
  @param mean = 0.0
  @param variogram = GaussianVariogram()
  @jparam rho = 0.7
  @global gpu = false
end
```

### Notes

This macro is not intended to be used directly, see other macros defined
below for estimation and simulation solvers.
"""
macro metasolver(solver, solvertype, body)
  # discard any content that doesn't start with @param or @global
  content = filter(arg -> arg isa Expr, body.args)

  # lines starting with @param refer to single variable parameters
  vparams = filter(p -> p.args[1] == Symbol("@param"), content)
  vparams = map(p -> p.args[3], vparams)

  # lines starting with @jparam refer to joint variable parameters
  jparams = filter(p -> p.args[1] == Symbol("@jparam"), content)
  jparams = map(p -> p.args[3], jparams)

  # lines starting with @global refer to global solver parameters
  gparams = filter(p -> p.args[1] == Symbol("@global"), content)
  gparams = map(p -> p.args[3], gparams)

  # add default value of `nothing` if necessary
  gparams = map(p -> p isa Symbol ? :($p = nothing) : p, gparams)

  # replace Expr(:=, a, 2) by Expr(:kw, a, 2) for valid kw args
  gparams = map(p -> Expr(:kw, p.args...), gparams)

  # keyword names
  gkeys = map(p -> p.args[1], gparams)

  # solver parameter type for single variable
  solvervparam = Symbol(solver,"Param")

  # solver parameter type for joint variables
  solverjparam = Symbol(solver,"JointParam")

  # variables are symbols or tuples of symbols
  vtype = Symbol
  jtype = NTuple{<:Any,Symbol}

  esc(quote
    $Parameters.@with_kw_noshow struct $solvervparam
      __dummy__ = nothing
      $(vparams...)
    end

    $Parameters.@with_kw_noshow struct $solverjparam
      __dummy__ = nothing
      $(jparams...)
    end

    @doc (@doc $solvervparam) (
    struct $solver <: $solvertype
      vparams::Dict{$vtype,$solvervparam}
      jparams::Dict{$jtype,$solverjparam}
      $(gkeys...)

      # auxiliary fields
      varnames::Vector{Symbol}
      adjacency::Matrix{Int}

      function $solver(vparams::Dict{$vtype,$solvervparam},
                       jparams::Dict{$jtype,$solverjparam},
                       $(gkeys...))
        svars = collect(keys(vparams))
        jvars = collect(keys(jparams))
        lens₁ = length.(jvars)
        lens₂ = length.(unique.(jvars))

        @assert all(lens₁ .== lens₂ .> 1) "invalid joint variable specification"

        varnames = svars ∪ Iterators.flatten(jvars)

        adjacency = zeros(Int, length(varnames), length(varnames))
        for (i, u) in enumerate(varnames)
          for vtuple in jvars
            if u ∈ vtuple
              for v in vtuple
                j = indexin([v], varnames)[1]
                i == j || (adjacency[i,j] = 1)
              end
            end
          end
        end

        new(vparams, jparams, $(gkeys...), varnames, adjacency)
      end
    end)

    function $solver(params...; $(gparams...))
      # build dictionaries for inner constructor
      vdict = Dict{$vtype,$solvervparam}()
      jdict = Dict{$jtype,$solverjparam}()

      # convert named tuples to solver parameters
      for (varname, varparams) in params
        kwargs = [k => v for (k,v) in zip(keys(varparams), varparams)]
        if varname isa Symbol
          push!(vdict, varname => $solvervparam(; kwargs...))
        else
          push!(jdict, varname => $solverjparam(; kwargs...))
        end
      end

      $solver(vdict, jdict, $(gkeys...))
    end

    function GeoStatsBase.covariables(var::Symbol, solver::$solver)
      vind = indexin([var], solver.varnames)[1]
      if vind ≠ nothing
        comp = GeoStatsBase.component(solver.adjacency, vind)
        vars = Tuple(solver.varnames[sort(comp)])
        params = []
        for v in vars
          push!(params, (v,) => solver.vparams[v])
        end
        for vtuple in keys(solver.jparams)
          if any(v ∈ vars for v in vtuple)
            push!(params, vtuple => solver.jparams[vtuple])
          end
        end
      else
        # default parameter for single variable
        vars = (var,)
        params = [(var,) => $solvervparam()]
      end

      (names=vars, params=Dict(params))
    end

    GeoStatsBase.variables(solver::$solver) = solver.varnames

    # ------------
    # IO methods
    # ------------
    function Base.show(io::IO, solver::$solver)
      print(io, $solver)
    end

    function Base.show(io::IO, ::MIME"text/plain", solver::$solver)
      println(io, solver)
      for (var, varparams) in merge(solver.vparams, solver.jparams)
        if var isa Symbol
          println(io, "  └─", var)
        else
          println(io, "  └─", join(var, "—"))
        end
        pnames = setdiff(fieldnames(typeof(varparams)), [:__dummy__])
        for pname in pnames
          pval = getfield(varparams, pname)
          if pval ≠ nothing
            print(io, "    └─", pname, " ⇨ ")
            show(IOContext(io, :compact => true), pval)
            println(io, "")
          end
        end
      end
    end
  end)
end

"""
    @estimsolver solver body

A helper macro to create a estimation solver named `solver` with parameters
specified in `body`. For examples, please check the documentation for
`@metasolver`.
"""
macro estimsolver(solver, body)
  esc(quote
    GeoStatsBase.@metasolver $solver GeoStatsBase.AbstractEstimationSolver $body
  end)
end

"""
    @estimsolver solver body

A helper macro to create a simulation solver named `solver` with parameters
specified in `body`. For examples, please check the documentation for
`@metasolver`.
"""
macro simsolver(solver, body)
  esc(quote
    GeoStatsBase.@metasolver $solver GeoStatsBase.AbstractSimulationSolver $body
  end)
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ProtoBuf.google.protobuf

mutable struct CollectionDef_NodeList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CollectionDef_NodeList(; kwargs...)
        obj = new(meta(CollectionDef_NodeList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CollectionDef_NodeList
const __meta_CollectionDef_NodeList = Ref{ProtoMeta}()
function meta(::Type{CollectionDef_NodeList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CollectionDef_NodeList)
            __meta_CollectionDef_NodeList[] = target = ProtoMeta(CollectionDef_NodeList)
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{AbstractString}]
            meta(target, CollectionDef_NodeList, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CollectionDef_NodeList[]
    end
end
function Base.getproperty(obj::CollectionDef_NodeList, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    else
        getfield(obj, name)
    end
end

mutable struct CollectionDef_BytesList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CollectionDef_BytesList(; kwargs...)
        obj = new(meta(CollectionDef_BytesList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CollectionDef_BytesList
const __meta_CollectionDef_BytesList = Ref{ProtoMeta}()
function meta(::Type{CollectionDef_BytesList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CollectionDef_BytesList)
            __meta_CollectionDef_BytesList[] = target = ProtoMeta(CollectionDef_BytesList)
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{Vector{UInt8}}]
            meta(target, CollectionDef_BytesList, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CollectionDef_BytesList[]
    end
end
function Base.getproperty(obj::CollectionDef_BytesList, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Vector{UInt8}}
    else
        getfield(obj, name)
    end
end

mutable struct CollectionDef_Int64List <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CollectionDef_Int64List(; kwargs...)
        obj = new(meta(CollectionDef_Int64List), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CollectionDef_Int64List
const __meta_CollectionDef_Int64List = Ref{ProtoMeta}()
function meta(::Type{CollectionDef_Int64List})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CollectionDef_Int64List)
            __meta_CollectionDef_Int64List[] = target = ProtoMeta(CollectionDef_Int64List)
            pack = Symbol[:value]
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{Int64}]
            meta(target, CollectionDef_Int64List, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CollectionDef_Int64List[]
    end
end
function Base.getproperty(obj::CollectionDef_Int64List, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    else
        getfield(obj, name)
    end
end

mutable struct CollectionDef_FloatList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CollectionDef_FloatList(; kwargs...)
        obj = new(meta(CollectionDef_FloatList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CollectionDef_FloatList
const __meta_CollectionDef_FloatList = Ref{ProtoMeta}()
function meta(::Type{CollectionDef_FloatList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CollectionDef_FloatList)
            __meta_CollectionDef_FloatList[] = target = ProtoMeta(CollectionDef_FloatList)
            pack = Symbol[:value]
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{Float32}]
            meta(target, CollectionDef_FloatList, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CollectionDef_FloatList[]
    end
end
function Base.getproperty(obj::CollectionDef_FloatList, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float32}
    else
        getfield(obj, name)
    end
end

mutable struct CollectionDef_AnyList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CollectionDef_AnyList(; kwargs...)
        obj = new(meta(CollectionDef_AnyList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CollectionDef_AnyList
const __meta_CollectionDef_AnyList = Ref{ProtoMeta}()
function meta(::Type{CollectionDef_AnyList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CollectionDef_AnyList)
            __meta_CollectionDef_AnyList[] = target = ProtoMeta(CollectionDef_AnyList)
            allflds = Pair{Symbol,Union{Type,String}}[:value => Base.Vector{ProtoBuf.google.protobuf._Any}]
            meta(target, CollectionDef_AnyList, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CollectionDef_AnyList[]
    end
end
function Base.getproperty(obj::CollectionDef_AnyList, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ProtoBuf.google.protobuf._Any}
    else
        getfield(obj, name)
    end
end

mutable struct CollectionDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CollectionDef(; kwargs...)
        obj = new(meta(CollectionDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CollectionDef
const __meta_CollectionDef = Ref{ProtoMeta}()
function meta(::Type{CollectionDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CollectionDef)
            __meta_CollectionDef[] = target = ProtoMeta(CollectionDef)
            allflds = Pair{Symbol,Union{Type,String}}[:node_list => CollectionDef_NodeList, :bytes_list => CollectionDef_BytesList, :int64_list => CollectionDef_Int64List, :float_list => CollectionDef_FloatList, :any_list => CollectionDef_AnyList]
            oneofs = Int[1,1,1,1,1]
            oneof_names = Symbol[Symbol("kind")]
            meta(target, CollectionDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_CollectionDef[]
    end
end
function Base.getproperty(obj::CollectionDef, name::Symbol)
    if name === :node_list
        return (obj.__protobuf_jl_internal_values[name])::CollectionDef_NodeList
    elseif name === :bytes_list
        return (obj.__protobuf_jl_internal_values[name])::CollectionDef_BytesList
    elseif name === :int64_list
        return (obj.__protobuf_jl_internal_values[name])::CollectionDef_Int64List
    elseif name === :float_list
        return (obj.__protobuf_jl_internal_values[name])::CollectionDef_FloatList
    elseif name === :any_list
        return (obj.__protobuf_jl_internal_values[name])::CollectionDef_AnyList
    else
        getfield(obj, name)
    end
end

mutable struct MetaGraphDef_CollectionDefEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function MetaGraphDef_CollectionDefEntry(; kwargs...)
        obj = new(meta(MetaGraphDef_CollectionDefEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct MetaGraphDef_CollectionDefEntry (mapentry)
const __meta_MetaGraphDef_CollectionDefEntry = Ref{ProtoMeta}()
function meta(::Type{MetaGraphDef_CollectionDefEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_MetaGraphDef_CollectionDefEntry)
            __meta_MetaGraphDef_CollectionDefEntry[] = target = ProtoMeta(MetaGraphDef_CollectionDefEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => CollectionDef]
            meta(target, MetaGraphDef_CollectionDefEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_MetaGraphDef_CollectionDefEntry[]
    end
end
function Base.getproperty(obj::MetaGraphDef_CollectionDefEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::CollectionDef
    else
        getfield(obj, name)
    end
end

mutable struct MetaGraphDef_SignatureDefEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function MetaGraphDef_SignatureDefEntry(; kwargs...)
        obj = new(meta(MetaGraphDef_SignatureDefEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct MetaGraphDef_SignatureDefEntry (mapentry) (has cyclic type dependency)
const __meta_MetaGraphDef_SignatureDefEntry = Ref{ProtoMeta}()
function meta(::Type{MetaGraphDef_SignatureDefEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_MetaGraphDef_SignatureDefEntry)
            __meta_MetaGraphDef_SignatureDefEntry[] = target = ProtoMeta(MetaGraphDef_SignatureDefEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => "SignatureDef"]
            meta(target, MetaGraphDef_SignatureDefEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_MetaGraphDef_SignatureDefEntry[]
    end
end
function Base.getproperty(obj::MetaGraphDef_SignatureDefEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

mutable struct MetaGraphDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function MetaGraphDef(; kwargs...)
        obj = new(meta(MetaGraphDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct MetaGraphDef (has cyclic type dependency)
const __meta_MetaGraphDef = Ref{ProtoMeta}()
function meta(::Type{MetaGraphDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_MetaGraphDef)
            __meta_MetaGraphDef[] = target = ProtoMeta(MetaGraphDef)
            allflds = Pair{Symbol,Union{Type,String}}[:meta_info_def => MetaGraphDef_MetaInfoDef, :graph_def => GraphDef, :saver_def => SaverDef, :collection_def => "Base.Dict{AbstractString,CollectionDef}", :signature_def => "Base.Dict{AbstractString,SignatureDef}", :asset_file_def => "Base.Vector{AssetFileDef}", :object_graph_def => SavedObjectGraph]
            meta(target, MetaGraphDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_MetaGraphDef[]
    end
end
function Base.getproperty(obj::MetaGraphDef, name::Symbol)
    if name === :meta_info_def
        return (obj.__protobuf_jl_internal_values[name])::MetaGraphDef_MetaInfoDef
    elseif name === :graph_def
        return (obj.__protobuf_jl_internal_values[name])::GraphDef
    elseif name === :saver_def
        return (obj.__protobuf_jl_internal_values[name])::SaverDef
    elseif name === :collection_def
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :signature_def
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :asset_file_def
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :object_graph_def
        return (obj.__protobuf_jl_internal_values[name])::SavedObjectGraph
    else
        getfield(obj, name)
    end
end

mutable struct TensorInfo_CompositeTensor <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TensorInfo_CompositeTensor(; kwargs...)
        obj = new(meta(TensorInfo_CompositeTensor), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorInfo_CompositeTensor (has cyclic type dependency)
const __meta_TensorInfo_CompositeTensor = Ref{ProtoMeta}()
function meta(::Type{TensorInfo_CompositeTensor})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorInfo_CompositeTensor)
            __meta_TensorInfo_CompositeTensor[] = target = ProtoMeta(TensorInfo_CompositeTensor)
            allflds = Pair{Symbol,Union{Type,String}}[:type_spec => TypeSpecProto, :components => "Base.Vector{TensorInfo}"]
            meta(target, TensorInfo_CompositeTensor, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TensorInfo_CompositeTensor[]
    end
end
function Base.getproperty(obj::TensorInfo_CompositeTensor, name::Symbol)
    if name === :type_spec
        return (obj.__protobuf_jl_internal_values[name])::TypeSpecProto
    elseif name === :components
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

mutable struct TensorInfo <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TensorInfo(; kwargs...)
        obj = new(meta(TensorInfo), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorInfo (has cyclic type dependency)
const __meta_TensorInfo = Ref{ProtoMeta}()
function meta(::Type{TensorInfo})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorInfo)
            __meta_TensorInfo[] = target = ProtoMeta(TensorInfo)
            fnum = Int[1,4,5,2,3]
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :coo_sparse => TensorInfo_CooSparse, :composite_tensor => TensorInfo_CompositeTensor, :dtype => Int32, :tensor_shape => TensorShapeProto]
            oneofs = Int[1,1,1,0,0]
            oneof_names = Symbol[Symbol("encoding")]
            meta(target, TensorInfo, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_TensorInfo[]
    end
end
function Base.getproperty(obj::TensorInfo, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :coo_sparse
        return (obj.__protobuf_jl_internal_values[name])::TensorInfo_CooSparse
    elseif name === :composite_tensor
        return (obj.__protobuf_jl_internal_values[name])::TensorInfo_CompositeTensor
    elseif name === :dtype
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :tensor_shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    else
        getfield(obj, name)
    end
end

mutable struct SignatureDef_InputsEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SignatureDef_InputsEntry(; kwargs...)
        obj = new(meta(SignatureDef_InputsEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SignatureDef_InputsEntry (mapentry) (has cyclic type dependency)
const __meta_SignatureDef_InputsEntry = Ref{ProtoMeta}()
function meta(::Type{SignatureDef_InputsEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SignatureDef_InputsEntry)
            __meta_SignatureDef_InputsEntry[] = target = ProtoMeta(SignatureDef_InputsEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => TensorInfo]
            meta(target, SignatureDef_InputsEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SignatureDef_InputsEntry[]
    end
end
function Base.getproperty(obj::SignatureDef_InputsEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::TensorInfo
    else
        getfield(obj, name)
    end
end

mutable struct SignatureDef_OutputsEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SignatureDef_OutputsEntry(; kwargs...)
        obj = new(meta(SignatureDef_OutputsEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SignatureDef_OutputsEntry (mapentry) (has cyclic type dependency)
const __meta_SignatureDef_OutputsEntry = Ref{ProtoMeta}()
function meta(::Type{SignatureDef_OutputsEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SignatureDef_OutputsEntry)
            __meta_SignatureDef_OutputsEntry[] = target = ProtoMeta(SignatureDef_OutputsEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => TensorInfo]
            meta(target, SignatureDef_OutputsEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SignatureDef_OutputsEntry[]
    end
end
function Base.getproperty(obj::SignatureDef_OutputsEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::TensorInfo
    else
        getfield(obj, name)
    end
end

mutable struct SignatureDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SignatureDef(; kwargs...)
        obj = new(meta(SignatureDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SignatureDef (has cyclic type dependency)
const __meta_SignatureDef = Ref{ProtoMeta}()
function meta(::Type{SignatureDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SignatureDef)
            __meta_SignatureDef[] = target = ProtoMeta(SignatureDef)
            allflds = Pair{Symbol,Union{Type,String}}[:inputs => "Base.Dict{AbstractString,TensorInfo}", :outputs => "Base.Dict{AbstractString,TensorInfo}", :method_name => AbstractString]
            meta(target, SignatureDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SignatureDef[]
    end
end
function Base.getproperty(obj::SignatureDef, name::Symbol)
    if name === :inputs
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :outputs
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :method_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct AssetFileDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AssetFileDef(; kwargs...)
        obj = new(meta(AssetFileDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AssetFileDef (has cyclic type dependency)
const __meta_AssetFileDef = Ref{ProtoMeta}()
function meta(::Type{AssetFileDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AssetFileDef)
            __meta_AssetFileDef[] = target = ProtoMeta(AssetFileDef)
            allflds = Pair{Symbol,Union{Type,String}}[:tensor_info => TensorInfo, :filename => AbstractString]
            meta(target, AssetFileDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AssetFileDef[]
    end
end
function Base.getproperty(obj::AssetFileDef, name::Symbol)
    if name === :tensor_info
        return (obj.__protobuf_jl_internal_values[name])::TensorInfo
    elseif name === :filename
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

export MetaGraphDef_MetaInfoDef_FunctionAliasesEntry, MetaGraphDef_MetaInfoDef, MetaGraphDef_CollectionDefEntry, MetaGraphDef_SignatureDefEntry, MetaGraphDef, CollectionDef_NodeList, CollectionDef_BytesList, CollectionDef_Int64List, CollectionDef_FloatList, CollectionDef_AnyList, CollectionDef, TensorInfo_CooSparse, TensorInfo_CompositeTensor, TensorInfo, SignatureDef_InputsEntry, SignatureDef_OutputsEntry, SignatureDef, AssetFileDef, MetaGraphDef_SignatureDefEntry, MetaGraphDef, TensorInfo_CompositeTensor, TensorInfo, SignatureDef_InputsEntry, SignatureDef_OutputsEntry, SignatureDef, AssetFileDef
# mapentries: "MetaGraphDef_MetaInfoDef_FunctionAliasesEntry" => ("AbstractString", "AbstractString"), "MetaGraphDef_SignatureDefEntry" => ("AbstractString", "SignatureDef"), "MetaGraphDef_CollectionDefEntry" => ("AbstractString", "CollectionDef"), "SignatureDef_InputsEntry" => ("AbstractString", "TensorInfo"), "SignatureDef_OutputsEntry" => ("AbstractString", "TensorInfo")


text:
using InfiniteOpt, JuMP, MathOptInterface, Distributions
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const JuMPC = JuMP.Containers
const IOTO = InfiniteOpt.TranscriptionOpt
using Test
using Random
const MEM = InfiniteOpt.MeasureEvalMethods
using FastGaussQuadrature
const FGQ = FastGaussQuadrature
using MutableArithmetics
const MA = MutableArithmetics

include("utilities.jl")

# Run unit tests
println("-----------------------------------------")
println("----------------UNIT TESTS---------------")
println("-----------------------------------------")
@time @testset "Datatypes" begin include("datatypes.jl") end
println("")
@time @testset "Utilities" begin include("utility_tests.jl") end
println("")
@time @testset "Infinite Set Methods" begin include("infinite_sets.jl") end
println("")
@time @testset "Parameter Methods" begin include("parameters.jl") end
println("")
@time @testset "Optimizer Setup Methods" begin include("optimizer_setup.jl") end
println("")
@time @testset "Variable Methods" begin
    @testset "Basic " begin include("variable_basics.jl") end
    @testset "Info" begin include("variable_info.jl") end
    @testset "Definition" begin include("variable_definition.jl") end
    @testset "Macros" begin include("variable_macros.jl") end
    @testset "Query" begin include("variable_queries.jl") end
    @testset "Reduction Variables" begin include("reduction_variables.jl") end
end
println("")
@time @testset "Operators" begin include("operators.jl") end
println("")
@time @testset "Expression Methods" begin include("expressions.jl") end
println("")
@time @testset "MutableArithmetics Methods" begin include("mutable_arithmetics.jl") end
println("")
@time @testset "Macro Expressions" begin include("macro_expressions.jl") end
println("")
@time @testset "Measure Evaluation Methods" begin
    include("MeasureEvalMethods/methods.jl")
end
println("")
@time @testset "Measure Methods" begin include("measures.jl") end
println("")
@time @testset "Objective Methods" begin include("objective.jl") end
println("")
@time @testset "Constraint Methods" begin include("constraints.jl") end
println("")
@time @testset "Deletion Methods" begin include("deletion.jl") end
println("")
@time @testset "Expansion Methods" begin include("measure_expansions.jl") end
println("")
@time @testset "Printing Methods" begin include("show.jl") end
println("")
@time @testset "TranscriptionOpt" begin
    @testset "Model" begin include("TranscriptionOpt/model.jl") end
    @testset "Transcribe" begin include("TranscriptionOpt/transcribe.jl") end
    @testset "Optimize" begin include("TranscriptionOpt/optimize.jl") end
end
println("")
@time @testset "Solution Methods" begin include("optimizer.jl") end
println("")
@time @testset "Solution Queries" begin include("results.jl") end
println("")
@time @testset "Extensions" begin include("extensions.jl") end
println("")

# TODO add involved deletion tests
println("-----------------------------------------")
println("-------------TESTING COMPLETE------------")
println("-----------------------------------------")