text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashbio_molecule3dviewer

"""
    dashbio_molecule3dviewer(;kwargs...)

A Molecule3dViewer component.
The Molecule3dViewer component is used to render schematic diagrams
of biomolecules. It can display ribbon-structure diagrams, or
render atoms in the molecule as sticks or spheres.
Read more about the component here:
https://github.com/Autodesk/molecule-3d-for-react
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in callbacks
- `atomLabelsShown` (Bool; optional): Property to either show or hide labels
- `backgroundColor` (String; optional): Property to change the background color of the molecule viewer
- `backgroundOpacity` (Real; optional): Property to change the background opacity - ranges from 0 to 1
- `labels` (Array of Dicts; optional): Labels corresponding to the atoms of the molecule.
Each label has a `text` field, a string containing the label content,
and can have many other styling fields as described in
https://3dmol.csb.pitt.edu/doc/types.html#LabelSpec
- `modelData` (optional): The data that will be used to display the molecule in 3D
The data will be in JSON format
and should have two main dictionaries - atoms, bonds. modelData has the following type: lists containing elements 'atoms', 'bonds'.
Those elements have the following types:
  - `atoms` (Array; optional)
  - `bonds` (Array; optional)
- `orbital` (optional): Add an isosurface from volumetric data provided in the `cube_file`. orbital has the following type: lists containing elements 'cube_file', 'iso_val', 'opacity', 'positiveVolumetricColor', 'negativeVolumetricColor'.
Those elements have the following types:
  - `cube_file` (String; optional): The filepath containing raw volumetric data for vertex coloring
  - `iso_val` (Real; optional): The isovalue to draw the surface at
  - `opacity` (Real; optional): Transparency of the surface, between 0 and 1
  - `positiveVolumetricColor` (String; optional): Color for the positive value of the isosurface orbital
  - `negativeVolumetricColor` (String; optional): Color for the negative value of the isosurface orbital
- `selectedAtomIds` (Array; optional): Property that stores a list of all selected atoms
- `selectionType` (a value equal to: 'atom', 'residue', 'chain'; optional): The selection type - may be atom, residue or chain
- `shapes` (Array of Dicts; optional): Add a predefined renderable shape objects to the molecule.
Valid shape types are Arrow, Sphere, and Cylinder.
- `styles` (optional): Property that can be used to change the representation of
the molecule. Options include sticks, cartoon and sphere. styles has the following type: Array of lists containing elements 'color', 'visualization_type'.
Those elements have the following types:
  - `color` (String; optional)
  - `visualization_type` (a value equal to: 'cartoon', 'sphere', 'stick'; optional)s
- `zoom` (optional): Zoom the current view by a constant factor, with optional parameters
to modify the duration and motion of the zoom animation.. zoom has the following type: lists containing elements 'factor', 'animationDuration', 'fixedPath'.
Those elements have the following types:
  - `factor` (Real; optional): Magnification factor. Values greater than 1 will zoom,
in, less than one will zoom out. Default 2.
  - `animationDuration` (Real; optional): An optional parameter that denotes the duration of a
zoom animation, in milliseconds.
  - `fixedPath` (Bool; optional): If true, animation is constrained to requested motion,
overriding updates that happen during the animation.
- `zoomTo` (optional): Zoom to center of atom selection.. zoomTo has the following type: lists containing elements 'sel', 'animationDuration', 'fixedPath'.
Those elements have the following types:
  - `sel` (optional): Selection specification specifying model and atom properties
to select. Default: all atoms in viewer.. sel has the following type: lists containing elements 'chain', 'res'.
Those elements have the following types:
  - `chain` (String; optional): Chain that the residue is located on.
  - `res` (Real; optional): The index value used to identify the residue;
residues are numbered sequentially starting from 1.
  - `animationDuration` (Real; optional): An optional parameter that denotes the duration of a zoom animation
, in milliseconds.
  - `fixedPath` (Bool; optional): If true, animation is constrained to requested motion,
overriding updates that happen during the animation.
"""
function dashbio_molecule3dviewer(; kwargs...)
        available_props = Symbol[:id, :atomLabelsShown, :backgroundColor, :backgroundOpacity, :labels, :modelData, :orbital, :selectedAtomIds, :selectionType, :shapes, :styles, :zoom, :zoomTo]
        wild_props = Symbol[]
        return Component("dashbio_molecule3dviewer", "Molecule3dViewer", "dash_bio", available_props, wild_props; kwargs...)
end


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module DSGE
    using ModelConstructors, SMC
    using Dates, Test, BenchmarkTools
    using Distributed, Distributions, FileIO, FredData, HDF5, JLD2, LinearAlgebra
    using Missings, Nullables, Optim, Printf, Random, RecipesBase, SparseArrays, SpecialFunctions
    using StateSpaceRoutines, StatsPlots
    using CSV, DataFrames, DataStructures, OrderedCollections
    using DataStructures: SortedDict, insert!, ForwardOrdering
    using Roots: fzero, ConvergenceFailed
    using StatsBase: sample, Weights
    using StatsFuns: chisqinvcdf
    using Statistics: std

    import Calculus
    import Base.isempty, Base.<, Base.min, Base.max
    import LinearAlgebra: rank
    import Optim: optimize, SecondOrderOptimizer, MultivariateOptimizationResults
    import StateSpaceRoutines: KalmanFilter, augment_states_with_shocks, solve_discrete_lyapunov
    import ModelConstructors
    import ModelConstructors: posterior!, posterior, <=, n_states,
                              n_states, n_states_augmented, n_shocks_exogenous,
                              n_shocks_expectational, n_observables, n_pseudo_observables,
                              n_equilibrium_conditions, n_parameters, n_parameters_steady_state,
                              n_parameters_free, SteadyStateParameterGrid
    import SMC: get_vals
    import Calculus, Missings, Nullables
    import StateSpaceRoutines: KalmanFilter
    import SparseArrays: sparse, spdiagm, spzeros

    export
        # defaults.jl
        default_settings!, default_test_settings!,

        # abstractdsgemodel.jl
        AbstractDSGEModel, AbstractRepModel, description,
        n_anticipated_shocks, n_anticipated_shocks_padding,
        date_presample_start, date_mainsample_start, date_zlb_start,
        date_presample_end, date_prezlb_end, date_mainsample_end, date_conditional_end,
        index_presample_start, index_mainsample_start, index_zlb_start, index_forecast_start,
        index_shockdec_start,
        n_presample_periods, n_prezlb_periods, n_zlb_periods, n_mainsample_periods,
        n_conditional_periods,
        inds_presample_periods, inds_prezlb_periods, inds_zlb_periods, inds_mainsample_periods,
        n_states, n_states_augmented, n_shocks_exogenous, n_shocks_expectational,
        n_equilibrium_conditions, n_observables, n_parameters, n_parameters_steady_state,
        n_parameters_free, n_pseudo_observables, get_dict, get_key,
        inds_states_no_ant, inds_shocks_no_ant, inds_obs_no_ant,
        spec, subspec, saveroot, dataroot,
        data_vintage, data_id, cond_vintage, cond_id, cond_full_names, cond_semi_names,
        use_population_forecast,
        use_parallel_workers,
        reoptimize, calculate_hessian, hessian_path, n_hessian_test_params,
        n_mh_blocks, n_mh_simulations, n_mh_burn, mh_thin,
        date_forecast_start, date_forecast_end,
        forecast_block_size, forecast_start_block,
        forecast_input_file_overrides, forecast_uncertainty_override,
        forecast_smoother, forecast_horizons,
        forecast_zlb_value, forecast_tdist_shocks, forecast_tdist_df_val,
        shockdec_startdate, date_shockdec_end,
        n_shockdec_periods, impulse_response_horizons,
        load_parameters_from_file, specify_mode!, specify_hessian!,
        logpath, workpath, rawpath, tablespath, figurespath, inpath,
        transform_to_model_space!, transform_to_real_line!,
        ShockGroup, alternative_policy,

        # statespace.jl
        Transition, Measurement, PseudoMeasurement, System, compute_system,

        # benchmark/
        print_all_benchmarks, construct_trial_group, write_ref_trial, write_ref_trial_group,

        # data/
        load_data, load_data_levels, load_cond_data_levels, load_fred_data,
        transform_data, save_data, get_data_filename,
        df_to_matrix, hpfilter, difflog, quartertodate, percapita, nominal_to_real,
        oneqtrpctchange, annualtoquarter, quartertoannual, quartertoannualpercent,
        loggrowthtopct_percapita, loggrowthtopct, logleveltopct_annualized,
        loggrowthtopct_annualized_percapita, loggrowthtopct_annualized,
        logleveltopct_annualized_percapita,
        logleveltopct_annualized_approx, loggrowthtopct_4q_approx, logleveltopct_4q_approx,
        parse_data_series, collect_data_transforms, reverse_transform,
        subtract_quarters, iterate_quarters,
        simulate_data, simulate_observables, simulate_states,

        # solve/
        gensys, solve, klein,

        # estimate/
        simulated_annealing, combined_optimizer, lbfgs,
        filter, filter_shocks, likelihood, posterior, posterior!,
        optimize!, csminwel, hessian!, estimate, proposal_distribution,
        metropolis_hastings, compute_parameter_covariance, prior, get_estimation_output_files,
        compute_moments, find_density_bands, mutation, resample, smc,
        mvnormal_mixture_draw, nearest_spd, marginal_data_density,
        initial_draw!, ParticleCloud, Particle, estimate_bma,

        # backwards_compatibility.jl
        smc2, old_to_new_cloud, # TO REMOVE

        # forecast/
        load_draws, forecast_one,
        smooth, forecast, shock_decompositions, deterministic_trends, trends, impulse_responses,
        compute_system, compute_system_function, add_requisite_output_vars, n_forecast_draws,
        get_forecast_input_file, get_forecast_output_files, get_forecast_filename,
        read_forecast_output,

        # analysis/
        find_density_bands, moment_tables, compute_meansbands, MeansBands,
        meansbands_to_matrix, read_mb, read_bdd_and_unbdd_mb,
        get_meansbands_input_file, get_meansbands_output_file, get_product, get_class,
        which_density_bands,
        prepare_meansbands_tables_timeseries, prepare_means_tables_shockdec,
        prepare_meansbands_table_irf,
        write_meansbands_tables_timeseries, write_means_tables_shockdec,
        prepare_meansbands_table_irf,
        write_meansbands_tables_all, construct_fcast_and_hist_dfs,
        df_to_table, load_posterior_moments, sample_λ, compute_Eλ,
        create_q4q4_mb,

        # decomp/
        decompose_forecast, decomposition_means,

        # decomp/
        decompose_forecast, decomposition_means,

        # altpolicy/
        AltPolicy, taylor93, taylor99,

        # scenarios/
        AbstractScenario, SingleScenario, Scenario, SwitchingScenario, ScenarioAggregate,
        n_targets, n_instruments, n_target_horizons, targets_to_data,
        compute_scenario_system, filter_shocks!, forecast_scenario, simulate_switching,
        scenario_means_bands,
        get_scenario_input_file, n_scenario_draws, get_scenario_filename,
        get_scenario_output_files,
        read_scenario_output, get_scenario_mb_input_file, get_scenario_mb_output_file,
        read_scenario_mb,
        count_scenario_draws,

        # packet/
        usual_settings!, usual_forecast, write_forecast_centric_packet,
        write_standard_packet, plot_standard_packet,

        # plot/
        plot_prior_posterior, plot_impulse_response, plot_history_and_forecast, hair_plot,
        plot_forecast_comparison, plot_shock_decomposition, plot_altpolicies, plot_scenario,
        plot_posterior_intervals, plot_posterior_interval_comparison, plot_forecast_decomposition,
        plot_forecast_sequence, date_to_float,

        # models/
        init_parameters!, steadystate!, init_observable_mappings!,
        init_pseudo_observable_mappings!,
        Model990, Model1002, Model1010, Model805, Model904, SmetsWouters, SmetsWoutersOrig, AnSchorfheide,
        PoolModel, eqcond, measurement, pseudo_measurement,
        shock_groupings, transition,

        # models/heterogeneous/
        KrusellSmith, BondLabor, RealBond, RealBondMkup, HetDSGE, HetDSGEGovDebt,
        RepDSGEGovDebt, HetDSGESimpleTaylor, HetDSGELag, Grid,

        #### Continuous time
        # models
        solve_hjb, solve_kfe, model_settings!, AbstractCTModel, KrusellSmithCT,
        SteadyStateParameterArray, OneAssetHANK, calibrate_pLH_pHL,

	    # TwoAssetHANK
        TwoAssetHANK,

        # solve/
        gensysct, gensysct!, new_divct, decomposition_svdct!, <,
        krylov_reduction, valueref_reduction, deflated_block_arnoldi, change_basis,
        oneDquad_spline, extend_to_nd, projection_for_subset, spline_basis,
        solve_static_conditions,

        # estimate
        hessizero, hess_diag_element, hess_offdiag_element, transform_transition_matrices,

        # estimate/ct_filters
        BlockKalmanFilter, init_stationary_states, block_kalman_filter, CTBlockKalmanFilter,
        ct_block_kalman_filter, ct_kalman_filter, forecast!

    const VERBOSITY = Dict(:none => 0, :low => 1, :high => 2)
    const DSGE_DATE_FORMAT = "yymmdd"
    const DSGE_DATASERIES_DELIM = "__"
    const DSGE_SHOCKDEC_DELIM = "__"

    include("abstractdsgemodel.jl")
    include("defaults.jl")
    include("models/poolmodel/poolmodel.jl")
    include("statespace.jl")
    include("util.jl")
    include("grids.jl")

    include("benchmark/util.jl")
    include("benchmark/benchmark.jl")
    include("benchmark/io.jl")

    include("data/load_data.jl")
    include("data/fred_data.jl")
    include("data/transformations.jl")
    include("data/transform_data.jl")
    include("data/reverse_transform.jl")
    include("data/util.jl")
    include("data/simulate_data.jl")

    include("solve/gensys.jl")
    include("solve/solve.jl")
    include("solve/klein.jl")

    include("estimate/util.jl")
    include("estimate/kalman.jl")
    include("estimate/filter.jl")
    include("estimate/posterior.jl")
    include("estimate/optimize.jl")
    include("estimate/csminwel.jl")
    include("estimate/hessian.jl")
    include("estimate/hessizero.jl")
    include("estimate/simulated_annealing.jl")
    include("estimate/combined_optimizer.jl")
    include("estimate/lbfgs.jl")
    include("estimate/nelder_mead.jl")
    include("estimate/marginal_data_density.jl")
    include("estimate/estimate.jl")
    include("estimate/estimate_bma.jl")
    include("estimate/nearest_spd.jl")

    include("estimate/metropolis_hastings.jl")

    include("estimate/smc/particle.jl")
    #include("estimate/smc/initialization.jl")
    #include("estimate/smc/helpers.jl")
    #include("estimate/smc/util.jl")
    #include("estimate/smc/mutation.jl")
    #include("estimate/smc/resample.jl")
    #include("estimate/smc/smc.jl")

    include("estimate/smc.jl")
    include("estimate/backwards_compatibility.jl")

    # CT HANK code
    include("estimate/filter_hank.jl")
    include("estimate/hessizero_hank.jl")
    include("estimate/transform_transition_matrices.jl")
    include("estimate/ct_filters/ct_kalman_filter.jl")
    include("estimate/ct_filters/block_kalman_filter.jl")
    # include("estimate/ct_filters/ct_block_kalman_filter.jl")

    include("analysis/moments.jl")
    include("analysis/meansbands.jl")
    include("analysis/compute_meansbands.jl")
    include("analysis/meansbands_to_matrix.jl")
    include("analysis/io.jl")
    include("analysis/util.jl")
    include("analysis/df_to_table.jl")
    include("analysis/create_q4q4_mb.jl")

    include("decomp/drivers.jl")
    include("decomp/io.jl")
    include("decomp/meansbands.jl")

    include("altpolicy/altpolicy.jl")
    include("altpolicy/taylor93.jl")
    include("altpolicy/taylor99.jl")

    include("scenarios/scenario.jl")
    include("scenarios/io.jl")
    include("scenarios/forecast.jl")
    include("scenarios/switching.jl")
    include("scenarios/transform.jl")

    include("packet/util.jl")
    include("packet/drivers.jl")
    include("packet/standard_packet.jl")
    include("packet/forecast_centric_packet.jl")

    include("plot/util.jl")
    include("plot/plot_posterior_intervals.jl")
    include("plot/plot_prior_posterior.jl")
    include("plot/plot_impulse_response.jl")
    include("plot/plot_history_and_forecast.jl")
    include("plot/hair_plot.jl")
    include("plot/plot_forecast_comparison.jl")
    include("plot/plot_shock_decomposition.jl")
    include("plot/plot_altpolicies.jl")
    include("plot/plot_scenario.jl")
    include("plot/plot_forecast_decomposition.jl")
    include("plot/plot_forecast_sequence.jl")

    # models/
    include("models/financial_frictions.jl")

    # Representative Agent Models
    include("models/representative/m990/m990.jl")
    include("models/representative/m990/subspecs.jl")
    include("models/representative/m990/eqcond.jl")
    include("models/representative/m990/observables.jl")
    include("models/representative/m990/measurement.jl")
    include("models/representative/m990/pseudo_observables.jl")
    include("models/representative/m990/pseudo_measurement.jl")
    include("models/representative/m990/augment_states.jl")

    include("models/representative/m1002/m1002.jl")
    include("models/representative/m1002/subspecs.jl")
    include("models/representative/m1002/eqcond.jl")
    include("models/representative/m1002/observables.jl")
    include("models/representative/m1002/measurement.jl")
    include("models/representative/m1002/pseudo_observables.jl")
    include("models/representative/m1002/pseudo_measurement.jl")
    include("models/representative/m1002/augment_states.jl")

    include("models/representative/m1010/m1010.jl")
    include("models/representative/m1010/subspecs.jl")
    include("models/representative/m1010/eqcond.jl")
    include("models/representative/m1010/observables.jl")
    include("models/representative/m1010/measurement.jl")
    include("models/representative/m1010/pseudo_observables.jl")
    include("models/representative/m1010/pseudo_measurement.jl")
    include("models/representative/m1010/augment_states.jl")

    include("models/representative/m805/m805.jl")
    include("models/representative/m805/subspecs.jl")
    include("models/representative/m805/eqcond.jl")
    include("models/representative/m805/observables.jl")
    include("models/representative/m805/measurement.jl")
    include("models/representative/m805/augment_states.jl")

    include("models/representative/m904/m904.jl")
    include("models/representative/m904/subspecs.jl")
    include("models/representative/m904/eqcond.jl")
    include("models/representative/m904/observables.jl")
    include("models/representative/m904/measurement.jl")
    include("models/representative/m904/pseudo_observables.jl")
    include("models/representative/m904/pseudo_measurement.jl")
    include("models/representative/m904/augment_states.jl")

    include("models/representative/smets_wouters/smets_wouters.jl")
    include("models/representative/smets_wouters/subspecs.jl")
    include("models/representative/smets_wouters/eqcond.jl")
    include("models/representative/smets_wouters/observables.jl")
    include("models/representative/smets_wouters/measurement.jl")
    include("models/representative/smets_wouters/augment_states.jl")

    include("models/representative/smets_wouters_orig/smets_wouters_orig.jl")
    include("models/representative/smets_wouters_orig/subspecs.jl")
    include("models/representative/smets_wouters_orig/eqcond.jl")
    include("models/representative/smets_wouters_orig/observables.jl")
    include("models/representative/smets_wouters_orig/measurement.jl")
    include("models/representative/smets_wouters_orig/augment_states.jl")

    include("models/representative/an_schorfheide/an_schorfheide.jl")
    include("models/representative/an_schorfheide/subspecs.jl")
    include("models/representative/an_schorfheide/eqcond.jl")
    include("models/representative/an_schorfheide/observables.jl")
    include("models/representative/an_schorfheide/measurement.jl")
    include("models/representative/an_schorfheide/pseudo_observables.jl")
    include("models/representative/an_schorfheide/pseudo_measurement.jl")
    include("models/representative/an_schorfheide/augment_states.jl")

    # PoolModel
    include("models/poolmodel/subspecs.jl")
    include("models/poolmodel/observables.jl")
    include("models/poolmodel/measurement.jl")
    include("models/poolmodel/transition.jl")

    # Heterogeneous Agent Models
    include("models/heterogeneous/util.jl")

    include("models/heterogeneous/krusell_smith/krusell_smith.jl")
    include("models/heterogeneous/krusell_smith/steady_state.jl")
    include("models/heterogeneous/krusell_smith/subspecs.jl")
    include("models/heterogeneous/krusell_smith/jacobian.jl")
    include("models/heterogeneous/krusell_smith/shock_loading.jl")
    include("models/heterogeneous/krusell_smith/observables.jl")
    include("models/heterogeneous/krusell_smith/measurement.jl")

    include("models/heterogeneous/bond_labor/bond_labor.jl")
    include("models/heterogeneous/bond_labor/steady_state.jl")
    include("models/heterogeneous/bond_labor/jacobian.jl")
    include("models/heterogeneous/bond_labor/shock_loading.jl")
    include("models/heterogeneous/bond_labor/observables.jl")
    include("models/heterogeneous/bond_labor/measurement.jl")

    include("models/heterogeneous/real_bond/real_bond.jl")
    include("models/heterogeneous/real_bond/steady_state.jl")
    include("models/heterogeneous/real_bond/jacobian.jl")
    include("models/heterogeneous/real_bond/shock_loading.jl")
    include("models/heterogeneous/real_bond/observables.jl")
    include("models/heterogeneous/real_bond/measurement.jl")

    include("models/heterogeneous/real_bond_mkup/real_bond_mkup.jl")
    include("models/heterogeneous/real_bond_mkup/steady_state.jl")
    include("models/heterogeneous/real_bond_mkup/subspecs.jl")
    include("models/heterogeneous/real_bond_mkup/jacobian.jl")
    include("models/heterogeneous/real_bond_mkup/shock_loading.jl")
    include("models/heterogeneous/real_bond_mkup/observables.jl")
    include("models/heterogeneous/real_bond_mkup/measurement.jl")
    include("models/heterogeneous/real_bond_mkup/augment_states.jl")

    include("models/heterogeneous/het_dsge/het_dsge.jl")
    include("models/heterogeneous/het_dsge/steady_state.jl")
    include("models/heterogeneous/het_dsge/subspecs.jl")
    include("models/heterogeneous/het_dsge/jacobian.jl")
    include("models/heterogeneous/het_dsge/shock_loading.jl")
    include("models/heterogeneous/het_dsge/observables.jl")
    include("models/heterogeneous/het_dsge/measurement.jl")
    include("models/heterogeneous/het_dsge/augment_states.jl")

    include("models/heterogeneous/het_dsge_gov_debt/util.jl")
    include("models/heterogeneous/het_dsge_gov_debt/het_dsge_gov_debt.jl")
    include("models/heterogeneous/het_dsge_gov_debt/steady_state.jl")
    include("models/heterogeneous/het_dsge_gov_debt/subspecs.jl")
    include("models/heterogeneous/het_dsge_gov_debt/jacobian.jl")
    include("models/heterogeneous/het_dsge_gov_debt/shock_loading.jl")
    include("models/heterogeneous/het_dsge_gov_debt/observables.jl")
    include("models/heterogeneous/het_dsge_gov_debt/measurement.jl")
    include("models/heterogeneous/het_dsge_gov_debt/augment_states.jl")

    include("models/representative/rep_dsge_gov_debt/rep_dsge_gov_debt.jl")
    include("models/representative/rep_dsge_gov_debt/subspecs.jl")
    include("models/representative/rep_dsge_gov_debt/jacobian.jl")
    include("models/representative/rep_dsge_gov_debt/shock_loading.jl")
    include("models/representative/rep_dsge_gov_debt/measurement.jl")
    include("models/representative/rep_dsge_gov_debt/augment_states.jl")

    include("models/heterogeneous/het_dsge_lag/het_dsge_lag.jl")
    include("models/heterogeneous/het_dsge_lag/steady_state.jl")
    include("models/heterogeneous/het_dsge_lag/subspecs.jl")
    include("models/heterogeneous/het_dsge_lag/jacobian.jl")
    include("models/heterogeneous/het_dsge_lag/shock_loading.jl")
    include("models/heterogeneous/het_dsge_lag/observables.jl")
    include("models/heterogeneous/het_dsge_lag/measurement.jl")

    #=
    include("models/heterogeneous/het_dsge_simple_taylor/het_dsge_simple_taylor.jl")
    include("models/heterogeneous/het_dsge_simple_taylor/steady_state.jl")
    include("models/heterogeneous/het_dsge_simple_taylor/subspecs.jl")
    include("models/heterogeneous/het_dsge_simple_taylor/jacobian.jl")
    =#

    # Continuous Time Heterogenous Agent Models
    include("solve/solve_hank.jl")
    include("solve/gensysct.jl")
    include("solve/reduction.jl")
    include("solve/sparse_reduction.jl")

    include("models/heterogeneous/solve_hjb.jl")
    include("models/heterogeneous/solve_kfe.jl")

    include("models/heterogeneous/krusell_smith_ct/krusell_smith_ct.jl")
    include("models/heterogeneous/krusell_smith_ct/measurement.jl")
    include("models/heterogeneous/krusell_smith_ct/eqcond.jl")

    include("models/heterogeneous/one_asset_hank/one_asset_hank.jl")
    include("models/heterogeneous/one_asset_hank/measurement.jl")
    include("models/heterogeneous/one_asset_hank/eqcond.jl")
    include("models/heterogeneous/one_asset_hank/helpers.jl")

    include("models/heterogeneous/two_asset_hank/income_transition.jl")
    include("models/heterogeneous/two_asset_hank/two_asset_hank.jl")
    #include("models/heterogeneous/two_asset_hank/measurement.jl")
    include("models/heterogeneous/two_asset_hank/eqcond.jl")
    include("models/heterogeneous/two_asset_hank/util.jl")
    include("models/heterogeneous/two_asset_hank/helpers.jl")
    include("models/heterogeneous/two_asset_hank/interp.jl")

    include("forecast/util.jl")
    include("forecast/io.jl")
    include("forecast/smooth.jl")
    include("forecast/forecast.jl")
    include("forecast/shock_decompositions.jl")
    include("forecast/impulse_responses.jl")
    include("forecast/drivers.jl")
end

ext: jl


text:
######################################## Model 2 ~ PMCMC/SMC2
# Parameter
μ = [-2., 2.]
σ = [1., 1.]
p = [.05, .95]
# Latent data
latent = rand(_rng, Categorical(p), N)
data = [rand(_rng, Normal(μ[iter], σ[iter])) for iter in latent]

# Create ModelWrapper struct, assuming we do not know latent
latent_init = rand(_rng, Categorical(p), N_SMC2)
data_init = data[1:N_SMC2]

myparameter = (;
    μ = Param(μ, [Normal(-2., 5), Normal(2., 5)]),
    σ = Param(σ, [Gamma(2.,2.), Gamma(2.,2.)]),
    p = Param(p, Dirichlet(2, 2)),
    latent = Param(latent_init, [Categorical(p) for _ in Base.OneTo(N_SMC2)]),
)
mymodel = ModelWrapper(myparameter)
myobjective = Objective(mymodel, data_init)

# Assign an objective for both a particle filter and an mcmc kernel:
myobjective_pf = Objective(mymodel, data_init, :latent)
myobjective_mcmc = Objective(mymodel, data_init, (:μ, :σ, :p))

@test length(myobjective_mcmc.model.val.latent) == N_SMC2

# Assign Model dynamics
function BaytesFilters.dynamics(objective::Objective{<:ModelWrapper{BaseModel}})
    @unpack model, data = objective
    @unpack μ, σ, p = model.val

    initial_latent = Categorical(p)
    transition_latent(particles, iter) = initial_latent
    transition_data(particles, iter) = Normal(μ[particles[iter]], σ[particles[iter]])

    return Markov(initial_latent, transition_latent, transition_data)
end
dynamics(myobjective)

# Assign log target
function (objective::Objective{<:ModelWrapper{BaseModel}})(θ::NamedTuple)
    @unpack model, data, tagged = objective
    @unpack μ, σ, p, latent = θ
## Prior -> a faster shortcut without initializing the priors again
    lprior = log_prior(tagged.info.constraint, ModelWrappers.subset(θ, tagged.parameter) )
##Likelihood
    dynamicsᵉ = [Normal(μ[iter], σ[iter]) for iter in eachindex(μ)]
    dynamicsˢ = Categorical(p)
    ll = 0.0
#FOR PMCMC ~ target p(θ ∣ latent_1:t, data_1:t)
    for iter in eachindex(data)
        ll += logpdf(dynamicsᵉ[latent[iter]], data[iter])
        ll += logpdf(dynamicsˢ, latent[iter] )
    end
#=
# FOR MCMC ~ target p(θ ∣ data_1:t) by integrating out latent_1:t
    for time in eachindex(data)
        ll += logsumexp(logpdf(dynamicsˢ, iter) + logpdf(dynamicsᵉ[iter], grab(data, time)) for iter in eachindex(dynamicsᵉ))
    end
=#
    return ll + lprior
end


text:
using JLD
using SparseArrays

const n_user = 943
const n_item = 1682

R = spzeros(n_user, n_item)

open("ml-100k/u.data", "r") do f
    for line in eachline(f)
        l = split(line, "\t")
        user, item, value = parse(Int, l[1]), parse(Int, l[2]), parse(Int, l[3])
        R[user, item] = value
    end
end

save("ml-100k.jld", "R", R)


text:
SUITE["Reconstructions"] = BenchmarkGroup()
SR = SUITE["Reconstructions"]

srand(1234)
N = 10000; D = 3; B = 3

data = Dataset(rand(N, B))
x = data[:, 1]

combinations = ["standard", "multitime", "multidim", "allmulti"]
taus = (3, [2,7,8], 3, hcat([2,7,8],[2, 5,9],[4,8,12]))

for i in 1:4
    SR[combinations[i]] = BenchmarkGroup()
    if i < 3
        SR[combinations[i]]["full"] = @benchmarkable reconstruct($x, $D, $(taus[i]))
        de = DelayEmbedding(D, taus[i])
        SR[combinations[i]]["vector"] = @benchmarkable $(de)($x, 1)
    else
        SR[combinations[i]]["full"] = @benchmarkable reconstruct($data, $D, $(taus[i]))
        de = MTDelayEmbedding(D, taus[i], B)
        SR[combinations[i]]["vector"] = @benchmarkable $(de)($data, 1)
    end
end


text:
# reference in tree version of CalibrateEmulateSample
prepend!(LOAD_PATH, [joinpath(@__DIR__, "..")])

using EnsembleKalmanProcesses,
    Documenter,
    Plots,  # so that Literate.jl does not capture precompilation output
    Literate

# Gotta set this environment variable when using the GR run-time on CI machines.
# This happens as examples will use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR = joinpath(@__DIR__, "src/literated")

examples_for_literation = ["LossMinimization/loss_minimization.jl", "AerosolActivation/aerosol_activation.jl"]

for example in examples_for_literation
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor())
end

#----------

api = [
    "EnsembleKalmanProcessModule" => "API/EnsembleKalmanProcessModule.md",
    "ParameterDistribution" => "API/ParameterDistribution.md",
    "Observations" => "API/Observations.md",
    "DataStorage" => "API/DataStorage.md",
]

examples = [
    "Template" => "examples/template_example.md",
    "Cloudy" => "examples/Cloudy_example.md",
    "Lorenz" => "examples/lorenz_example.md",
    "Minimization Loss" => "literated/loss_minimization.md",
    "Aerosol activation" => "literated/aerosol_activation.md",
    "HPC interfacing example: ClimateMachine" => "examples/ClimateMachine_example.md",
]

pages = [
    "Home" => "index.md",
    "Installation instructions" => "installation_instructions.md",
    "Prior distributions" => "parameter_distributions.md",
    "Observations" => "observations.md",
    "Ensemble Kalman Inversion" => "ensemble_kalman_inversion.md",
    "Ensemble Kalman Sampler" => "ensemble_kalman_sampler.md",
    "Unscented Kalman Inversion" => "unscented_kalman_inversion.md",
    "Examples" => examples,
    "API" => api,
    "Glossary" => "glossary.md",
]

#----------

format = Documenter.HTML(collapselevel = 1, prettyurls = !isempty(get(ENV, "CI", "")))

makedocs(
    sitename = "EnsembleKalmanProcesses.jl",
    authors = "CliMA Contributors",
    format = format,
    pages = pages,
    modules = [EnsembleKalmanProcesses],
    doctest = true,
    strict = true,
    clean = true,
    checkdocs = :none,
)

if !isempty(get(ENV, "CI", ""))
    deploydocs(
        repo = "github.com/CliMA/EnsembleKalmanProcesses.jl.git",
        versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
        push_preview = true,
        devbranch = "main",
    )
end


text:
"""
    solver_name(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.SolverName`](@ref) attribute is implemented for `model`.
"""
function solver_name(model::MOI.ModelLike, config::TestConfig)
    if config.solve
        @test MOI.get(model, MOI.SolverName()) isa AbstractString
    end
end
unittests["solver_name"] = solver_name

"""
    silent(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.Silent`](@ref) attribute is implemented for `model`.
"""
function silent(model::MOI.ModelLike, config::TestConfig)
    if config.solve
        @test MOI.supports(model, MOI.Silent())
        # Get the current value to restore it at the end of the test
        value = MOI.get(model, MOI.Silent())
        MOI.set(model, MOI.Silent(), !value)
        @test !value == MOI.get(model, MOI.Silent())
        # Check that `set` does not just take `!` of the current value
        MOI.set(model, MOI.Silent(), !value)
        @test !value == MOI.get(model, MOI.Silent())
        MOI.set(model, MOI.Silent(), value)
        @test value == MOI.get(model, MOI.Silent())
    end
end
unittests["silent"] = silent

"""
    time_limit_sec(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.TimeLimitSec`](@ref) attribute is implemented for `model`.
"""
function time_limit_sec(model::MOI.ModelLike, config::TestConfig)
    if config.solve
        @test MOI.supports(model, MOI.TimeLimitSec())
        # Get the current value to restore it at the end of the test
        value = MOI.get(model, MOI.TimeLimitSec())
        MOI.set(model, MOI.TimeLimitSec(), 0.0)
        @test MOI.get(model, MOI.TimeLimitSec()) == 0.0
        MOI.set(model, MOI.TimeLimitSec(), 1.0)
        @test MOI.get(model, MOI.TimeLimitSec()) == 1.0
        MOI.set(model, MOI.TimeLimitSec(), value)
        @test value == MOI.get(model, MOI.TimeLimitSec()) # Equality should hold
    end
end
unittests["time_limit_sec"] = time_limit_sec

"""
    number_threads(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.NumberOfThreads`](@ref) attribute is implemented for `model`.
"""
function number_threads(model::MOI.ModelLike, config::TestConfig)
    if config.solve
        @test MOI.supports(model, MOI.NumberOfThreads())
        # Get the current value to restore it at the end of the test
        value = MOI.get(model, MOI.NumberOfThreads())
        MOI.set(model, MOI.NumberOfThreads(), 1)
        @test MOI.get(model, MOI.NumberOfThreads()) == 1
        MOI.set(model, MOI.NumberOfThreads(), 3)
        @test MOI.get(model, MOI.NumberOfThreads()) == 3
        MOI.set(model, MOI.NumberOfThreads(), value)
        @test value == MOI.get(model, MOI.NumberOfThreads())
    end
end
unittests["number_threads"] = number_threads

"""
    raw_status_string(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.RawStatusString`](@ref) attribute is implemented for
`model`.
"""
function raw_status_string(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(),
            MOI.SingleVariable(x))
    test_model_solution(model, config, objective_value = 0.0,
                        variable_primal = [(x, 0.0)])
    if config.solve
        @test MOI.get(model, MOI.RawStatusString()) isa AbstractString
    end
end
unittests["raw_status_string"] = raw_status_string

"""
    solve_time(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.SolveTime`](@ref) attribute is implemented for `model`.
"""
function solve_time(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(),
            MOI.SingleVariable(x))
    test_model_solution(model, config, objective_value = 0.0,
                        variable_primal = [(x, 0.0)])
    if config.solve
        time = MOI.get(model, MOI.SolveTime())
        @test time ≥ 0.0
    end
end
unittests["solve_time"] = solve_time


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""DEPRECATED - This group version of ControllerRevision is deprecated by apps/v1beta2/ControllerRevision. See the release notes for more information. ControllerRevision implements an immutable snapshot of state data. Clients are responsible for serializing and deserializing the objects that contain their internal state. Once a ControllerRevision has been successfully created, it can not be updated. The API Server will fail validation of all requests that attempt to mutate the Data field. ControllerRevisions may, however, be deleted. Note that, due to its use by both the DaemonSet and StatefulSet controllers for update and rollback, this object is beta. However, it may be subject to name and representation changes in future releases, and clients should not depend on its stability. It is primarily for internal use by controllers.

    IoK8sApiAppsV1beta1ControllerRevision(;
        apiVersion=nothing,
        data=nothing,
        kind=nothing,
        metadata=nothing,
        revision=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - data::IoK8sApimachineryPkgRuntimeRawExtension : Data is the serialized representation of the state.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#metadata
    - revision::Int64 : Revision indicates the revision of the state represented by Data.
"""
mutable struct IoK8sApiAppsV1beta1ControllerRevision <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    data::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgRuntimeRawExtension } # spec name: data
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    revision::Any # spec type: Union{ Nothing, Int64 } # spec name: revision

    function IoK8sApiAppsV1beta1ControllerRevision(;apiVersion=nothing, data=nothing, kind=nothing, metadata=nothing, revision=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta1ControllerRevision, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAppsV1beta1ControllerRevision, Symbol("data"), data)
        setfield!(o, Symbol("data"), data)
        validate_property(IoK8sApiAppsV1beta1ControllerRevision, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAppsV1beta1ControllerRevision, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiAppsV1beta1ControllerRevision, Symbol("revision"), revision)
        setfield!(o, Symbol("revision"), revision)
        o
    end
end # type IoK8sApiAppsV1beta1ControllerRevision

const _property_map_IoK8sApiAppsV1beta1ControllerRevision = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("data")=>Symbol("data"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("revision")=>Symbol("revision"))
const _property_types_IoK8sApiAppsV1beta1ControllerRevision = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("data")=>"IoK8sApimachineryPkgRuntimeRawExtension", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("revision")=>"Int64")
Base.propertynames(::Type{ IoK8sApiAppsV1beta1ControllerRevision }) = collect(keys(_property_map_IoK8sApiAppsV1beta1ControllerRevision))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta1ControllerRevision }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta1ControllerRevision[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta1ControllerRevision }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta1ControllerRevision[property_name]

function check_required(o::IoK8sApiAppsV1beta1ControllerRevision)
    (getproperty(o, Symbol("revision")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta1ControllerRevision }, name::Symbol, val)
end


text:
mutable struct AccumulativeVariables

    TFLUX_CONV :: AbstractArray
    SFLUX_CONV :: AbstractArray

    TFLUX_DEN_x :: AbstractArray
    TFLUX_DEN_y :: AbstractArray
    TFLUX_DEN_z :: AbstractArray

    SFLUX_DEN_x :: AbstractArray
    SFLUX_DEN_y :: AbstractArray
    SFLUX_DEN_z :: AbstractArray

    ∇∇T      :: AbstractArray
    ∇∇S      :: AbstractArray

    dTdt_ent :: AbstractArray
    dSdt_ent :: AbstractArray

    TFLUX_bot       :: AbstractArray
    SFLUX_bot       :: AbstractArray
    SFLUX_top       :: AbstractArray

    function AccumulativeVariables(Nx, Ny, Nz)
        return new(
            zeros(Nz, Nx, Ny),
            zeros(Nz, Nx, Ny),
            zeros(Nz, Nx, Ny),
            zeros(Nz, Nx, Ny+1),
            zeros(Nz+1, Nx, Ny),
            zeros(Nz, Nx, Ny),
            zeros(Nz, Nx, Ny+1),
            zeros(Nz+1, Nx, Ny),
            zeros(Nz, Nx, Ny),
            zeros(Nz, Nx, Ny),
            zeros(Nx, Ny),
            zeros(Nx, Ny),
            zeros(Nx, Ny),
            zeros(Nx, Ny),
            zeros(Nx, Ny),
        )
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


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

mutable struct TrialResult
    tour::Tour
    insertion::Power
    removal::Power
    noise::Power
    score::Float64
end

"""
Select a removal and an insertion method using powers, and then perform
removal followed by insertion on tour.  Operation done in place.
"""
function remove_insert(
    RNGs::Array{MersenneTwister, 1},
    current::Tour,
    best::Tour,
    dist::Array{Int64, 2},
    member::Array{Int64, 1},
    setdist::Distsv,
    sets::Array{Any, 1},
    order_constraints::Array{Constraints, 1},
    powers,
    param::Dict{Symbol, Any},
    phase::Symbol,
    start_set::Int64,
)
    # make a new tour to perform the insertion and deletion on
    trial = Tour(copy(current.tour), current.cost)
    pivot_tour!(RNGs, trial.tour)
    num_removals = rand(RNGs[Threads.threadid()], param[:min_removals]:param[:max_removals])

    removal = power_select(RNGs, powers["removals"], powers["removal_total"], phase)
    if removal.name == "distance"
        sets_to_insert = distance_removal!(
            RNGs,
            trial.tour,
            dist,
            num_removals,
            member,
            removal.value,
            start_set,
        )
    elseif removal.name == "worst"
        sets_to_insert = worst_removal!(
            RNGs,
            trial.tour,
            dist,
            num_removals,
            member,
            removal.value,
            start_set,
        )
    else
        sets_to_insert = segment_removal!(RNGs, trial.tour, num_removals, member, start_set)
    end

    randomize_sets!(RNGs, sets, sets_to_insert)
    # then perform insertion
    insertion = power_select(RNGs, powers["insertions"], powers["insertion_total"], phase)
    noise = power_select(RNGs, powers["noise"], powers["noise_total"], phase)

    if insertion.name == "cheapest"
        cheapest_insertion!(
            trial.tour,
            sets_to_insert,
            dist,
            setdist,
            sets,
            order_constraints,
            member,
            start_set,
        )
    else
        randpdf_insertion!(
            RNGs,
            trial.tour,
            sets_to_insert,
            dist,
            setdist,
            sets,
            order_constraints,
            member,
            insertion.value,
            noise,
            start_set,
        )
    end
    rand(RNGs[Threads.threadid()]) < param[:prob_reopt] && opt_cycle!(
        RNGs,
        trial,
        dist,
        sets,
        order_constraints,
        member,
        param,
        setdist,
        "partial",
        start_set,
    )

    score = 100 * max(current.cost - trial.cost, 0) / current.cost

    return trial, insertion, removal, noise, score
end

"""
Select an integer between 1 and num according to
and exponential distribution with lambda = power
# goes from left of array if power is positive
# and right of array if it is negative
"""
function select_k(RNGs::Array{MersenneTwister, 1}, num::Int64, power::Float64)
    base = (1 / 2)^abs(power)
    # (1 - base^num)/(1 - base) is sum of geometric series
    rand_select = (1 - base^num) / (1 - base) * rand(RNGs[Threads.threadid()])
    bin = 1
    for k in 1:num
        if rand_select < bin
            return (power >= 0 ? (num - k + 1) : k)
        end
        rand_select -= bin
        bin *= base
    end
    return (power >= 0 ? num : 1)
end

"""
selecting a random k in 1 to length(weights) according to power
and then selecting the kth smallest element in weights
"""
function pdf_select(
    RNGs::Array{MersenneTwister, 1},
    weights::Array{Int64, 1},
    power::Float64,
)
    power == 0.0 && return rand(RNGs[Threads.threadid()], 1:length(weights))
    power > 9.0 && return rand_select(RNGs, weights, maximum(weights))
    power < -9.0 && return rand_select(RNGs, weights, minimum(weights))

    # select kth smallest.  If 1 or length(weights), simply return
    k = select_k(RNGs, length(weights), power)
    k == 1 && return rand_select(RNGs, weights, minimum(weights))
    k == length(weights) && return rand_select(RNGs, weights, maximum(weights))
    val = partialsort(weights, k)

    return rand_select(RNGs, weights, val)
end

"""
choose set with pdf_select, and then insert in best place with noise
"""
function randpdf_insertion!(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    sets_to_insert::Array{Int64, 1},
    dist::Array{Int64, 2},
    setdist::Distsv,
    sets::Array{Any, 1},
    order_constraints::Array{Constraints, 1},
    member::Array{Int64, 1},
    power::Float64,
    noise::Power,
    start_set::Int64,
)
    mindist = [typemax(Int64) for i in 1:length(sets_to_insert)]
    for i in 1:length(sets_to_insert)
        set = sets_to_insert[i]
        for vertex in tour
            if setdist.min_sv[set, vertex] < mindist[i]
                mindist[i] = setdist.min_sv[set, vertex]
            end
        end
    end
    new_vertex_in_tour = 0

    while length(sets_to_insert) > 0
        if new_vertex_in_tour != 0
            for i in 1:length(sets_to_insert)
                set = sets_to_insert[i]
                if setdist.min_sv[set, new_vertex_in_tour] < mindist[i]
                    mindist[i] = setdist.min_sv[set, new_vertex_in_tour]
                end
            end
        end
        set_index = pdf_select(RNGs, mindist, power) # select set to insert from pdf
        # find the closest vertex and the best insertion in that vertex
        nearest_set = sets_to_insert[set_index]

        min_insert_idx, max_insert_idx =
            calc_bounds(tour, nearest_set, order_constraints, member)

        if nearest_set != start_set
            min_insert_idx = max(min_insert_idx, 2)
            if max_insert_idx < min_insert_idx
                max_insert_idx = min_insert_idx
            end
        end

        if noise.name == "subset"
            bestv, bestpos = insert_subset_lb(
                RNGs,
                tour,
                dist,
                sets[nearest_set],
                nearest_set,
                setdist,
                noise.value,
            )
        else
            bestv, bestpos = insert_lb(
                RNGs,
                tour,
                dist,
                sets[nearest_set],
                nearest_set,
                min_insert_idx,
                max_insert_idx,
                setdist,
                noise.value,
            )
        end
        insert!(tour, bestpos, bestv)  # perform the insertion
        new_vertex_in_tour = bestv
        # remove the inserted set from data structures
        splice!(sets_to_insert, set_index)
        splice!(mindist, set_index)
    end
end

"""
choose vertex that can be inserted most cheaply, and insert it in that position
"""
function cheapest_insertion!(
    tour::Array{Int64, 1},
    sets_to_insert::Array{Int64, 1},
    dist::Array{Int64, 2},
    setdist::Distsv,
    sets::Array{Any, 1},
    order_constraints::Array{Constraints, 1},
    member::Array{Int64, 1},
    start_set::Int64,
)
    while length(sets_to_insert) > 0
        best_cost = typemax(Int64)
        curr_v = 0
        curr_pos = 0
        best_set = 0
        best_v = 0
        best_pos = 0
        for i in 1:length(sets_to_insert)
            set_ind = sets_to_insert[i]

            min_insert_idx, max_insert_idx =
                calc_bounds(tour, set_ind, order_constraints, member)

            if set_ind != start_set
                min_insert_idx = max(min_insert_idx, 2)
                if max_insert_idx < min_insert_idx
                    max_insert_idx = min_insert_idx
                end
            end

            # find the best place to insert the vertex
            curr_v, curr_pos, cost = insert_cost_lb(
                tour,
                dist,
                sets[set_ind],
                set_ind,
                setdist,
                min_insert_idx,
                max_insert_idx,
                curr_v,
                curr_pos,
                best_cost,
            )
            if cost < best_cost
                best_set = i
                best_cost = cost
                best_v = curr_v
                best_pos = curr_pos
            end
        end

        # now, perform the insertion
        insert!(tour, best_pos, best_v)
        # remove the inserted set from data structures
        splice!(sets_to_insert, best_set)
    end
end

"""
Given a tour and a set, this function finds the vertex in the set with minimum
insertion cost, along with the position of this insertion in the tour.  If
best_position is i, then vertex should be inserted between tour[i-1] and tour[i].
"""
@inline function insert_lb(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    dist::Array{Int64, 2},
    set::Array{Int64, 1},
    setind::Int,
    min_insert_idx::Int,
    max_insert_idx::Int,
    setdist::Distsv,
    noise::Float64,
)
    best_cost = typemax(Int64)
    bestv = 0
    bestpos = 0

    if min_insert_idx == max_insert_idx && max_insert_idx == length(tour) + 1
        v1 = tour[length(tour)]
        for v in set
            insert_cost = dist[v1, v]
            noise > 0.0 && (
                insert_cost +=
                    round(Int64, noise * rand(RNGs[Threads.threadid()]) * abs(insert_cost))
            )
            if insert_cost < best_cost
                best_cost = insert_cost
                bestv = v
            end
        end
        bestpos = max_insert_idx
        return bestv, bestpos
    end

    for i in min_insert_idx:max_insert_idx
        v1 = prev_tour(tour, i)
        lb =
            setdist.vert_set[v1, setind] + setdist.set_vert[setind, tour[i]] -
            dist[v1, tour[i]]
        lb > best_cost && continue

        for v in set
            insert_cost = dist[v1, v] + dist[v, tour[i]] - dist[v1, tour[i]]
            noise > 0.0 && (
                insert_cost +=
                    round(Int64, noise * rand(RNGs[Threads.threadid()]) * abs(insert_cost))
            )
            if insert_cost < best_cost
                best_cost = insert_cost
                bestv = v
                bestpos = i
            end
        end
    end
    return bestv, bestpos
end

@inline function insert_subset_lb(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    dist::Array{Int64, 2},
    set::Array{Int64, 1},
    setind::Int,
    setdist::Distsv,
    noise::Float64,
)
    best_cost = typemax(Int64)
    bestv = 0
    bestpos = 0
    tour_inds = collect(1:length(tour))

    for i in 1:ceil(Int64, length(tour) * noise)
        i = incremental_shuffle!(RNGs, tour_inds, i)
        v1 = prev_tour(tour, i)
        lb =
            setdist.vert_set[v1, setind] + setdist.set_vert[setind, tour[i]] -
            dist[v1, tour[i]]
        lb > best_cost && continue

        for v in set
            insert_cost = dist[v1, v] + dist[v, tour[i]] - dist[v1, tour[i]]
            if insert_cost < best_cost
                best_cost = insert_cost
                bestv = v
                bestpos = i
            end
        end
    end
    return bestv, bestpos
end

############ Initial Tour Construction ##########################

"""build tour from scratch on a cold restart"""
function initial_tour!(
    RNGs::Array{MersenneTwister, 1},
    lowest::Tour,
    dist::Array{Int64, 2},
    sets::Array{Any, 1},
    order_constraints::Array{Constraints, 1},
    setdist::Distsv,
    trial_num::Int64,
    member::Array{Int64, 1},
    param::Dict{Symbol, Any},
    start_set::Int64,
)
    sets_to_insert = collect(1:param[:num_sets])
    best = Tour(Int64[], typemax(Int64))

    # compute random initial tour only if past first trial
    # in this case, randomly choose between random and insertion tour.
    if param[:init_tour] == "rand" &&
       (trial_num > 1) &&
       (rand(RNGs[Threads.threadid()]) < 0.5)
        random_initial_tour!(
            RNGs,
            best.tour,
            sets_to_insert,
            dist,
            sets,
            order_constraints,
            member,
            start_set,
        )
    else
        random_insertion!(
            RNGs,
            best.tour,
            sets_to_insert,
            dist,
            sets,
            order_constraints,
            setdist,
            member,
            start_set,
        )
    end
    best.cost = tour_cost(best.tour, dist)
    lowest.cost > best.cost && (lowest = best)
    return best
end

"""
Randomly shuffle the sets, and then insert the best vertex from each set back into
the tour where sets are considered in shuffled order.
"""
function random_insertion!(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    sets_to_insert::Array{Int64, 1},
    dist::Array{Int64, 2},
    sets::Array{Any, 1},
    order_constraints::Array{Constraints, 1},
    setdist::Distsv,
    member::Array{Int64, 1},
    start_set::Int64,
)
    # Randomly permute the sets.
    shuffle!(RNGs[Threads.threadid()], sets_to_insert)
    pushed_cnt = 0
    to_insert = copy(sets_to_insert)
    while !isempty(to_insert)
        pushed_cnt += 1

        if pushed_cnt == 1
            rnd_set_idx = findfirst(x -> x == start_set, to_insert)
            rnd_set = to_insert[rnd_set_idx]
        else
            rnd_set_idx = rand(RNGs[Threads.threadid()], 1:length(to_insert))
            rnd_set = to_insert[rnd_set_idx]
        end

        min_insert_idx, max_insert_idx =
            calc_bounds(tour, rnd_set, order_constraints, member)

        if pushed_cnt > 1
            min_insert_idx = max(min_insert_idx, 2)
            if max_insert_idx < min_insert_idx
                max_insert_idx = min_insert_idx
            end
        end

        # Only have to compute the insert cost for the changed portion of the tour.
        if min_insert_idx > length(tour)
            best_vertex = rand(RNGs[Threads.threadid()], sets[rnd_set])
            best_position = length(tour) + 1
        elseif isempty(tour)
            best_vertex = rand(RNGs[Threads.threadid()], sets[rnd_set])
            best_position = 1
        else
            best_vertex, best_position = insert_lb(
                RNGs,
                tour,
                dist,
                sets[rnd_set],
                rnd_set,
                min_insert_idx,
                max_insert_idx,
                setdist,
                0.75,
            )
        end

        # Now, perform the insertion.
        insert!(tour, best_position, best_vertex)
        deleteat!(to_insert, rnd_set_idx)
    end
end

"""
Randomly shuffle the sets, and then insert the best vertex from each set back into
the tour where sets are considered in shuffled order.
"""
function random_initial_tour!(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    sets_to_insert::Array{Int64, 1},
    dist::Array{Int64, 2},
    sets::Array{Any, 1},
    order_constraints::Array{Constraints, 1},
    member::Array{Int64, 1},
    start_set::Int64,
)
    shuffle!(RNGs[Threads.threadid()], sets_to_insert)
    pushed_cnt = 0
    to_insert = copy(sets_to_insert)
    while !isempty(to_insert)
        pushed_cnt += 1

        if pushed_cnt == 1
            rnd_set_idx = findfirst(x -> x == start_set, to_insert)
            rnd_set = to_insert[rnd_set_idx]
        else
            rnd_set_idx = rand(RNGs[Threads.threadid()], 1:length(to_insert))
            rnd_set = to_insert[rnd_set_idx]
        end

        min_insert_idx, max_insert_idx =
            calc_bounds(tour, rnd_set, order_constraints, member)

        if pushed_cnt > 1
            min_insert_idx = max(min_insert_idx, 2)
            if max_insert_idx < min_insert_idx
                max_insert_idx = min_insert_idx
            end
        end

        idx_to_insert = rand(RNGs[Threads.threadid()], min_insert_idx:max_insert_idx)
        vert_to_insert = rand(RNGs[Threads.threadid()], sets[rnd_set])
        insert!(tour, idx_to_insert, vert_to_insert)
        deleteat!(to_insert, rnd_set_idx)
    end
end

######################### Removals ################################
"""
Remove the vertices randomly, but biased towards those that add the most length to the
tour.  Bias is based on the power input.  Vertices are then selected via pdf select.
"""
function worst_removal!(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    dist::Array{Int64, 2},
    num_to_remove::Int64,
    member::Array{Int64, 1},
    power::Float64,
    start_set::Int64,
)
    deleted_sets = Array{Int}(undef, 0)
    while length(deleted_sets) < num_to_remove
        removal_costs = worst_vertices(tour, dist)
        ind = pdf_select(RNGs, removal_costs[2:length(removal_costs)], power) + 1
        set_to_delete = member[tour[ind]]

        # perform the deletion
        push!(deleted_sets, set_to_delete)
        splice!(tour, ind)
    end
    return deleted_sets
end

""" removing a single continuos segment of the tour of size num_remove """
function segment_removal!(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    num_to_remove::Int64,
    member::Array{Int64, 1},
    start_set::Int64,
)
    i = rand(RNGs[Threads.threadid()], 2:length(tour))
    deleted_sets = Array{Int}(undef, 0)
    while length(deleted_sets) < num_to_remove
        i > length(tour) && (i = 2)
        set_to_delete = member[tour[i]]
        if set_to_delete == start_set
            continue
        end
        push!(deleted_sets, set_to_delete)
        splice!(tour, i)
    end
    return deleted_sets
end

"""  pick a random vertex, and delete its closest neighbors  """
function distance_removal!(
    RNGs::Array{MersenneTwister, 1},
    tour::Array{Int64, 1},
    dist::Array{Int64, 2},
    num_to_remove::Int64,
    member::Array{Int64, 1},
    power::Float64,
    start_set::Int64,
)
    deleted_sets = Array{Int}(undef, 0)
    deleted_vertices = Array{Int}(undef, 0)

    seed_index = rand(RNGs[Threads.threadid()], 2:length(tour))
    push!(deleted_sets, member[tour[seed_index]])
    push!(deleted_vertices, tour[seed_index])
    splice!(tour, seed_index)

    while length(deleted_sets) < num_to_remove
        # pick a random vertex from the set of deleted vertices
        seed_vertex = rand(RNGs[Threads.threadid()], deleted_vertices)
        # find closest vertex to the seed vertex that's still in the tour
        mindist = zeros(Int64, length(tour) - 1)
        for i in 2:length(tour)
            mindist[i - 1] = min(dist[seed_vertex, tour[i]], dist[tour[i], seed_vertex])
        end
        del_index = pdf_select(RNGs, mindist, power) + 1
        set_to_delete = member[tour[del_index]]
        if set_to_delete == start_set
            continue
        end
        push!(deleted_sets, set_to_delete)
        push!(deleted_vertices, tour[del_index])
        splice!(tour, del_index)
    end

    return deleted_sets
end

"""
determine the cost of removing each vertex from the tour, given that all others remain.
"""
function worst_vertices(tour::Array{Int64, 1}, dist::Array{Int64, 2})
    removal_cost = zeros(Int64, length(tour))
    for i in 2:length(tour)
        if i == length(tour)
            removal_cost[i] =
                dist[tour[i - 1], tour[i]] + dist[tour[i], tour[1]] -
                dist[tour[i - 1], tour[1]]
        else
            removal_cost[i] =
                dist[tour[i - 1], tour[i]] + dist[tour[i], tour[i + 1]] -
                dist[tour[i - 1], tour[i + 1]]
        end
    end
    return removal_cost
end

max_stars_repo_name: AndreiKud/PCGLNS

max_issues_repo_name: AndreiKud/PCGLNS

max_forks_repo_name: AndreiKud/PCGLNS


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module PolicyDefinitionsClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_PolicyDefinitionsApi.jl")

# export models
export convert, CloudError
export convert, ErrorAdditionalInfo
export convert, ErrorResponse
export convert, ParameterDefinitions
export convert, ParameterDefinitionsValue
export convert, PolicyDefinition
export convert, PolicyDefinitionListResult
export convert, PolicyDefinitionProperties

# export operations
export convert, PolicyDefinitionsApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
module JavaVM

export init, destroy

using JavaCall.JNI
using JavaCall.InitOptions
using JavaCall.Utils

const JAVA_HOME_CANDIDATES = [
    "/usr/lib/jvm/default-java/",
    "/usr/lib/jvm/default/"
]

init() = init(defaultopts())

init(opts::JavaCallInitOptions) = init(Val(opts.fromcurrentvm), opts)

function init(::Val{false}, opts::JavaCallInitOptions)
    @verify !jvmloaded() "JVM already initialised."
    fullopts = []
    push!(fullopts, "-Djava.class.path="*join(opts.jvmclasspath, ":"))
    push!(fullopts, opts.jvmopts...)
    JNI.init_new_vm(findjvm(opts.javahome), fullopts)
end

function init(::Val{true}, opts::JavaCallInitOptions)
    @verify !jvmloaded() "JVM already initialised."
    JNI.init_current_vm(findjvm(opts.javahome))
end

function destroy()
    @verify jvmloaded() "JVM not initialised."
    JNI.destroy_vm()
end

jvmloaded() = JNI.is_jni_loaded() && JNI.is_env_loaded()

function findjvm(javahome::Union{String, Nothing})
    try
        return possible_javahomes(javahome) |>
            l -> flatmap(libpathdirectories, l) |>
            # Include pwd in directories to search for libpath
            l -> chain((pwd(),), l) |>
            l -> flatmap(libpathfromdirectory, l)  |>
            # Return only first path
            first
    catch
        errormsg = [
            "Cannot find java library $(libfile)\n",
            "Search Path:"
        ]
        for path in possible_javahomes(javahome) |> l -> flatmap(libpathdirectories, l) |> l -> chain((pwd(),), l)
            push!(errormsg,"\n   $path")
        end
        error(reduce(*,errormsg))
    end
end

possible_javahomes(javahome::String) = [javahome]
possible_javahomes(::Nothing) = possible_javahomes()

@static if Sys.isunix()
    function possible_javahomes()
        javahomes = Any[]
        if haskey(ENV, "JAVA_HOME")
            push!(javahomes, ENV["JAVA_HOME"])
        else
            # Find default javahome by checking location of the java command
            try
                javapath = chomp(read(`which java`,String))
                while(islink(javapath))
                    javapath = readlink(javapath)
                end
                javapath = dirname(javapath)
                javapath = match(r"(.*)((/jre/bin)|(/bin))+",javapath)[1]
                push!(javahomes,javapath)
            catch err
                @debug "JavaCall could not determine javapath from `which java`" err
            end
        end
        isfile("/usr/libexec/java_home") && push!(javahomes, chomp(read(`/usr/libexec/java_home`, String)))

        for fname ∈ JAVA_HOME_CANDIDATES
            isdir(fname) && push!(javahomes, fname)
        end
        javahomes
    end

    @static if Sys.isapple()
        libfile = "libjvm.dylib"

        function libpathdirectories(java_home::AbstractString)::Vector{AbstractString}
            [
                joinpath(java_home, "jre", "lib", "server"),
                joinpath(java_home, "lib", "server")
            ]
        end

        function libpathfromdirectory(directory::AbstractString)::Vector{AbstractString}
            libpath = joinpath(directory, libfile)
            return (isfile(libpath) ? [libpath] : [])
        end
    else
        libfile = "libjvm.so"

        function libpathdirectories(java_home::AbstractString)::Vector{AbstractString}
            libpaths = []
            if Sys.WORD_SIZE==64
                push!(libpaths, joinpath(java_home, "jre", "lib", "amd64", "server"))
                push!(libpaths, joinpath(java_home, "lib", "amd64", "server"))
            elseif Sys.WORD_SIZE==32
                push!(libpaths, joinpath(java_home, "jre", "lib", "i386", "server"))

                push!(libpaths, joinpath(java_home, "lib", "i386", "server"))
            end
            push!(libpaths, joinpath(java_home, "jre", "lib", "server"))
            push!(libpaths, joinpath(java_home, "lib", "server"))
            libpaths
        end

        function libpathfromdirectory(directory::AbstractString)::Vector{AbstractString}
            libpath = joinpath(directory,libfile)
            return (isfile(libpath) ? [libpath] : [])
        end
    end
else
    libfile = "jvm.dll"

    function javahome_winreg()
        keys = ["SOFTWARE\\JavaSoft\\Java Runtime Environment", "SOFTWARE\\JavaSoft\\Java Development Kit", "SOFTWARE\\JavaSoft\\JRE", "SOFTWARE\\JavaSoft\\JDK"]

        for key in keys
            try
                value = querykey(WinReg.HKEY_LOCAL_MACHINE, key, "CurrentVersion")
                key *= "\\" * value
                return querykey(WinReg.HKEY_LOCAL_MACHINE, key, "JavaHome")
            catch
                # Try the next value in the loop.
            end
        end

        error("Cannot find an installation of Java in the Windows Registry. Please install a JRE/JDK, or set the JAVA_HOME environment variable if one is already installed.")
    end

    function possible_javahomes()
        javahomes = Any[]
        if haskey(ENV, "JAVA_HOME")
            push!(javahomes, ENV["JAVA_HOME"])
        else
            ENV["JAVA_HOME"] = javahome_winreg()
            push!(javahomes, ENV["JAVA_HOME"])
        end
        isfile("/usr/libexec/java_home") && push!(javahomes, chomp(read(`/usr/libexec/java_home`, String)))

        for fname ∈ JAVA_HOME_CANDIDATES
            isdir(fname) && push!(javahomes, fname)
        end
        javahomes
    end

    function libpathdirectories(java_home::AbstractString)::Vector{AbstractString}
        [
            joinpath(java_home, "bin", "server"),
            joinpath(java_home, "jre", "bin", "server"),
            joinpath(java_home, "bin", "client"),
            joinpath(java_home, "jre", "lib", "server"),
            joinpath(java_home, "lib", "server")
        ]
    end

    function libpathfromdirectory(directory::AbstractString)::Vector{AbstractString}
        libpath = joinpath(directory, libfile)
        if isfile(libpath)
            bindir = dirname(dirname(libpath))
            m = filter(x -> occursin(r"msvcr(?:.*).dll",x), readdir(bindir))
            if !isempty(m)
                return [joinpath(bindir,m[1]),libpath]
            end
        end
        return []
    end
end

end


text:
using DataStructures: Deque

@testset "Deque" begin
a = Deque{Float64}()

@test isempty(a)
@test length(a) == 0

push!(a, 13)

@test !isempty(a)
@test length(a) == 1

end


text:
# Testing the model-building interface

# All models should have:
# 1. An initial distribution, μ0(θ)
# 2. A transition distribution, f(x,t,θ)
# 3. An observation distribution g(x,t,θ)
# 4. A proposal distribution q(x,y,t,θ)
# 5. A weight function ξ(x,y,t,θ)

δ = 1.0
φ = 0.95
ση = 1.0
σε = 1.0

θ0 = [δ;atanh(φ);ση;σε]

β = [0.0;1.0]

N = 1000

@testset "LinearStateSpaceModel" begin
    @testset "Generic model interface" begin
        arh = ARH(β)
        @test typeof(arh) == LinearStateSpaceModel

        μ0 = initial_distribution(arh,θ0)

        @test mean(μ0) == [δ/(1-φ)]
        @test cov(μ0) == fill(ση^2/(1-φ^2),1,1)

        x0 = rand(μ0)

        u1 = [1.0]

        f  = transition_distribution(arh,x0,u1,1,θ0)

        @test mean(f) == φ*x0 .+ δ
        @test cov(f) == fill(ση^2,1,1)

        x1 = rand(f)

        g  = observation_distribution(arh,x1,u1,1,θ0)

        @test mean(g) == x1
        @test cov(g) == fill(σε^2,1,1)

        y1 = rand(g)

        # TODO: Proposal distributions for LinearStateSpaceModels (#4)
        q = proposal_distribution(arh,x0,y1,u1,1,θ0)

        @test mean(q) == φ*x0 .+ δ
        @test cov(q) == fill(ση^2,1,1)

        ξ = weight_function(arh,x0,y1,u1,1,θ0)

        @test ξ == 0.0
    end

    @testset "Linear model interface" begin
        arh = ARH([0.0;1.0])

        μ = initial_mean(arh,θ0)

        @test μ == fill(δ/(1-φ),1)

        Σ0 = initial_covariance(arh,θ0)

        @test Σ0 == fill(ση^2/(1-φ^2),1,1)

        F = state_transition_matrix(arh,θ0)

        @test F == fill(φ,1,1)

        Q = state_covariance(arh,θ0)

        @test Q == fill(ση^2,1,1)

        G = state_input_matrix(arh,θ0)

        @test G == fill(δ,1,1)

        H = observation_matrix(arh,θ0)

        @test H == fill(β[2],1,1)

        R = observation_covariance(arh,θ0)

        @test R == fill(σε^2,1,1)

        Γ = observation_input_matrix(arh,θ0)

        @test Γ == fill(β[1],1,1)

        Θ = state_noise_transformation_matrix(arh,θ0)

        @test Θ == fill(1.0,1,1)
    end
end

@testset "Kalman filtering and smoothing" begin
    arh = ARH(β)

    X,Y = Filtering.simulate(arh,θ0,1000,fill([1.0],N))

    @test length(X) == N+1
    @test length(Y) == N

    kf = kalman_filter(arh,θ0,Y,fill([1.0],length(Y)))

    ks = kalman_smoother(arh,θ0,kf)
end

@testset "Particle filtering and smoothing" begin
    arh = ARH(β)

    X,Y = Filtering.simulate(arh,θ0,N,fill([1.0],N))

    pf = particle_filter(arh,θ0,Y,100,fill([1.0],N))

    B = 2
    ps = particle_smoother(arh,θ0,Y,pf,B,fill([1.0],N))

    @test size(ps) == (1,N+1,B)
end


text:
using LinearAlgebra, SparseArrays, SuiteSparse
mutable struct LinSolveFactorize{F}
  factorization::F
  A
end
LinSolveFactorize(factorization) = LinSolveFactorize(factorization,nothing)
function (p::LinSolveFactorize)(x,A,b,update_matrix=false;kwargs...)
  if update_matrix
    p.A = p.factorization(A)
  end
  if typeof(p.A) <: SuiteSparse.UMFPACK.UmfpackLU || typeof(p.factorization) <: typeof(lu)
    ldiv!(x,p.A,b) # No 2-arg form for SparseArrays!
  else
    x .= b
    ldiv!(p.A,x)
  end
end
function (p::LinSolveFactorize)(::Type{Val{:init}},f,u0_prototype)
  LinSolveFactorize(p.factorization,nothing)
end

mutable struct LinSolveGPUFactorize{F,T}
  factorization::F
  A
  x_cache::T
end
LinSolveGPUFactorize(factorization=qr) = LinSolveGPUFactorize(factorization,nothing,nothing)
function (p::LinSolveGPUFactorize)(x,A,b,update_matrix=false;kwargs...)
  if update_matrix
    p.A = p.factorization(cuify(A))
  end
  ldiv!(p.x_cache,p.A,cuify(b))
  x .= Array(p.x_cache)
end
function (p::LinSolveGPUFactorize)(::Type{Val{:init}},f,u0_prototype)
  LinSolveGPUFactorize(p.factorization,nothing,cuify(u0_prototype))
end

### Default Linsolve

# Try to be as smart as possible
# lu! if Matrix
# lu if sparse
# gmres if operator

mutable struct DefaultLinSolve
  A
  iterable
end
DefaultLinSolve() = DefaultLinSolve(nothing, nothing)

function (p::DefaultLinSolve)(x,A,b,update_matrix=false;tol=nothing, kwargs...)
  if p.iterable isa Vector && eltype(p.iterable) <: LinearAlgebra.BlasInt # `iterable` here is the pivoting vector
    F = LU{eltype(A)}(A, p.iterable, zero(LinearAlgebra.BlasInt))
    ldiv!(x, F, b)
    return nothing
  end
  if update_matrix
    if typeof(A) <: Matrix
      blasvendor = BLAS.vendor()
      if (blasvendor === :openblas || blasvendor === :openblas64) && size(A,1) <= 500 && ArrayInterface.can_setindex(x) # if the user doesn't use OpenBLAS, we assume that is a much better BLAS implementation like MKL
        p.A = RecursiveFactorization.lu!(A)
      else
        p.A = lu!(A)
      end
    elseif typeof(A) <: Tridiagonal
      p.A = lu!(A)
    elseif typeof(A) <: Union{SymTridiagonal}
      p.A = ldlt!(A)
    elseif typeof(A) <: Union{Symmetric,Hermitian}
      p.A = bunchkaufman!(A)
    elseif typeof(A) <: SparseMatrixCSC
      p.A = lu(A)
    elseif ArrayInterface.isstructured(A)
      p.A = factorize(A)
    elseif !(typeof(A) <: AbstractDiffEqOperator)
      # Most likely QR is the one that is overloaded
      # Works on things like CuArrays
      p.A = qr(A)
    end
  end

  if typeof(A) <: Union{Matrix,SymTridiagonal,Tridiagonal,Symmetric,Hermitian} # No 2-arg form for SparseArrays!
    x .= b
    ldiv!(p.A,x)
  # Missing a little bit of efficiency in a rare case
  #elseif typeof(A) <: DiffEqArrayOperator
  #  ldiv!(x,p.A,b)
  elseif ArrayInterface.isstructured(A) || A isa SparseMatrixCSC
    ldiv!(x,p.A,b)
  elseif typeof(A) <: AbstractDiffEqOperator
    # No good starting guess, so guess zero
    if p.iterable === nothing
      p.iterable = IterativeSolvers.gmres_iterable!(x,A,b;initially_zero=true,restart=5,maxiter=5,tol=1e-16,kwargs...)
      p.iterable.reltol = tol
    end
    x .= false
    iter = p.iterable
    purge_history!(iter, x, b)

    for residual in iter
    end
  else
    ldiv!(x,p.A,b)
  end
  return nothing
end

function (p::DefaultLinSolve)(::Type{Val{:init}},f,u0_prototype)
  if has_Wfact(f) || has_Wfact_t(f)
    piv = collect(one(LinearAlgebra.BlasInt):convert(LinearAlgebra.BlasInt, length(u0_prototype))) # pivoting vector
    DefaultLinSolve(f, piv)
  else
    DefaultLinSolve()
  end
end

const DEFAULT_LINSOLVE = DefaultLinSolve()

### Default GMRES

# Easily change to GMRES

mutable struct LinSolveIterativeSolvers{F,PL,PR,AR,A}
  generate_iterator::F
  iterable
  Pl::PL
  Pr::PR
  args::AR
  kwargs::A
end
LinSolveIterativeSolvers(generate_iterator,args...;
                         Pl=IterativeSolvers.Identity(),
                         Pr=IterativeSolvers.Identity(),
                         kwargs...) =
                         LinSolveIterativeSolvers(generate_iterator, nothing, Pl, Pr, args, kwargs)

LinSolveGMRES(args...;kwargs...) = LinSolveIterativeSolvers(IterativeSolvers.gmres_iterable!,args...;kwargs...)
LinSolveCG(args...;kwargs...) = LinSolveIterativeSolvers(IterativeSolvers.cg_iterator!,args...;kwargs...)
LinSolveBiCGStabl(args...;kwargs...) = LinSolveIterativeSolvers(IterativeSolvers.bicgstabl_iterator!,args...;kwargs...)
LinSolveChebyshev(args...;kwargs...) = LinSolveIterativeSolvers(IterativeSolvers.chebyshev_iterable!,args...;kwargs...)
LinSolveMINRES(args...;kwargs...) = LinSolveIterativeSolvers(IterativeSolvers.minres_iterable!,args...;kwargs...)

function (f::LinSolveIterativeSolvers)(x,A,b,update_matrix=false; Pl=nothing, Pr=nothing, tol=nothing, kwargs...)
  if f.iterable === nothing
    Pl = ComposePreconditioner(f.Pl, Pl, true)
    Pr = ComposePreconditioner(f.Pr, Pr, false)
    f.iterable = f.generate_iterator(x,A,b,f.args...;
                                     initially_zero=true,restart=5,
                                     maxiter=5,tol=1e-16,Pl=Pl,Pr=Pr,
                                     f.kwargs...,kwargs...)
    tol′ = get(f.kwargs, :tol, nothing)
    tol′ !== nothing && (tol = tol′)
    f.iterable.reltol = tol
  end
  x .= false
  iter = f.iterable
  purge_history!(iter, x, b)

  for residual in iter
  end
  return nothing
end

function (p::LinSolveIterativeSolvers)(::Type{Val{:init}},f,u0_prototype)
  LinSolveIterativeSolvers(p.generate_iterator, nothing, p.Pl, p.Pr, p.args, p.kwargs)
end

# scaling for iterative solvers
struct ScaleVector{T}
  x::T
  isleft::Bool
end
function LinearAlgebra.ldiv!(v::ScaleVector, x)
  vx = vec(v.x)
  if v.isleft
    return @.. x = x * vx
  else
    return @.. x = x / vx
  end
end
function LinearAlgebra.ldiv!(y, v::ScaleVector, x)
  vx = vec(v.x)
  if v.isleft
    return @.. y = x * vx
  else
    return @.. y = x / vx
  end
end

struct ComposePreconditioner{P, S<:Union{ScaleVector,Nothing}}
  P::P
  scale::S
  isleft::Bool
end

function LinearAlgebra.ldiv!(v::ComposePreconditioner, x)
  isid = v.P isa IterativeSolvers.Identity
  isid || ldiv!(v.P, x)
  s = v.scale
  s === nothing && return x
  ldiv!(s, x)
  return x
end

function LinearAlgebra.ldiv!(y, v::ComposePreconditioner, x)
  isid = v.P isa IterativeSolvers.Identity
  isid || ldiv!(y, v.P, x)
  s = v.scale
  s === nothing && return x
  if isid
    ldiv!(y, s, x)
  else
    if v.isleft
      @.. y = y * s.x
    else
      @.. y = y / s.x
    end
  end
  return y
end

function purge_history!(iter::IterativeSolvers.GMRESIterable, x, b)
  iter.k = 1
  iter.x  = x
  iter.b  = b

  iter.residual.current = IterativeSolvers.init!(iter.arnoldi, iter.x, iter.b, iter.Pl, iter.Ax, initially_zero = true)
  IterativeSolvers.init_residual!(iter.residual, iter.residual.current)
  iter.β = iter.residual.current
  nothing
end


text:
using DashBootstrapComponents, DashHtmlComponents

button_group = html_div(
    [
        dbc_radioitems(
            id = "radios",
            className = "btn-group",
            inputClassName = "btn-check",
            labelClassName = "btn btn-outline-primary",
            labelCheckedClassName = "active",
            options = [
                Dict("label" => "Option 1", "value" => 1),
                Dict("label" => "Option 2", "value" => 2),
                Dict("label" => "Option 3", "value" => 3),
            ],
            value = 1,
        ),
        html_div(id = "output"),
    ],
    className = "radio-group",
);

callback!(app, Output("output", "children"), Input("radios", "value")) do value
    return "Selected value: $value"
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module LyapunovExponents

export LEProblem, ContinuousLEProblem, DiscreteLEProblem, lyapunov_exponents,
    phase_state, report

# Re-export methods from DifferentialEquations extended here:
export init, solve, solve!, step!
import DiffEqBase: init, step!
using DiffEqBase: solve!, set_u!
# solve and solve! are imported/extended in stages.jl
using LinearAlgebra
using Printf
using Markdown

include("documents.jl")
include("stages.jl")
include("utils.jl")
include("online_stats.jl")
include("types.jl")
include("core.jl")
include("terminators.jl")
include("coevolve.jl")
include("continuous_exponents.jl")
include("discrete_exponents.jl")
include("rode.jl")
include("tangent_utils.jl")
include("clv/covariant_vectors.jl")
include("interface.jl")
include("implementations.jl")
include("reports.jl")
include("examples/examples.jl")
include("test.jl")
include("diffeq_hack.jl")
include("plotting.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Imath"
version = v"3.0.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/AcademySoftwareFoundation/Imath/archive/refs/tags/v3.0.1.tar.gz", "9cd984bb6b0a9572dd4a373b1fab60bc4c992a52ec5c68328fe0f48f194ba3c0")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Imath*
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DBUILD_TESTING=OFF \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
make -j${nproc}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libImath-3_0", :libImath)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
extend(f::typeof(+), ::Type{Val{:name}}) = nameof(f)

function Base.getproperty(f::typeof(+), name::Symbol)
    extend(f, Val{name})
end


text:
# Create the agent, model, and grid types.
"AbstractModel type for the Schelling Model"
mutable struct SchellingModel{T<:Integer, Y<:AbstractArray, Z<:AbstractSpace} <: AbstractModel
  # Object should always be a subtype of AbstractModel.
  "A field of the model for a space object, always a subtype of AbstractSpace."
  space::Z
  "A list of agents."
  agents::Y
  "A field for the scheduler function."
  scheduler::Function
  "The minimum number of neighbors for agent to be happy."
  min_to_be_happy::T
end

"AbstractAgent type for the Schelling Agent"
mutable struct SchellingAgent{T<:Integer} <: AbstractAgent # Object should always be a subtype of AbstractAgent
  "The identifier number of the agent."
  id::T
  "The x, y location of the agent."
  pos::Tuple{T, T}
  "Whether or not the agent is happy with cell (true is 'happy' and false is 'unhappy')."
  mood::Bool
  "The group of the agent, determines mood as it interacts with neighbors."
  group::T
end

"The space of the experiment."
mutable struct SGrid{T<:Integer, Y<:AbstractArray} <: AbstractSpace
  "Dimensions of the grid."
  dimensions::Tuple{T, T}
  "The space type."
  space::SimpleGraph
  "An array of arrays for each grid node."
  agent_positions::Y
end

"Function to instantiate the model."
function instantiate_modelS(;numagents=320, griddims=(20, 20), min_to_be_happy=3)
  # Create an empty Array of Arrays.
  agent_positions = [Int64[] for i in 1:gridsize(griddims)]
  # Use SGrid to create a grid from griddims and agent_positions using the grid function.
  # Create a 2D grid with nodes have a max of 8 neighbors.
  mygrid = SGrid(griddims, grid(griddims, false, true), agent_positions)
  # Instantiate the model using mygrid, the SchellingAgent type, the random_activation function from Agents.jl
  # and the argument min_to_be_happy.
  model = SchellingModel(mygrid, SchellingAgent[], random_activation, min_to_be_happy)

  # Create a 1-dimension list of agents, balanced evenly between group 0 and group 1.
  agents = vcat(
    [SchellingAgent(Int(i), (1,1), false, 0) for i in 1:(numagents/2)],
    [SchellingAgent(Int(i), (1,1), false, 1) for i in (numagents/2)+1:numagents]
  )

  # Add the agents to the model.
  for agent in agents
    # Use add_agent_single (from Agents.jl) to add the agents to the grid at random locations.
    add_agent_single!(agent, model)
  end
  return model
end

"Move a single agent until a satisfactory location is found."
function agent_step!(agent, model)
  if agent.mood == true
    return
  end
  while agent.mood == false
    neighbor_cells = node_neighbors(agent, model)
    count_neighbors_same_group = 0

	# For each neighbor, get group and compare to current agent's group...
    # ...and increment count_neighbors_same_group as appropriately.
    for neighbor_cell in neighbor_cells
      node_contents = get_node_contents(neighbor_cell, model)
      # Skip iteration if the node is empty.
      if length(node_contents) == 0
        continue
      else
        # Otherwise, get the first agent in the node...
        node_contents = node_contents[1]
      end
      # ...and increment count_neighbors_same_group if the neighbor's group is the same.
      neighbor_agent_group = model.agents[node_contents].group
      if neighbor_agent_group == agent.group
        count_neighbors_same_group += 1
      end
    end

    # After evaluating and adding up the groups of the neighbors, decide whether or not to move the agent.
    # If count_neighbors_same_group is at least the min_to_be_happy, set the mood to true. Otherwise,
	# move the agent using move_agent_single.
    if count_neighbors_same_group >= model.min_to_be_happy
      agent.mood = true
    else
      move_agent_single!(agent, model)
    end
  end
end


text:
using DataFrames
using Random
using LinearAlgebra
using FixedEffectModels

Random.seed!(123)
N = 2
T = 100_000
x = randn(T,N)
β = [1 0; -1 2; 0.5 0.6 ; 0.3 0; 0.2 1]
y = zeros(Float64, T)
for t in 1:T,l in 0:min(4,t-1)
    y[t] += dot(β[l+1, :], x[t-l, :]) + 0.01*randn()
end

# Make the dataframe
df = DataFrame(y=y)
df[!, :x1] = x[:, 1]
df[!, :x2] = x[:, 2]

# Make Lags
L!(df, :x1, 1:5)
L!(df, :x2, 1:5)

Term.(vcat([:x1, :x2], [Symbol("L",i,"x1") for i in 1:5], [Symbol("L",i,"x2") for i in 1:5]))
f = Term(:y) ~ foldl(+, Term.(vcat([:x1, :x2], [Symbol("L",i,"x1") for i in 1:5], [Symbol("L",i,"x2") for i in 1:5])))

lp = LocalProjection(;df=df, horizons=1:10, formula = f, regstyle=:levels)

makedepvars(lp)
runregs(lp)

@test all(round.(LocalProjections.coef(lp, :x1); digits=1)[1:4] .== β[2:end,1])
@test all(round.(LocalProjections.coef(lp, :x2); digits=1)[1:4] .== β[2:end,2])