text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# ╔═╡ 53c06a7e-4f02-11ec-1109-65224850d083
begin
	using Pkg
	Pkg.add( url="/home/graham_s/julia/vw/ScottishTaxBenefitModel" )
end

# ╔═╡ 3dce1b1e-9756-48f2-a7dc-45363d1e3995
begin
	using BudgetConstraints
	using SurveyDataWeighting
	using PovertyAndInequalityMeasures
	using Plots,CSV,DataFrames,PlutoUI
end

# ╔═╡ b6b22dd4-bbb3-4c29-a27f-521cf14fc1d6
begin
	using ScottishTaxBenefitModel
	import .ExampleHouseholdGetter
	using .STBParameters
	using .BCCalcs
	using .ModelHousehold
	import .FRSHouseholdGetter
	using .Utils
	using .Definitions
	using .SingleHouseholdCalculations
	using .RunSettings
end

# ╔═╡ 9ba57d28-9f19-43c6-b64b-79416ca8a12a
begin

	PlutoUI.TableOfContents(aside=true)

end

# ╔═╡ 38beba12-681c-4d02-a247-8388e0b7df9b
md"""

# A Tax Benefit Model for Scotland

Graham Stark [gks56@open.ac.uk](gks56@open.ac.uk)/[graham.stark@virtual-worlds.biz](PlutoUI.TableOfContents(aside=true))

"""

# ╔═╡ 9250ad90-8270-468c-9b51-87fc08451e08

begin
	# housekeeping stuff
	const DEFAULT_NUM_TYPE = Float64
	settings = RunSettings.Settings()
	settings.requested_threads = 4
	settings.data_dir = "/home/graham_s/julia/vw/ScottishTaxBenefitModel/data"

	function init_data(; reset :: Bool = false )
	   nhh = FRSHouseholdGetter.get_num_households()
	   num_people = -1
	   if( nhh == 0 ) || reset
		  @time nhh, num_people,nhh2 = FRSHouseholdGetter.initialise( settings )
	   end
	   (nhh,num_people)
	end

	function load_system()::TaxBenefitSystem
		sys = load_file( joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021_22.jl" ))
		#
		# Note that as of Budget21 removing these doesn't actually happen till May 2022.
		#
		load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021-uplift-removed.jl"))
		# uc taper to 55
		load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "budget_2021_uc_changes.jl"))
		weeklyise!( sys )
		return sys
	end

md"""
"""

end

# ╔═╡ e8fa154d-e60b-407e-84ea-2b247141374d
begin
	ExampleHouseholdGetter.initialise(settings)
	nhhs,npeople = init_data()
	md"  "
end

# ╔═╡ b46666a3-3a24-485d-a1b6-41250fc21613
begin
hh = FRSHouseholdGetter.get_household(1)
for i in 1:nhhs
	hh =  FRSHouseholdGetter.get_household(i)

end
end

# ╔═╡ 2917cd0f-2869-41f9-8b65-62cb75f3d278
md"""

## What?

Scotben is an open source microsimulation model tax-benefit of Scotland. Brand new - only brought to a useable state in the last fortnight. This is its first outing.

"""

# ╔═╡ ed1ca53b-0b7c-47e6-b4c8-9167731a49fb
md"""
## Why?

* Independence vs the Union: what would the budget of an independent Scotland budget look like?
* Relationship with rUK is complicated;
* Standard of Debate in Scotland not high!
* I was bored! This was my lockdown distraction. I'd been thinking about building this for several years, though.

"""

# ╔═╡ 067fd91d-77c6-436e-880b-45772b27d17a
md"""
## How?
* Conventional Tax Benefit Model - similar in outline to e.g. IFS Taxben, IPPR
* Highly modular design: bolt things together to make e.g. forecasting models, social care sims etc.
* Static model - no non-takeup, labour supply etc. But get the structure right and these things are much easier;
* [Test first development](https://testing.googleblog.com/2007/01/introducing-testing-on-toilet.html) principle - write unit tests and then just enough code to make them pass - the [model test suite](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/tree/master/test) is almost as big as the model itself.

"""

# ╔═╡ 0dfaf005-b39a-44c3-8b51-005feac8d78b
md"""
## Julia

* [Julia](https://julialang.org) bridges the gap between statistics packages and high-level programming languages;
* you can do detailed, accurate programming with good support for [rich types](https://docs.julialang.org/en/v1/manual/types/) (vital in a tax benefit model);
* but also use it for [data science/econometrics](https://juliastats.org/) (all the model regressions and graphics are native julia);
* *very* fast - [outpacing C/Fortran in some cases](https://julialang.org/benchmarks/);
* [Julia has a huge collection of contributed packages](https://juliahub.com) - including standard econometrics routines, data handling, graphics, differential equations and much else;
* This presentation is itself a Julia program, written using [Pluto](https://github.com/fonsp/Pluto.jl). The model and its data are actually loaded directly in to this presentation and we'll interact with them very briefly in a minute.
"""

# ╔═╡ 2e9567f8-10f6-46d0-b32d-39df2038aaf7
begin
	x = 99
	x + 1

end

# ╔═╡ 444b957f-0bcb-4703-9477-dcfcdc320ca6
md"""
## Open Source

* all the program code is available on [GitHub](https://github.com/grahamstark/ScottishTaxBenefitModel.jl) under a [permissive licence](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/LICENSE)
* I can't hide! Much harder to fudge things this way.

"""

# ╔═╡ 9512f4de-46ac-4f85-84c9-7a37ab391143
md"""
## Structure

Divided into [packages](https://julialang.org/packages/) and [modules](https://docs.julialang.org/en/v1/manual/modules/). A package is a high level generic chunk of code that can be downloaded and used independently. A module is a namespace - a small chunk of package code in which you can hide messy details from the rest of the program.

#### High Level generic packages:

* [Budget Constraints](https://github.com/grahamstark/BudgetConstraints.jl);
* [Data Weighting](https://github.com/grahamstark/SurveyDataWeighting.jl);
* [Poverty & Inequality](https://github.com/grahamstark/PovertyAndInequalityMeasures.jl).

#### The Model

[The Model is itself a package](https://github.com/grahamstark/ScottishTaxBenefitModel.jl). Intermally, it's broken down into a collection of semi-independent modules, for example:

* [a household](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/src/ModelHousehold.jl);
* [the fiscal system parameters](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/src/STBParameters.jl);
* [means-tested benefits](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/src/LegacyMeansTestedBenefits.jl);
* [income tax](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/src/IncomeTaxCalculations.jl)

.. [and so on](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/tree/master/src)

Some of the modules (e.g. [Equivalence Scales](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/src/EquivalenceScales.jl)) may eventually be moved out into generic packages.

"""

# ╔═╡ 2334e293-57f7-4b80-ad50-f889c8f6bab0
md"""
## Data

* uses pooled (2015-2018) Scottish Households from [Family Resources Survey](https://www.gov.uk/government/collections/family-resources-survey--2) from the [UK Data Service](https://ukdataservice.ac.uk/);
* Data from the [Scottish Household Survey](https://www.gov.scot/collections/scottish-household-survey/) is [matched in](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/tree/master/matching). This allows much more accurate modelling of local taxes and (in the future) health, housing conditions;
* Since we have a [Calmar-like weighting system built-in](https://github.com/grahamstark/SurveyDataWeighting.jl), we can [weight to Scottish Population, Employment Levels, and so on](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/blob/master/src/Weighting.jl) very accurately and easily.
"""

# ╔═╡ d10d14c7-2b1a-4dde-8ea4-e5c8dcd86887
begin
	hhs = CSV.File( "$(settings.data_dir)/model_households_scotland.tab" )|>DataFrame
	people = CSV.File( "$(settings.data_dir)/model_people_scotland.tab" )|>DataFrame
end

# ╔═╡ abc6d492-562f-4186-abf7-70aa825648fe
md"""
## The Model In Action

No point if people can't use it. I've no plans to build a full web UI. [Been there, done that, boring and rarely used](https://www.virtual-worlds.scot/demonstrations/).

Instead, simple single-page UIs illustrating one point of interest. So far:

* [Budget Constraints](https://stb.virtual-worlds.scot/bcd/) - the weird, kinky world of the UK Tax Benefit System
* [A Simple Budget For Scotland](https://stb.virtual-worlds.scot/scotbudg/)

[These visualisations]() use [Dash](https://dash-julia.plotly.com/).

For full detailed interactions, this Pluto notebook system and 'lab assistant' programs such as [Dr Watson](https://juliadynamics.github.io/DrWatson.jl/dev/).

"""

# ╔═╡ 43eaea11-1764-4dee-af79-05aa264937c4
md"""

## Known Problems

This is a very new model and full results are just emerging over the last 2 weeks. The test suite gives me confidence that the low-level calculations are accurate. However:

* Modelled [Inequality is too low](https://data.gov.scot/poverty/) compared to official figures - possibly related to 100% takeup
* Revenue estimates for Income Tax ~£1bn too high compared to [official forecasts](https://www.fiscalcommission.scot/publications/scotlands-economic-and-fiscal-forecasts-august-2021/):
  - possibly employers pension contributions;
  - or something to do with the pandemic.

"""

# ╔═╡ 47148425-e812-4ec3-ab40-49603b5d53ee
md"""

## Next Steps

Scotben's clean interfaces and modular structure makes adding features easy. Julia also working as a statistics package makes adding behavioural features much easier.

* [Benefit non-takeup corrections](https://ifs.org.uk/publications/1954)  - this seems important for e.g. inequality estimates;
* long term projections (e.g. the [Scottish Growth Commission](https://www.sustainablegrowthcommission.scot/)- long term projection is technically quite easy with the components we have. [An earlier exercise is available](https://github.com/grahamstark/scottish_child_poverty_projections);
* consumption and spending data - either using estimated Engel Curves or matching in [Living Costs and Food Survey](https://www.ons.gov.uk/peoplepopulationandcommunity/personalandhouseholdfinances/incomeandwealth/methodologies/livingcostsandfoodsurvey) (LCF) data.

And much else.

"""

# ╔═╡ 285cdbf5-1e73-428a-bc1d-88ea8a731ce0
md"""
## To Find Out More

* **Tax Benefit Models**: [A short introduction to microsimulation and tax benefit models](https://stb.virtual-worlds.scot/intro.html). Originally written for the Open University, it covers all the essential ideas. | I've also written [the most boring blog ever about the Model](https://stb-blog.virtual-worlds.scot/);

* **Poverty and Inequality**: [My Notes](https://stb.virtual-worlds.scot/poverty.html) | [World Bank Handbook](http://documents.worldbank.org/curated/en/488081468157174849/Handbook-on-poverty-and-inequality) | [Official Figures for Scotland](https://data.gov.scot/poverty/);

* **Scotland's Finances**: [Scottish Fiscal Commission](https://www.fiscalcommission.scot/publications/scotlands-economic-and-fiscal-forecasts-august-2021/) | [Scottish Government Budget Documents](https://www.gov.scot/budget/).

I'd very much welcome contributions and suggestions. If you spot anything odd or if you have any ideas for how this can be improved, you can:

* [Open an issue on GitHub](https://github.com/grahamstark/ScottishTaxBenefitModel.jl/issues); or
* [email me](mailto:graham.stark@virtual-worlds.biz).

finally, you can download this presentation from [https://virtual-worlds.scot/ima2021/pres.zip]

"""

# ╔═╡ Cell order:
# ╟─53c06a7e-4f02-11ec-1109-65224850d083
# ╟─3dce1b1e-9756-48f2-a7dc-45363d1e3995
# ╟─b6b22dd4-bbb3-4c29-a27f-521cf14fc1d6
# ╟─9ba57d28-9f19-43c6-b64b-79416ca8a12a
# ╟─38beba12-681c-4d02-a247-8388e0b7df9b
# ╟─9250ad90-8270-468c-9b51-87fc08451e08
# ╟─e8fa154d-e60b-407e-84ea-2b247141374d
# ╟─b46666a3-3a24-485d-a1b6-41250fc21613
# ╟─2917cd0f-2869-41f9-8b65-62cb75f3d278
# ╟─ed1ca53b-0b7c-47e6-b4c8-9167731a49fb
# ╟─067fd91d-77c6-436e-880b-45772b27d17a
# ╟─0dfaf005-b39a-44c3-8b51-005feac8d78b
# ╟─2e9567f8-10f6-46d0-b32d-39df2038aaf7
# ╟─444b957f-0bcb-4703-9477-dcfcdc320ca6
# ╟─9512f4de-46ac-4f85-84c9-7a37ab391143
# ╟─2334e293-57f7-4b80-ad50-f889c8f6bab0
# ╠═d10d14c7-2b1a-4dde-8ea4-e5c8dcd86887
# ╟─abc6d492-562f-4186-abf7-70aa825648fe
# ╟─43eaea11-1764-4dee-af79-05aa264937c4
# ╟─47148425-e812-4ec3-ab40-49603b5d53ee
# ╟─285cdbf5-1e73-428a-bc1d-88ea8a731ce0

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function device!(
    x,
    output_ode,
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    ix_range::UnitRange{Int64},
    ode_range::UnitRange{Int64},
    device::DynG,
    sys::PSY.System,
) where {DynG <: PSY.DynamicGenerator}
    #Obtain local device states
    n_states = PSY.get_n_states(device)
    device_states = @view x[ix_range]

    #Obtain references
    sys_Sbase = PSY.get_base_power(sys)
    sys_f0 = PSY.get_frequency(sys)
    sys_ω = get_ω_sys(sys)

    #Update Voltage data
    get_inner_vars(device)[VR_gen_var] = voltage_r[1]
    get_inner_vars(device)[VI_gen_var] = voltage_i[1]

    #Obtain ODEs and Mechanical Power for Turbine Governor
    mdl_tg_ode!(device_states, view(output_ode, ode_range), device)

    #Obtain ODEs for AVR
    mdl_pss_ode!(device_states, view(output_ode, ode_range), sys_ω, device)

    #Obtain ODEs for AVR
    mdl_avr_ode!(device_states, view(output_ode, ode_range), device)

    #Obtain ODEs for Machine
    mdl_machine_ode!(
        device_states,
        view(output_ode, ode_range),
        current_r,
        current_i,
        sys_Sbase,
        sys_f0,
        device,
    )

    #Obtain ODEs for PSY.Shaft
    mdl_shaft_ode!(device_states, view(output_ode, ode_range), sys_f0, sys_ω, device)

    return

end

function device!(
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    device::PSY.Source,
    sys::PSY.System,
)
    mdl_source!(voltage_r, voltage_i, current_r, current_i, device, sys)
    return
end

function device!(
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    device::PSY.PowerLoad,
    sys::PSY.System,
)
    mdl_Zload!(voltage_r, voltage_i, current_r, current_i, device, sys)
    return
end

function device!(
    x,
    output_ode::Vector{T},
    voltage_r,
    voltage_i,
    current_r,
    current_i,
    ix_range::UnitRange{Int64},
    ode_range::UnitRange{Int64},
    device::DynI,
    sys::PSY.System,
) where {DynI <: PSY.DynamicInverter, T <: Real}

    #Obtain local device states
    n_states = PSY.get_n_states(device)
    device_states = @view x[ix_range]

    #Obtain references
    Sbase = PSY.get_base_power(sys)
    sys_f0 = PSY.get_frequency(sys)
    sys_ω = get_ω_sys(sys)

    #Update Voltage data
    get_inner_vars(device)[VR_inv_var] = voltage_r[1]
    get_inner_vars(device)[VI_inv_var] = voltage_i[1]

    #Update V_ref
    V_ref = PSY.get_ext(device)[CONTROL_REFS][V_ref_index]
    get_inner_vars(device)[V_oc_var] = V_ref

    #Obtain ODES for DC side
    mdl_DCside_ode!(device)

    #Obtain ODEs for PLL
    mdl_freq_estimator_ode!(
        device_states,
        view(output_ode, ode_range),
        sys_f0,
        sys_ω,
        device,
    )

    #Obtain ODEs for OuterLoop
    mdl_outer_ode!(device_states, view(output_ode, ode_range), sys_f0, sys_ω, device)

    #Obtain inner controller ODEs and modulation commands
    mdl_inner_ode!(device_states, view(output_ode, ode_range), device)

    #Obtain converter relations
    mdl_converter_ode!(device)

    #Obtain ODEs for output filter
    mdl_filter_ode!(
        device_states,
        view(output_ode, ode_range),
        current_r,
        current_i,
        Sbase,
        sys_f0,
        sys_ω,
        device,
    )

    return
end

ext: jl


text:
import Base: size, ndims, Array, length
using AffineTransforms
using Grid

type NeuroImage{T, P, N} <: AbstractArray{T, N}
    values :: AbstractArray{T, N}
    coordsys :: CoordSystem
    transformations :: CoordmapSet
    function NeuroImage(values :: AbstractArray, coordsys :: CoordSystem, transformations :: CoordmapSet)
        if P > N
            error("Number of physical dimensions $P is too high for $N-dimensional array")
        end
        new(values, coordsys, transformations)
    end
end
NeuroImage{T, N}(values :: AbstractArray{T, N}, phys_ndims :: Int, coordsys :: CoordSystem, transformations :: CoordmapSet) = NeuroImage{T, phys_ndims, N}(values, coordsys, transformations)
NeuroImage(values :: AbstractArray, phys_ndims :: Int, coordsys :: CoordSystem) = NeuroImage(values, phys_ndims, coordsys, CoordmapSet())
NeuroImage(values, phys_ndims :: Int) = NeuroImage(values, phys_ndims, voxel())
NeuroImage(values :: AbstractArray) = NeuroImage(values, min(3, ndims(values)))
Array(img :: NeuroImage) = img.values
coordsys(img :: NeuroImage) = img.coordsys

size(img :: NeuroImage) = size(img.values)
size(img :: NeuroImage, d) = size(img.values, d)
ndims{T, P, N}(img :: NeuroImage{T, P, N}) = N
length(img :: NeuroImage) = length(img.values)
endof(img :: NeuroImage) = endof(img.values)

phys_ndims{T, P, N}(img :: NeuroImage{T, P, N}) = P
phys_size(img :: NeuroImage) = size(img)[1:phys_ndims(img)]

typealias RangeIndex Union(Int, Range{Int}, UnitRange{Int}, Colon)
function getindex{T, P, N}(img :: NeuroImage{T, P, N}, inds :: RangeIndex...)
    new_arr = img.values[inds...]
    if all(dim -> dim > ndims(new_arr) ? true : size(new_arr, dim) == 1, 1:phys_ndims(img))
        return new_arr
    end
    if ndims(new_arr) < P
        new_arr = reshape(new_arr, ntuple(P, ix -> ix <= ndims(new_arr) ? size(new_arr, ix) : 1))
    end
    new_coordsys = voxel()
    new_trans = union(img.transformations, CoordmapSet([Coordmap(new_coordsys, img.coordsys, inds_to_transform(phys_ndims(img), inds[1:phys_ndims(img)]...))]))
    NeuroImage{T, P, ndims(new_arr)}(new_arr, new_coordsys, new_trans)
end

setindex!(img :: NeuroImage, values, inds :: RangeIndex...) = img.values[inds...] = values

find_transformation(from_img :: NeuroImage, target_coordsys :: CoordSystem) = find_transformation(from_img.transformations, from_img.coordsys, target_coordsys, phys_ndims(from_img))
find_transformation(from_coordsys :: CoordSystem, target_img :: NeuroImage) = find_transformation(target_img.transformations, from_coordsys, target_img.coordsys, phys_ndims(target_img))
find_transformation(from_img :: NeuroImage, target_img :: NeuroImage) = find_transformation(union(from_img.transformations, target_img.transformations), from_img.coordsys, target_img._coordsys, phys_ndims(from_img))

location(img :: NeuroImage, coordsys :: CoordSystem, inds :: Integer...) = find_transformation(img, coordsys) * [dim > size(inds, 1) ? 1 : inds[dim] for dim=1:phys_ndims(img)]
location(img :: NeuroImage, inds :: Integer...) = [dim > size(inds, 1) ? 1 : inds[dim] for dim=1:phys_ndims(img)]

function physical_interpolator{T, P, N}(img :: NeuroImage{T, P, N}, interpolator=InterpLinear, boundary_condition=BCnil)
    non_phys_size = ntuple(N - P, dim -> size(img, dim + P))
    arr = reshape(Array(img), ntuple(P + 1, dim -> dim > P ? prod(non_phys_size) : size(img, dim)))
    grid_size = size(arr)
    grid_strides = strides(arr)
    strd = stride(arr, P+1)
    interp_invert!(arr, boundary_condition, interpolator, 1:P)
    ic = InterpGridCoefs(eltype(arr), interpolator, grid_size[1:P], grid_strides[1:P])
    function return_arr(index :: Vector)
        set_position(ic, boundary_condition, false, false, index)
        spec = Array(Float64, non_phys_size)
        index = 1
        for i = 1:prod(non_phys_size)
            spec[i] = interp(ic, arr, index)
            index += strd
        end
        spec
    end
    return_arr
end

function tuple_index{T <: Integer, S, N}(A :: AbstractArray{S, N}, index :: T)
    Astrides = strides(A)
    res = Array(T, N)
    for dim in sortperm(Astrides)[end:-1,1]
        res[dim] = 1 + ifloor((index - 1) / Astrides[dim])
        index = index - res[dim] * Astrides[dim]
    end
    res
end

function ndindex{N}(sz :: NTuple{N, Integer})
    if N == 0
        produce(())
    else
        for subset in @task ndindex(sz[2:end])
            for index in 1:sz[1]
                produce(tuple(index, subset...))
            end
        end
    end
end

function ndindex{T, P, N}(img :: NeuroImage{T, P, N})
    addition = ntuple(N - P, dim -> :)
    for index in @task ndindex{P}(phys_size(image))
        produce(tuple(index..., addition...))
    end
end

function transform_to{T, P, N}(from_img :: NeuroImage{T, P, N}, target_coordsys :: CoordSystem, phys_size :: NTuple{P, Integer}, transformations :: CoordmapSet=CoordmapSet([]), interpolator :: Grid.DataType=InterpLinear)
    sz = tuple(phys_size..., size(from_img)[P+1:N]...)
    new_img = NeuroImage(Array(T, sz), P, target_coordsys, transformations)
    trans = find_transformation(union(from_img.transformations, transformations), target_coordsys, from_img.coordsys, P)
    func = physical_interpolator(from_img, interpolator, BCnan)
    nvalues = prod(sz[P+1:N])
    for index in @task ndindex(sz[1:P])
        new_img[index..., 1:end] = func(trans * [index...])
    end
    new_img
end

transform_to{T, S, P, N, M}(from_img :: NeuroImage{T, P, N}, target_img :: NeuroImage{S, P, M}) = transform_to(from_img, target_img.coordsys, phys_size(target_img), target_img.transformations)

push_transformation!(img :: NeuroImage, transformation, to_coordsys :: CoordSystem, from_coordsys :: CoordSystem) = push!(img.transformations, Coordmap(from_coordsys, to_coordsys, transformation))
push_transformation!(img :: NeuroImage, transformation, to_coordsys :: CoordSystem) = push_transformation!(img, transformation, to_coordsys, img.coordsys)


text:
using Test: @testset, @test
import ConvolutionalRBMs as ConvRBMs
import Zygote
import FiniteDifferences
using FiniteDifferences: central_fdm

@testset "conv n=$n C=$C M=$M" for n in 1:3, C in 1:3, B in 1:3, M=1:3
    kernel = (rand(1:4, n)...,)
    output = (rand(1:4, n)...,)
    inputs = output .+ kernel .- 1
    w = randn(C, kernel..., M)
    v = randn(C, inputs..., B)
    h = randn(M, output..., B)
    I_v2h = zeros(M, output..., B)
    I_h2v = zeros(C, inputs..., B)
    for b in 1:B, k in CartesianIndices(output), μ = 1:M, j in CartesianIndices(kernel), c in 1:C
        i = CartesianIndex(Tuple(j) .+ Tuple(k) .- 1)
        @assert i ∈ CartesianIndices(inputs)
        I_v2h[μ,k,b] += w[c,j,μ] * v[c,i,b]
        I_h2v[c,i,b] += w[c,j,μ] * h[μ,k,b]
    end
    # need rtol here because the (non-compensated) sum loop just above accumulates some error
    @test ConvRBMs.conv_v2h(w, v) ≈ I_v2h rtol=1e-6
    @test ConvRBMs.conv_h2v(w, h) ≈ I_h2v rtol=1e-6

    ∂w, ∂v = Zygote.gradient(w, v) do w, v
        sum(ConvRBMs.conv_v2h(w, v))
    end
    Δw, Δv = FiniteDifferences.grad(central_fdm(5,1), (w,v) -> sum(ConvRBMs.conv_v2h(w,v)), w,v)
    @test ∂w ≈ Δw rtol=1e-6
    @test ∂v ≈ Δv rtol=1e-6

    ∂w, ∂h = Zygote.gradient(w, h) do w, h
        sum(ConvRBMs.conv_h2v(w, h))
    end
    Δw, Δh = FiniteDifferences.grad(central_fdm(5,1), (w,h) -> sum(ConvRBMs.conv_h2v(w,h)), w,h)
    @test ∂w ≈ Δw rtol=1e-6
    @test ∂h ≈ Δh rtol=1e-6
end


text:
module Assembly

using LinearAlgebra
using SparseArrays

using GreenDeltaStar

include("assemPotentialKandV.jl")
export assemPotentialKandV

include("assemPotentialKandVwithDuffy.jl")
export assemPotentialKandVwithDuffy

include("evalPsiFromBEM.jl")
export evalPsiFromBEM

include("rectangleCoil.jl")
export evalPsiFromRectangleCoil
export assemblyRhsRectangleCoil

end


text:
module Bracketing

abstract type BracketingMethod end

struct Bracket
    left::Real
    right::Real

    function Bracket(a, b)
        if a <= b
            new(a, b)
        else
            new(b, a)
        end
    end
end

function search(params::BracketingMethod, f, x_0 = 0)::Real
    bracket = find_bracket(f, x_0)
    search_bracket(params, f, bracket).right
end

function find_bracket(f, x; step = 1e-2, factor = 2.0, max_iter = 1_000_000)::Bracket
    """
    Algorithm 3.1

    `max_iter` is a very high number since upon reaching this limit the algorithm
    will crash with error. This is because a failure to find a bracket means there's
    something wrong with either the algorithm, or the objective function.
    """
    a, ya = x, f(x)
    b, yb = a + step, f(a + step)
    if yb > ya
        a, b = b, a
        ya, yb = yb, ya
        step = step > 0 ? -step : step
    end
    for _ = 1:max_iter
        c, yc = b + step, f(b + step)
        if yc > yb
            return Bracket(a, c)
        end
        a, ya, b, yb = b, yb, c, yc
        step *= factor
    end
    error("Bracketing.find_bracket: max number of iterations reached")
end

struct GoldenSection <: BracketingMethod
    ϵ
    max_iter
    GoldenSection(; ϵ = eps(), max_iter = 100) = new(ϵ, max_iter)
end

function search_bracket(params::GoldenSection, f, bracket::Bracket)::Bracket
    """
    Golden section search
    Algorithm 3.3
    """
    a, b = bracket.left, bracket.right
    max_eval = (b - a) / (params.ϵ * log(Base.MathConstants.golden))
    n = params.max_iter > max_eval ? max_eval : params.max_iter

    ρ = Base.MathConstants.golden - 1
    d = ρ * b + (1 - ρ) * a
    yd = f(d)
    for _ = 1:n - 1
        if abs(a - b) < params.ϵ
            return Bracket(a, b)
        end
        c = ρ * a + (1 - ρ) * b
        yc = f(c)
        if yc < yd
            b, d, yd = d, c, yc
        else
            a, b = b, c
        end
    end
    Bracket(a, b)
end

end # module


text:
using Yeppp
using Base.Test

for T in (Float32, Float64)
    x::Vector{T} = rand(1000)
    y::Vector{T} = rand(1000)
    z::Vector{T} = randn(1000)
    res::Vector{T} = rand(1000)

    @test Yeppp.add!(res, x, y) == x .+ y
    @test Yeppp.multiply!(res, x, y) == x .* y
    @test Yeppp.subtract!(res, x, y) == x .- y
    @test Yeppp.multiply(x, y) == x .* y
    @test Yeppp.subtract(x, y) == x .- y
    @test Yeppp.min!(res, x, y) == min(x, y)
    @test Yeppp.max!(res, x, y) == max(x, y)
    @test Yeppp.min(x, y) == min(x, y)
    @test Yeppp.max(x, y) == max(x, y)

    @test_approx_eq Yeppp.sum(x) sum(x)
    @test_approx_eq Yeppp.sumabs2(x) sumabs2(x)
    @test_approx_eq Yeppp.dot(x, y) dot(x, y)

end

x = rand(1000)
y = rand(1000)
z = randn(1000)
res = rand(1000)

@test_approx_eq Yeppp.exp!(res, x) exp(x)
@test_approx_eq Yeppp.log!(res, x) log(x)
@test_approx_eq Yeppp.sin!(res, x) sin(x)
@test_approx_eq Yeppp.cos!(res, x) cos(x)
@test_approx_eq Yeppp.tan!(res, x) tan(x)

@test_approx_eq Yeppp.sum(x) sum(x)
@test_approx_eq Yeppp.sumabs(x) sumabs(x)
@test_approx_eq Yeppp.sumabs2(x) sumabs2(x)


text:
using MCMCChains
using Distributions
using Random
using Test

Random.seed!(20)

@testset "ess tests" begin
    # Consider the model
    # ```math
    # θ ~ Beta(1, 1)
    # y ~ Bernoulli(θ)
    # ```
    # for variable ``y`` of successes and failures. Given 10 observations of 2 successes
    # and 8 failures, sample four chains with 1000 samples each from the posterior
    # ```math
    # θ ~ Beta(3, 9)
    # ```
    chain = Chains(rand(Beta(3, 9), 1000, 1, 4))

    # Compute the expected sample size
    ess_chain = ess(chain)

    @test all(x -> 1000 ≤ x ≤ 2000, ess_chain[:ess])
    @test all(x -> isapprox(x , 1.0, atol=0.1), ess_chain[:r_hat])
end


text:
using ConstructiveGeometry

r=20.
s1 = sphere(20)
s2 = cube(20)
k=50

u =union(
  [ k, 0, 0]+color("cyan")*(s1 ∪ s2),
	[ k, k, 0]+color("pink")*(s1 ∩ s2),
	[ 0, 0, 0]+color("lime")*(s1 \ s2),
	[ 0, k, 0]+color("gold")*hull(s1, s2)
)
save("sphere_cube.png", u)


text:
function null_error_handler(error_code::Cint,
                  mod::Char,
                  func::Char,
                  eh_data::Ptr{Cvoid})
  nothing
end


text:
using Oceananigans.Grids
using Oceananigans.Grids: x_domain, y_domain
using Oceananigans.Solvers
using Oceananigans.Operators
using Oceananigans.Architectures
using Oceananigans.Fields: ReducedField
using Statistics

import Oceananigans.Solvers: solve!

struct FFTImplicitFreeSurfaceSolver{S, G3, G2, R}
    fft_based_poisson_solver :: S
    three_dimensional_grid :: G3
    horizontal_grid :: G2
    right_hand_side :: R
end

"""
    function FFTImplicitFreeSurfaceSolver(arch, grid, settings)

Return a solver based on the fast Fourier transform for the elliptic equation

```math
[∇² - Az / (g H Δt²)] ηⁿ⁺¹ = (∇ʰ ⋅ Q★ - Az ηⁿ / Δt) / (g H Δt) ,
```

representing an implicit time discretization of the linear free surface evolution equation
for a fluid with constant depth `H`, horizontal areas `Az`, barotropic volume flux `Q★`, time
step `Δt`, gravitational acceleration `g`, and free surface at time-step `n`, `ηⁿ`.
"""
function FFTImplicitFreeSurfaceSolver(arch, grid, gravitational_acceleration::Number, settings)

    grid isa RegRectilinearGrid || grid isa HRegRectilinearGrid ||
        throw(ArgumentError("FFTImplicitFreeSurfaceSolver requires horizontally-regular rectilinear grids."))

    # Construct a "horizontal grid". We support either x or y being Flat, but not both.
    TX, TY, TZ = topology(grid)
    sz = Nx, Ny = (grid.Nx, grid.Ny)
    halo = (grid.Hx, grid.Hy)

    domain = (x = x_domain(grid),
              y = y_domain(grid))

    # Reduce kwargs.
    # Either [1, 2], [1], or [2]
    nonflat_dims = findall(T -> !(T() isa Flat), (TX, TY))

    sz = Tuple(sz[i] for i in nonflat_dims)
    halo = Tuple(halo[i] for i in nonflat_dims)
    domain = NamedTuple((:x, :y)[i] => domain[i] for i in nonflat_dims)

    # Build a "horizontal grid" with a Flat vertical direction.
    # Even if the three dimensional grid is vertically stretched, we can only use
    # FFTImplicitFreeSurfaceSolver with grids that are regularly spaced in the
    # horizontal direction.
    horizontal_grid = RectilinearGrid(arch; topology = (TX, TY, Flat),
                                               size = sz,
                                               halo = halo,
                                               domain...)

    solver = FFTBasedPoissonSolver(arch, horizontal_grid)
    right_hand_side = solver.storage

    return FFTImplicitFreeSurfaceSolver(solver, grid, horizontal_grid, right_hand_side)
end

build_implicit_step_solver(::Val{:FastFourierTransform}, arch, grid, gravitational_acceleration, settings) =
    FFTImplicitFreeSurfaceSolver(arch, grid, gravitational_acceleration, settings)

#####
##### Solve...
#####

function solve!(η, implicit_free_surface_solver::FFTImplicitFreeSurfaceSolver, rhs, g, Δt)
    solver = implicit_free_surface_solver.fft_based_poisson_solver
    grid = implicit_free_surface_solver.three_dimensional_grid
    H = grid.Lz

    m = - 1 / (g * H * Δt^2) # units L⁻²

    # solve! is blocking:
    solve!(η, solver, rhs, m)

    return nothing
end

function compute_implicit_free_surface_right_hand_side!(rhs,
                                                        implicit_solver::FFTImplicitFreeSurfaceSolver,
                                                        g, Δt, ∫ᶻQ, η)

    solver = implicit_solver.fft_based_poisson_solver
    arch = solver.architecture
    grid = implicit_solver.three_dimensional_grid
    H = grid.Lz

    event = launch!(arch, grid, :xy,
                    fft_based_implicit_free_surface_right_hand_side!,
                    rhs, grid, g, H, Δt, ∫ᶻQ, η,
                    dependencies = device_event(arch))

    return event
end

@kernel function fft_based_implicit_free_surface_right_hand_side!(rhs, grid, g, H, Δt, ∫ᶻQ, η)
    i, j = @index(Global, NTuple)
    Az = Azᶜᶜᵃ(i, j, 1, grid)
    δ_Q = flux_div_xyᶜᶜᵃ(i, j, 1, grid, ∫ᶻQ.u, ∫ᶻQ.v)
    @inbounds rhs[i, j, 1] = (δ_Q - Az * η[i, j, 1] / Δt) / (g * H * Δt * Az)
end


text:
#

@doc raw"

Module `_HelpQuEST`

Functions for reading members of QuEST's C structs.

## Exports
* numQubits()
...
"
module _HelpQuEST
export numQubits, isDensityMatrix

import ..QuEST
using ..QuEST: Qreal, Qureg

function numQubits(qr ::Qureg) ::Int
    local n = ccall(:helpQuEST__numQubits, Cint, (Qureg,),
                    qr)
    @assert 0 < n < 200
    return n
end

function isDensityMatrix(qr ::Qureg) ::Bool
    local yes❔ = ccall(:helpQuEST__isDensityMatrix, Cint, (Qureg,),
                      qr)
    @assert yes❔ ∈ 0:1
    return yes❔==1
end

## Julia Module Init #----------------------------------------------------------
#
# Init
#
#-------------------------------------------------------------------------------

using Libdl: dlopen, RTLD_LAZY, RTLD_DEEPBIND, RTLD_GLOBAL

function __init__()
    dlopen("libQuEST",RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)

    @assert Qreal ∈ [ Cdouble, Cfloat ]

    if Qreal == Cdouble
        dlopen(joinpath(@__DIR__,"C","libhelpQuEST_d.so"),
               RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    else Qreal == Cfloat
        dlopen(joinpath(@__DIR__,"C","libhelpQuEST_f.so"),
               RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    end
end

end #^ module _HelpQuEST
#EOF


text:
export MorseBond

"""
    MorseBond(; D, α, r0)

A Morse potential bond between two atoms.
"""
struct MorseBond{T, A, R} <: SpecificInteraction
    D::T
    α::A
    r0::R
end

MorseBond(; D, α, r0) = MorseBond{typeof(D), typeof(α), typeof(r0)}(D, α, r0)

@inline @inbounds function force(b::MorseBond, coord_i, coord_j, box_size)
    dr = vector(coord_i, coord_j, box_size)
    r = norm(dr)
    ralp = exp(-b.α * (r - b.r0))
    c = 2 * b.D * b.α * (1 - ralp) * ralp
    f = c * normalize(dr)
    return SpecificForce2Atoms(f, -f)
end

@inline @inbounds function potential_energy(b::MorseBond, coord_i, coord_j, box_size)
    dr = vector(coord_i, coord_j, box_size)
    r = norm(dr)
    ralp = exp(-b.α * (r - b.r0))
    return b.D * (1 - ralp)^2
end


text:
"""
```
M = imfiltermtx(imgsize, kernel, [border, value])
```

creates a sparse matrix such that matrix-vector multiplication is identical to
calling `Images.imfilter` with the same options:

`imfiltermtx(size(img),kernel,border,value) * vec(img)  ≈  vec(imfilter(img,kernel,border,value))`

See `Images.imfilter` for an explanation of the available options.
Note that boundary condition `"value"` does not support any `value`s other than 0.
"""

function imfiltermtx{T,N}(szimg::NTuple{N,Int}, k::AbstractArray{T,N}, border::AbstractString="replicate", value=zero(T))
  border != "value" || convert(T,value) == zero(T) || throw(ArgumentError("Only value '0' supported for border condition 'value'."))
  npixels = prod(szimg)
  idximg = reshape(1:npixels,szimg)
  idximg = padforfilter(idximg, size(k), border, value)

  cols::Matrix{Int} = im2col(idximg,size(k),"sliding")
  nelems,ncols  = size(cols)
  I::Vector{Int} = vec(repmat(reshape(1:ncols,1,ncols),nelems))
  J::Vector{Int} = vec(cols)
  V::Vector{T}   = repmat(vec(k),ncols)
  if border == "value"
    valid = find(J)
    I,J,V = map(x->x[valid],(I,J,V))
  end

  sparse(I,J,V,ncols,npixels)
end

"""
```
otf = psf2otf(psf, [outsize])
```

creates the optical transfer function (`otf`) from the provided point-spread function (`psf`).
The tuple `outsize` can be used to specify the size of the `otf` array; by default, `outsize = size(psf)`.

The convolution theorem relates convolution with the `psf` in the spatial domain to
multiplication with the `otf` in the frequency domain:

`Images.imfilter(x,Images.reflect(psf),"circular") ≈ real(ifft( fft(x) .* psf2otf(psf,size(x)) ))`

This function is similar to Matlab's `psf2otf`.
See http://www.mathworks.com/help/images/ref/psf2otf.html
"""

function psf2otf{T,N}(psf::AbstractArray{T,N}, outsz::NTuple{N,Int}=size(psf))
  psfsz = size(psf)
  if psfsz != outsz
    pad = map(-,outsz,psfsz)
    all(x -> x >= 0, pad) || throw(DimensionMismatch("psf too large for outsz."))
    psf = padarray(psf,tuple(zeros(Int,N)...),pad,"value",T(0))
  end
  shift = map(x -> -floor(Int,x/2), psfsz)
  fft(circshift(psf,shift))
end

"""
```
imgpad = padforfilter(img, kernelsize, [border, value])
```

returns a padded copy of array `img` such that the following holds true:

`imfilter(padforfilter(img,size(kernel),border,value),kernel,"inner") == imfilter(img,kernel,border,value)`

See `Images.imfilter` for an explanation of options `border` and `value`.
"""

function padforfilter{T,N}(A::AbstractArray{T,N}, sz::NTuple{N,Int}, border::AbstractString="replicate", value=zero(T))
  if border in ("replicate", "circular", "reflect", "symmetric", "value")
    prepad  = ntuple(d->div(sz[d]-1, 2), N)
    postpad = ntuple(d->div(sz[d],   2), N)
    A = padarray(A, prepad, postpad, border, convert(T,value))
  elseif border == "inner" # nothing to do
  else
    throw(ArgumentError("Border condition '$border' not supported."))
  end
  A
end


text:
module ChainRulesCore
using Base.Broadcast: materialize, materialize!, broadcasted, Broadcasted, broadcastable

export AbstractRule, Rule, frule, rrule
export @scalar_rule, @thunk
export extern, cast, store!, Wirtinger, Zero, One, Casted, DNE, Thunk, DNERule

include("differentials.jl")
include("differential_arithmetic.jl")
include("rule_types.jl")
include("rules.jl")
include("rule_definition_tools.jl")
end # module


text:
using Test

@testset "TestEnsemblearn" begin
    include("test_decisiontree.jl")
    include("test_randomforest.jl")
end


text:
using SubspacePreservingSparsification
using SubspacePreservingSparsification: near_zero_row_col
using BenchmarkTools
using MAT
using SparseArrays
using DelimitedFiles
using LinearAlgebra

function init_script()
    ROOT_DIR = realpath(dirname(@__FILE__) * "/..")
    cd(ROOT_DIR)
end

init_script()

sq = "square100"
file = matopen(joinpath("test", "test_data", sq, "untransformed_matrices", "centrosymmetric.mat"))
A = read(file, "A")
close(file)

max_num_bin = 50

# print("Computation\n")
# @benchmark sparsify($A, 0.4, 2, max_num_bin, true)


text:
function read_dots_folds(input_file)
    coords = []
    folds = []
    for line in readlines(input_file)
        if ',' in line
            x,y = map(x->parse(Int,x), split(line, ','))
            push!(coords, [x,y])
        elseif '=' in line
            axis, value = split(split(line)[3], '=')
            value = parse(Int, value)
            push!(folds, [axis, value])
        end
    end
    return coords, folds
end

function fold(input_coords, fold_axis, fold_value)
    output_coords = []
    for coord in input_coords
        if fold_axis == "y" #fold up
            if coord[2] < fold_value
                push!(output_coords, coord)
            else
                folded_coord = [coord[1], coord[2*fold_value - coord[2]]]
                push!(output_coords, folded_coord)
            end
        else # fold left
            if coord[1] < fold_value
                push!(output_coords, coord)
            else
                folded_coord = [2*fold_value - coord[1], coord[2]]
                push!(output_coords, folded_coord)
            end
        end
    end
    # deduplicate by putting into a set
    return Set(output_coords)
end

input_dots, folds = read_dots_folds("input.txt")

output_dots = fold(input_dots, folds[1][1], folds[1][2])

println("number of output dots ", length(output_dots))


text:
module Chips

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export chip

register_normal_element("q__chip", context = @__MODULE__)

function chip(args...; kwargs...)
  q__chip(args...; kw(kwargs)...)
end

end


text:
struct CompressibleEuler <: AbstractTerm end
struct LinearCompressibleEuler <: AbstractTerm end

@inline calc_component!(flux, ::Nothing, _...) = nothing
@inline calc_component!(flux, ::AbstractTerm, _...) = nothing

@inline function calc_component!(flux, ::CompressibleEuler, state, aux, physics)
    ρ = state.ρ
    ρu = state.ρu
    ρe = state.ρe
    ρq = state.ρq
    eos = physics.eos
    parameters = physics.parameters

    p = calc_pressure(eos, state, aux, parameters)
    u = ρu / ρ

    flux.ρ += ρu
    flux.ρu += ρu ⊗ u + p * I
    flux.ρe += (ρe + p) * u
    flux.ρq += ρq * u

    nothing
end

@inline function calc_component!(flux, ::LinearCompressibleEuler, state, aux, physics)
    ρu = state.ρu
    ρᵣ = aux.ref_state.ρ
    pᵣ = aux.ref_state.p
    ρeᵣ = aux.ref_state.ρe

    p = calc_linear_pressure(eos, state, aux, parameters)

    flux.ρ += ρu + p * I
    flux.ρe += (ρeᵣ + pᵣ) * ρu / ρᵣ

    nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
function isHit(phraseWords::Vector{SubString{String}}, words::Vector{SubString{String}})::Bool
    if length(phraseWords) > length(words)
		return false
	end
	m = length(phraseWords)
	n = length(words) - length(phraseWords) + 1
	for idxWord in 1:n
		if words[idxWord] != phraseWords[1]
			continue
        end
		matched = true
		for i in 1:m - 1
			if words[idxWord + i] != phraseWords[1 + i]
				matched = false
				break
            end
		end
		if matched
			return true
        end
    end
	return false
end

function countMAGHits(phrasePairs::Vector{Tuple{String,String}})::Vector{Dict{Int,Int}}
    numPairs = length(phrasePairs)
    pairHits = [Dict{Int,Int}() for i in 1:numPairs]
    pairTokens = [(split(pair[1], "-"), split(pair[2], "-")) for pair in phrasePairs]

    numCPUs = length(Sys.cpu_info())
    jobs = Channel{Vector{Tuple{Int,String}}}(numCPUs)
    jobOutputs = Channel{Vector{Dict{Int,Int}}}(numCPUs)

    function runJob(iCPU::Int)
        myPairHits = [Dict{Int,Int}() for i in 1:numPairs]
        for job in jobs
            for yearTitle in job
                year = yearTitle[1]
                title = yearTitle[2]
                titleTokens = split(title)
                for pairID in 1:numPairs
                    phraseWords1 = pairTokens[pairID][1]
                    phraseWords2 = pairTokens[pairID][2]
                    isHit1 = isHit(phraseWords1, titleTokens)
                    isHit2 = isHit(phraseWords2, titleTokens)
                    if isHit1 && isHit2
                        myPairHits[pairID][year] = get(myPairHits[pairID], year, 0) + 1
                    end
                end
            end
        end
        put!(jobOutputs, myPairHits)
    end

    function makeJobs(batchSize::Int)
        for magID in 0:106
            fileName = "TokenizedMAGTitles/MAGTitles.$magID.txt"
            file = open(fileName)
            job = Vector{Tuple{Int,String}}()
            for line in readlines(file)
                fields = split(line, ". ")
                if length(fields) == 2
                    push!(job, (parse(Int, fields[1]), fields[2]))
                    if length(job) == batchSize
                        put!(jobs, job)
                        job = Vector{Tuple{Int,String}}()
                    end
                end
            end
            if length(job) > 0
                put!(jobs, job)
            end
            close(file)
        end
    end

    bind(jobs, @async makeJobs(1000))
    for iCPU in 1:numCPUs
        Threads.@spawn runJob(iCPU)
    end
    for iCPU in 1:numCPUs
        localPairHits = take!(jobOutputs)
        for pairID in 1:numPairs
            for (year, hits) in localPairHits[pairID]
                pairHits[pairID][year] = get(pairHits[pairID], year, 0) + hits
            end
        end
    end

    pairHits
end

function loadPhrasePairs(fileName::String)::Vector{Tuple{String,String}}
    phrasePairs = Vector{Tuple{String,String}}()
    file = open(fileName)
    for line in readlines(file)
        fields = split(line)
        if length(fields) == 8
            push!(phrasePairs, (fields[3], fields[4]))
            push!(phrasePairs, (fields[5], fields[6]))
            push!(phrasePairs, (fields[7], fields[8]))
        end
    end
    close(file)
    phrasePairs
end

function saveMAGHits(fileName::String, phrasePairs::Vector{Tuple{String,String}}, magHits::Vector{Dict{Int,Int}})
    numPairs = length(phrasePairs)
    file = open(fileName, "w")
    for pairID in 1:numPairs
        phrase1 = phrasePairs[pairID][1]
        phrase2 = phrasePairs[pairID][2]
        print(file, "$phrase1 $phrase2")
        pairYearHits = sort([yearHits for yearHits in magHits[pairID]])
        for yearHits in pairYearHits
            year = yearHits[1]
            hits = yearHits[2]
            print(file, " $year:$hits")
        end
        print(file, "\n")
    end
    close(file)
end

phrasePairs = loadPhrasePairs("DualAspectTopicRepresentation/workspace/NIPS.Top100Topics.txt")
magHits = countMAGHits(phrasePairs)
saveMAGHits("DualAspectTopicRepresentation/workspace/NIPS.Top100Topics.MAGHits.txt", phrasePairs, magHits)

max_stars_repo_path: src/07 evaluate_topics.jl

max_issues_repo_path: src/07 evaluate_topics.jl

max_forks_repo_path: src/07 evaluate_topics.jl


text:
@testset verbose=true "map2salm" begin
    # NOTE: Float16 irfft returns Float32, which leads to type conflicts, so we just don't support
    # map2salm on Float16

    function sYlm(s::Int, ell::Int, m::Int, theta::T, phi::T) where {T<:Real}
        # Eqs. (II.7) and (II.8) of https://arxiv.org/abs/0709.0093v3
        # Note their weird definition w.r.t. `-s`
        k_min = max(0, m + s)
        k_max = min(ell + m, ell + s)
        sin_half_theta, cos_half_theta = sincos(theta / 2)
        return (-1)^(-s) * sqrt((2 * ell + 1) / (4 * T(π))) *
            T(sum(
                (-1) ^ (k)
                * sqrt(factorial(big(ell + m)) * factorial(big(ell - m)) * factorial(big(ell - s)) * factorial(big(ell + s)))
                * (cos_half_theta ^ (2 * ell + m + s - 2 * k))
                * (sin_half_theta ^ (2 * k - s - m))
                / (factorial(big(ell + m - k)) * factorial(big(ell + s - k)) * factorial(big(k)) * factorial(big(k - s - m)))
                for k in k_min:k_max
            )) *
            exp(im * m * phi)
    end

    # Eqs. (II.9) through (II.13) of https://arxiv.org/abs/0709.0093v3
    m2Y22(iota::T, phi) where T = sqrt(5 / (64 * T(π))) * (1 + cos(iota)) ^ 2 * exp(im * 2phi)
    m2Y21(iota::T, phi) where T = sqrt(5 / (16 * T(π))) * sin(iota) * (1 + cos(iota)) * exp(im * phi)
    m2Y20(iota::T, phi) where T = sqrt(15 / (32 * T(π))) * sin(iota) ^ 2
    m2Y2m1(iota::T, phi) where T = sqrt(5 / (16 * T(π))) * sin(iota) * (1 - cos(iota)) * exp(im * -1phi)
    m2Y2m2(iota::T, phi) where T = sqrt(5 / (64 * T(π))) * (1 - cos(iota)) ^ 2 * exp(im * -2phi)

    @testset "Input expressions $T" for T in [BigFloat, Float64, Float32]
        # These are just internal consistency tests of the sYlm function
        # above, against the explicit expressions `mY2.`
        s = -2
        ℓ = 2
        Nϑ = 17
        Nφ = 18
        for (m, m2Y2m) in [(2, m2Y22), (1, m2Y21), (0, m2Y20), (-1, m2Y2m1), (-2, m2Y2m2)]
            f1 = mapslices(ϕθ -> sYlm(s, ℓ, m, ϕθ[2], ϕθ[1]), phi_theta(Nφ, Nϑ, T), dims=[3])
            f2 = mapslices(ϕθ -> m2Y2m(ϕθ[2], ϕθ[1]), phi_theta(Nφ, Nϑ, T), dims=[3])
            @test f1 ≈ f2 atol=10eps(T) rtol=10eps(T)
        end
    end

    @testset "map2salm $T" for T in [BigFloat, Float64, Float32]
        # These test the ability of map2salm to precisely decompose the results of `sYlm`.
        ℓmax = 7
        Nϑ = 2ℓmax + 1
        Nφ = 2ℓmax + 2
        for s in -2:2
            for ℓmin in 0:abs(s)
                for ℓ in abs(s):ℓmax
                    for m in -ℓ:ℓ
                        f = mapslices(ϕθ -> sYlm(s, ℓ, m, ϕθ[2], ϕθ[1]), phi_theta(Nφ, Nϑ, T), dims=[3])
                        computed = map2salm(f, s, ℓmax; ℓmin)
                        expected = zeros(Complex{T}, size(computed))
                        expected[SphericalFunctions.Yindex(ℓ, m, ℓmin)] = one(T)
                        if ≉(computed, expected, atol=30eps(T), rtol=30eps(T))
                            @show T
                            @show ℓmax
                            @show Nϑ
                            @show Nφ
                            @show s
                            @show ℓmin
                            @show ℓ
                            @show m
                            println("computed = $computed")
                            println("expected = $expected")
                            println("max_diff = ", maximum(abs, computed .- expected), ";")
                            println()
                        end
                        @test computed ≈ expected atol=30eps(T) rtol=30eps(T)
                    end
                end
            end
        end
    end

end


text:
# Newton's method

x0 = 1  # The initial guess
f(x) = x^2 - 2  # The function whose root we are trying to find
fprime(x) = 2 * x  # The derivative of the function
tolerance = 10^(-7)  # 7 digit accuracy is desired
epsilon = 10^(-14)  # Do not divide by a number smaller than this
maxIterations = 20  # Do not allow the iterations to continue indefinitely
solutionFound = false  # Have not converged to a solution yet

for i = 1:maxIterations
  y = f(x0)
  yprime = fprime(x0)

  if abs(yprime) < epsilon  # Stop if the denominator is too small
    break
  end

  global x1 = x0 - y/yprime  # Do Newton's computation

  if abs(x1 - x0) <= tolerance  # Stop when the result is within the desired tolerance
    global solutionFound = true
    break
  end

  global x0 = x1  # Update x0 to start the process again
end

if solutionFound
  println("Solution: ", x1)  # x1 is a solution within tolerance and maximum number of iterations
else
  println("Did not converge")  # Newton's method did not converge
end

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Newton Method/NewtonMethod.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
module Utils # Bukdu

"""
    Utils.read_stdout(f)
"""
function read_stdout(f)
    oldout = stdout
    rdout, wrout = redirect_stdout()
    out = @async read(rdout, String)
    f()
    redirect_stdout(oldout)
    close(wrout)
    rstrip(fetch(out))
end

end # module Bukdu.Utils


text:
using Distributions
using LinearAlgebra
using Random
using Test

using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions
using EnsembleKalmanProcesses.Localizers
import EnsembleKalmanProcesses: construct_mean, construct_cov, construct_sigma_ensemble
const EKP = EnsembleKalmanProcesses

@testset "EnsembleKalmanProcess" begin

    # Seed for pseudo-random number generator
    rng_seed = 42
    rng = Random.MersenneTwister(rng_seed)

    ### sanity check on rng:
    d = Parameterized(Normal(0, 1))
    u = ParameterDistribution(Dict("distribution" => d, "constraint" => no_constraint(), "name" => "test"))
    draw_1 = construct_initial_ensemble(u, 1)
    draw_2 = construct_initial_ensemble(u, 1)
    @test !isapprox(draw_1, draw_2)

    # Re-seed rng
    rng = Random.MersenneTwister(rng_seed)

    ### Generate data from a linear model: a regression problem with n_par parameters
    ### and 1 observation of G(u) = A \times u, where A : R^n_par -> R^n_obs
    n_obs = 30                  # dimension of synthetic observation from G(u)
    n_par = 2                  # Number of parameteres
    u_star = [-1.0, 2.0]          # True parameters
    noise_level = 0.1             # Defining the observation noise level (std)
    # Test different AbstractMatrices
    Γy_vec =
        [noise_level^2 * I, noise_level^2 * Matrix(I, n_obs, n_obs), noise_level^2 * Diagonal(Matrix(I, n_obs, n_obs))]
    # Test different localizers
    loc_methods = [RBF(2.0), Delta(), NoLocalization()]

    noise = MvNormal(zeros(n_obs), Γy_vec[1])
    C = [1 -0.9; -0.9 1]          # Correlation structure for linear operator
    A = rand(rng, MvNormal(zeros(2), C), n_obs)'    # Linear operator in R^{n_par x n_obs}

    @test size(A) == (n_obs, n_par)

    #### Define linear model
    function G(u)
        A * u
    end
    # Define nonlinear model
    function G₁(u)
        [sqrt((u[1] - u_star[1])^2 + (u[2] - u_star[2])^2)]
    end

    y_star = G(u_star)
    y_obs = y_star .+ rand(rng, noise)

    @test size(y_obs) == (n_obs,)

    # sum(y-G)^2 ~ n_obs*noise_level^2
    @test isapprox(norm(y_obs .- G(u_star))^2 - n_obs * noise_level^2, 0; atol = 0.05)

    #### Define prior information on parameters
    prior_u1 = Dict("distribution" => Parameterized(Normal(0.0, 0.5)), "constraint" => no_constraint(), "name" => "u1")
    prior_u2 = Dict("distribution" => Parameterized(Normal(3.0, 0.5)), "constraint" => no_constraint(), "name" => "u2")
    prior = ParameterDistribution([prior_u1, prior_u2])

    prior_mean = mean(prior)

    # Assuming independence of u1 and u2
    prior_cov = cov(prior) #convert(Array, Diagonal([sqrt(2.), sqrt(2.)]))

    @testset "EnsembleKalmanSampler" begin
        # Seed for pseudo-random number generator
        rng_seed = 42
        rng = Random.MersenneTwister(rng_seed)

        N_ens = 50 # number of ensemble members (NB for @test throws, make different to N_ens)
        N_iter = 20

        initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ens)
        @test size(initial_ensemble) == (n_par, N_ens)

        # Global scope to compare against EKI
        global eks_final_result = nothing
        global eksobj = nothing
        for Γy in Γy_vec
            eksobj = EKP.EnsembleKalmanProcess(initial_ensemble, y_obs, Γy, Sampler(prior_mean, prior_cov); rng = rng)

            g_ens = G(get_u_final(eksobj))
            @test size(g_ens) == (n_obs, N_ens)
            # as the columns of g are the data, this should throw an error

            # EKS iterations
            for i in 1:N_iter
                params_i = get_u_final(eksobj)
                g_ens = G(params_i)
                if i == 1
                    g_ens_t = permutedims(g_ens, (2, 1))
                    @test_throws DimensionMismatch EKP.update_ensemble!(eksobj, g_ens_t)
                end

                EKP.update_ensemble!(eksobj, g_ens)
            end

            eks_final_result = vec(mean(get_u_final(eksobj), dims = 2))
        end

        # Plot evolution of the EKS particles
        if TEST_PLOT_OUTPUT
            gr()
            p = plot(
                get_u_prior(eksobj)[1, :],
                get_u_prior(eksobj)[2, :],
                seriestype = :scatter,
                label = "Initial ensemble",
            )
            plot!(get_u_final(eksobj)[1, :], get_u_final(eksobj)[2, :], seriestype = :scatter, label = "Final ensemble")
            plot!(
                [u_star[1]],
                xaxis = "u1",
                yaxis = "u2",
                seriestype = "vline",
                linestyle = :dash,
                linecolor = :red,
                label = :none,
            )
            plot!([u_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = :none)
            savefig(p, "EKS_test.png")
        end
    end

    @testset "EnsembleKalmanInversion" begin
        # Seed for pseudo-random number generator
        rng_seed = 42
        rng = Random.MersenneTwister(rng_seed)

        N_ens = 50 # number of ensemble members (NB for @test throws, make different to N_ens)
        N_iter = 20
        iters_with_failure = [5, 8, 9, 15]
        initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ens)

        ekiobj = nothing
        eki_final_result = nothing
        for (Γy, loc_method) in zip(Γy_vec, loc_methods)
            ekiobj = EKP.EnsembleKalmanProcess(
                initial_ensemble,
                y_obs,
                Γy,
                Inversion();
                rng = rng,
                failure_handler_method = SampleSuccGauss(),
                localization_method = loc_method,
            )
            ekiobj_unsafe = EKP.EnsembleKalmanProcess(
                initial_ensemble,
                y_obs,
                Γy,
                Inversion();
                rng = rng,
                failure_handler_method = IgnoreFailures(),
                localization_method = loc_method,
            )

            # some checks
            g_ens = G(get_u_final(ekiobj))
            @test size(g_ens) == (n_obs, N_ens)
            # as the columns of g are the data, this should throw an error
            g_ens_t = permutedims(g_ens, (2, 1))
            @test_throws DimensionMismatch find_ekp_stepsize(ekiobj, g_ens_t)
            Δ = find_ekp_stepsize(ekiobj, g_ens)
            # huge collapse for linear problem so should find timestep should < 1
            if isa(loc_method, NoLocalization)
                @test Δ < 1
            end
            # NOTE We don't use this info, this is just for the test.

            # EKI iterations
            params_i_vec = []
            g_ens_vec = []
            for i in 1:N_iter
                # Check SampleSuccGauss handler
                params_i = get_u_final(ekiobj)
                push!(params_i_vec, params_i)
                g_ens = G(params_i)
                # Add random failures
                if i in iters_with_failure
                    g_ens[:, 1] .= NaN
                end

                EKP.update_ensemble!(ekiobj, g_ens)
                push!(g_ens_vec, g_ens)
                if i == 1
                    g_ens_t = permutedims(g_ens, (2, 1))
                    @test_throws DimensionMismatch EKP.update_ensemble!(ekiobj, g_ens_t)
                end
                @test !any(isnan.(params_i))

                # Check IgnoreFailures handler
                if i <= iters_with_failure[1]
                    params_i_unsafe = get_u_final(ekiobj_unsafe)
                    g_ens_unsafe = G(params_i_unsafe)
                    if i < iters_with_failure[1]
                        EKP.update_ensemble!(ekiobj_unsafe, g_ens_unsafe)
                    elseif i == iters_with_failure[1]
                        g_ens_unsafe[:, 1] .= NaN
                        EKP.update_ensemble!(ekiobj_unsafe, g_ens_unsafe)
                        u_unsafe = get_u_final(ekiobj_unsafe)
                        @test any(isnan.(u_unsafe))
                    end
                end
            end
            push!(params_i_vec, get_u_final(ekiobj))

            @test get_u_prior(ekiobj) == params_i_vec[1]
            @test get_u(ekiobj) == params_i_vec
            @test isequal(get_g(ekiobj), g_ens_vec)
            @test isequal(get_g_final(ekiobj), g_ens_vec[end])
            @test isequal(get_error(ekiobj), ekiobj.err)

            # EKI results: Test if ensemble has collapsed toward the true parameter
            # values
            eki_init_result = vec(mean(get_u_prior(ekiobj), dims = 2))
            eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))
            if isa(loc_method, NoLocalization)
                @test norm(u_star - eki_final_result) < norm(u_star - eki_init_result)
            end
        end

        # Plot evolution of the EKI particles
        #eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))
        if TEST_PLOT_OUTPUT
            gr()
            p = plot(
                get_u_prior(ekiobj)[1, :],
                get_u_prior(ekiobj)[2, :],
                seriestype = :scatter,
                label = "Initial ensemble",
            )
            plot!(get_u_final(ekiobj)[1, :], get_u_final(ekiobj)[2, :], seriestype = :scatter, label = "Final ensemble")
            plot!(
                [u_star[1]],
                xaxis = "u1",
                yaxis = "u2",
                seriestype = "vline",
                linestyle = :dash,
                linecolor = :red,
                label = :none,
            )
            plot!([u_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = :none)
            savefig(p, "EKI_test.png")
        end

        for Γy in Γy_vec
            posterior_cov_inv = (A' * (Γy \ A) + 1 * Matrix(I, n_par, n_par) / prior_cov)
            ols_mean = (A' * (Γy \ A)) \ (A' * (Γy \ y_obs))
            posterior_mean = posterior_cov_inv \ ((A' * (Γy \ A)) * ols_mean + (prior_cov \ prior_mean))

            #### This tests correspond to:
            # NOTE THESE ARE VERY SENSITIVE TO THE CORRESPONDING N_iter
            # EKI provides a solution closer to the ordinary Least Squares estimate
            @test norm(ols_mean - eki_final_result) < norm(ols_mean - eks_final_result)
            # EKS provides a solution closer to the posterior mean
            @test norm(posterior_mean - eks_final_result) < norm(posterior_mean - eki_final_result)

            ##### I expect this test to make sense:
            # In words: the ensemble covariance is still a bit ill-dispersed since the
            # algorithm employed still does not include the correction term for finite-sized
            # ensembles.
            @test abs(sum(diag(posterior_cov_inv \ cov(get_u_final(eksobj), dims = 2))) - n_par) > 1e-5
        end
    end

    @testset "UnscentedKalmanInversion" begin
        # Seed for pseudo-random number generator
        rng_seed = 42
        rng = Random.MersenneTwister(rng_seed)

        N_iter = 20 # number of UKI iterations
        α_reg = 1.0
        update_freq = 0
        process = Unscented(prior_mean, prior_cov; α_reg = α_reg, update_freq = update_freq)
        iters_with_failure = [5, 8, 9, 15]
        failed_particle_index = [1, 2, 3, 1]
        Γy = Γy_vec[3]
        ukiobj = EKP.EnsembleKalmanProcess(y_star, Γy, process; rng = rng, failure_handler_method = SampleSuccGauss())
        ukiobj_unsafe =
            EKP.EnsembleKalmanProcess(y_star, Γy, process; rng = rng, failure_handler_method = IgnoreFailures())
        # UKI iterations
        params_i_vec = []
        g_ens_vec = []
        failed_index = 1
        for i in 1:N_iter
            # Check SampleSuccGauss handler
            params_i = get_u_final(ukiobj)
            push!(params_i_vec, params_i)
            g_ens = G(params_i)
            # Add random failures
            if i in iters_with_failure
                g_ens[:, failed_particle_index[failed_index]] .= NaN
                failed_index += 1
            end

            EKP.update_ensemble!(ukiobj, g_ens)
            push!(g_ens_vec, g_ens)
            if i == 1
                g_ens_t = permutedims(g_ens, (2, 1))
                @test_throws DimensionMismatch EKP.update_ensemble!(ukiobj, g_ens_t)
            end
            @test !any(isnan.(params_i))

            # Check IgnoreFailures handler
            if i <= iters_with_failure[1]
                params_i_unsafe = get_u_final(ukiobj_unsafe)
                g_ens_unsafe = G(params_i_unsafe)
                if i < iters_with_failure[1]
                    EKP.update_ensemble!(ukiobj_unsafe, g_ens_unsafe)
                elseif i == iters_with_failure[1]
                    g_ens_unsafe[:, 1] .= NaN
                    EKP.update_ensemble!(ukiobj_unsafe, g_ens_unsafe)
                    u_unsafe = get_u_final(ukiobj_unsafe)
                    @test any(isnan.(u_unsafe))
                end
            end

        end
        push!(params_i_vec, get_u_final(ukiobj))

        @test get_u_prior(ukiobj) == params_i_vec[1]
        @test get_u(ukiobj) == params_i_vec
        @test isequal(get_g(ukiobj), g_ens_vec)
        @test isequal(get_g_final(ukiobj), g_ens_vec[end])
        @test isequal(get_error(ukiobj), ukiobj.err)

        @test isa(construct_mean(ukiobj, rand(rng, 2 * n_par + 1)), Float64)
        @test isa(construct_mean(ukiobj, rand(rng, 5, 2 * n_par + 1)), Vector{Float64})
        @test isa(construct_cov(ukiobj, rand(rng, 2 * n_par + 1)), Float64)
        @test isa(construct_cov(ukiobj, rand(rng, 5, 2 * n_par + 1)), Matrix{Float64})
        @test isposdef(construct_cov(ukiobj, construct_sigma_ensemble(ukiobj.process, [0.0; 0.0], [1.0 0; 0 0])))

        # UKI results: Test if ensemble has collapsed toward the true parameter
        # values
        uki_init_result = vec(mean(get_u_prior(ukiobj), dims = 2))
        uki_final_result = get_u_mean_final(ukiobj)
        @test norm(u_star - uki_final_result) < norm(u_star - uki_init_result)
        # end

        if TEST_PLOT_OUTPUT
            gr()
            θ_mean_arr = hcat(ukiobj.process.u_mean...)
            N_θ = length(ukiobj.process.u_mean[1])
            θθ_std_arr = zeros(Float64, (N_θ, N_iter + 1))
            for i in 1:(N_iter + 1)
                for j in 1:N_θ
                    θθ_std_arr[j, i] = sqrt(ukiobj.process.uu_cov[i][j, j])
                end
            end

            ites = Array(LinRange(1, N_iter + 1, N_iter + 1))
            p = plot(ites, grid = false, θ_mean_arr[1, :], yerror = 3.0 * θθ_std_arr[1, :], label = "u1")
            plot!(ites, fill(u_star[1], N_iter + 1), linestyle = :dash, linecolor = :grey, label = :none)
            plot!(
                ites,
                grid = false,
                θ_mean_arr[2, :],
                yerror = 3.0 * θθ_std_arr[2, :],
                label = "u2",
                xaxis = "Iterations",
            )
            plot!(ites, fill(u_star[2], N_iter + 1), linestyle = :dash, linecolor = :grey, label = :none)
            savefig(p, "UKI_test.png")
        end

        ### Generate data from G₁(u) with a sparse u
        n_obs = 1                  # Number of synthetic observations from G₁(u)
        n_par = 2                  # Number of parameteres
        u_star = [1.0, 0.0]          # True parameters
        noise_level = 1e-3            # Defining the observation noise level
        Γy = noise_level * Matrix(I, n_obs, n_obs) # Independent noise for synthetic observations
        noise = MvNormal(zeros(n_obs), Γy)

        y_star = G₁(u_star)
        y_obs = y_star + rand(rng, noise)

        @test size(y_star) == (n_obs,)

        #### Define prior information on parameters
        prior_u1 =
            Dict("distribution" => Parameterized(Normal(0.0, 2)), "constraint" => no_constraint(), "name" => "u1")
        prior_u2 =
            Dict("distribution" => Parameterized(Normal(3.0, 2)), "constraint" => no_constraint(), "name" => "u2")
        prior = ParameterDistribution([prior_u1, prior_u2])
    end

    ###
    ###  Calibrate (4): Sparse Ensemble Kalman Inversion
    ###
    @testset "SparseInversion" begin
        # Seed for pseudo-random number generator
        rng_seed = 42
        rng = Random.MersenneTwister(rng_seed)

        n_obs = 1
        n_par = 2
        N_ens = 20 # number of ensemble members
        N_iter = 5 # number of EKI iterations
        iters_with_failure = [1, 3]
        Γy_vec = [noise_level * Matrix(I, n_obs, n_obs), noise_level * I]
        loc_methods = [NoLocalization(), RBF(2.0)]
        # Sparse EKI parameters
        γ = 1.0
        regs = [1e-4, 1e-3]
        uc_idxs = [[1, 2], :]

        initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ens)
        @test size(initial_ensemble) == (n_par, N_ens)

        threshold_values = [0, 1e-2]
        test_names = ["test", "test_thresholded"]

        for (threshold_value, reg, uc_idx, test_name, Γy, loc_method) in
            zip(threshold_values, regs, uc_idxs, test_names, Γy_vec, loc_methods)
            process = SparseInversion(γ, threshold_value, uc_idx, reg)

            ekiobj = EKP.EnsembleKalmanProcess(
                initial_ensemble,
                y_obs,
                Γy,
                process;
                rng = rng,
                failure_handler_method = SampleSuccGauss(),
                localization_method = loc_method,
            )
            ekiobj_unsafe = EKP.EnsembleKalmanProcess(
                initial_ensemble,
                y_obs,
                Γy,
                process;
                rng = rng,
                failure_handler_method = IgnoreFailures(),
            )

            # EKI iterations
            params_i_vec = []
            g_ens_vec = []
            for i in 1:N_iter
                # Check SammpleSuccGauss handler
                params_i = get_u_final(ekiobj)
                push!(params_i_vec, params_i)
                g_ens = hcat([G₁(params_i[:, i]) for i in 1:N_ens]...)
                # Add random failures
                if i in iters_with_failure
                    g_ens[:, 1] .= NaN
                end

                push!(g_ens_vec, g_ens)
                if i == 1
                    g_ens_t = permutedims(g_ens, (2, 1))
                    @test_throws DimensionMismatch EKP.update_ensemble!(ekiobj, g_ens_t)
                    EKP.update_ensemble!(ekiobj, g_ens)
                else
                    EKP.update_ensemble!(ekiobj, g_ens, Δt_new = ekiobj.Δt[1])
                end
                @test !any(isnan.(params_i))
            end
            push!(params_i_vec, get_u_final(ekiobj))

            @test get_u_prior(ekiobj) == params_i_vec[1]
            @test get_u(ekiobj) == params_i_vec
            @test isequal(get_g(ekiobj), g_ens_vec)
            @test isequal(get_g_final(ekiobj), g_ens_vec[end])
            @test isequal(get_error(ekiobj), ekiobj.err)

            # EKI results: Test if ensemble has collapsed toward the true parameter
            # values
            eki_init_result = vec(mean(get_u_prior(ekiobj), dims = 2))
            eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))
            @test norm(u_star - eki_final_result) < norm(u_star - eki_init_result)
            @test sum(eki_final_result .> 0.05) < size(eki_final_result)[1]

            # Plot evolution of the EKI particles
            eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))

            if TEST_PLOT_OUTPUT
                gr()
                p = plot(get_u_prior(ekiobj)[1, :], get_u_prior(ekiobj)[2, :], seriestype = :scatter)
                plot!(get_u_final(ekiobj)[1, :], get_u_final(ekiobj)[2, :], seriestype = :scatter)
                plot!(
                    [u_star[1]],
                    xaxis = "u1",
                    yaxis = "u2",
                    seriestype = "vline",
                    linestyle = :dash,
                    linecolor = :red,
                )
                plot!([u_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red)
                savefig(p, string("SparseEKI_", test_name, ".png"))
            end

            # Test other constructors
            @test isa(SparseInversion(γ), SparseInversion)

        end
    end

    @testset "EnsembleKalmanProcess utils" begin
        g = rand(5, 10)
        g[:, 1] .= NaN
        successful_ens, failed_ens = split_indices_by_success(g)
        @test length(failed_ens) == 1
        @test length(successful_ens) == size(g, 2) - length(failed_ens)
        for i in 2:7
            g[:, i] .= NaN
        end
        @test_logs (:warn,) split_indices_by_success(g)

        u = rand(10, 4)
        @test_logs (:warn,) sample_empirical_gaussian(u, 2)
        @test_throws PosDefException sample_empirical_gaussian(u, 2, inflation = 0.0)
    end

end


text:
include(joinpath("..", "src", "simulation.jl"))

run_resume!()
readdir()
run_batch(reverse(configbatch), batch_name="BAv2")

barabasi_albert(1000, 500)

convert_results()

batchrun = Simulation[]

for file in readdir("results")
    raw = load(joinpath("results", file))
    push!(batchrun, raw[first(keys(raw))])
end

resultcomparison = DataFrame(
    config = String[],
    Agentcount = Int64[],
    Edgecount = Int64[],
    OutdegreeAVG = Float64[],
    OutdegreeSD = Float64[],
    OutdegreeMAX = Int64[],
    IndegreeAVG = Float64[],
    IndegreeSD = Float64[],
    IndegreeMAX = Int64[],
    CentralityAVG = Float64[],
    CentralityMAX = Float64[],
    CentralityMIN = Float64[],
    CCAVG = Float64[],
    CCMAX = Float64[],
    OpinionAVG = Float64[],
    OpinionSD = Float64[]
)

for (index, simulation) in enumerate(batchrun)
    componentsizes = ""
    for i in 1:length(connected_components(simulation.final_state[1]))
        componentsizes = componentsizes * "Component No.$i, $(length(connected_components(simulation.final_state[1])[i])) Agents. "
    end

    if simulation.config.simulation.addfriends == ""
        addfriends = "hybrid"
    elseif simulation.config.simulation.addfriends == "neighborsofneighbors"
        addfriends = "NoN"
    else
        addfriends = simulation.config.simulation.addfriends
    end

    config = (
        string(simulation.config.network.agent_count)
        * "_"
        * addfriends
        * "_"
        * string(simulation.config.opinion_threshs.unfriend)
    )

    push!(
        resultcomparison,
        (
            config,
            nv(simulation.final_state[1]),
            ne(simulation.final_state[1]),
            mean(outdegree(simulation.final_state[1])),
            sd(outdegree(simulation.final_state[1])),
            maximum(outdegree(simulation.final_state[1])),
            mean(indegree(simulation.final_state[1])),
            sd(indegree(simulation.final_state[1])),
            maximum(indegree(simulation.final_state[1])),
            mean(closeness_centrality(simulation.final_state[1])),
            maximum(closeness_centrality(simulation.final_state[1])),
            minimum(closeness_centrality(simulation.final_state[1])),
            mean(local_clustering_coefficient(simulation.final_state[1])),
            maximum(local_clustering_coefficient(simulation.final_state[1])),
            mean([agent.opinion for agent in simulation.final_state[2]]),
            std([agent.opinion for agent in simulation.final_state[2]])
        )
    )

    # println(
    #     "Simulation Run $index \n",
    #     "Configuration: $(simulation.config.network.agent_count) Agents, Friend Mode: $(simulation.config.simulation.addfriends), Unfriend Threshold: $(simulation.config.opinion_threshs.unfriend) \n" ,
    #     "Nodes: $(nv(simulation.final_state[1])), Edges: $(ne(simulation.final_state[1])) \n",
    #     "AVG Outdegree: $(mean(outdegree(simulation.final_state[1]))) \n",
    #     "Components of final Graph: $componentsizes \n"
    # )
    # histogram(outdegree(simulation.final_state[1]))
end

using Statistics
using RCall

@rput resultcomparison
R"save(resultcomparison,file=\"results.Rda\")"
degree
histogram(mean(indegree(batchrun[18].final_state[1])))

filter(row -> row[:TickNr] == 1000, BArun01.agent_log)[:Opinion]

batchrun = Simulation[]
for i in 1:3
    temp = load(joinpath("results", readdir("results")[i]))
    push!(batchrun, temp[first(keys(temp))])
end

using Plots
histogram([agent.opinion for agent in batchrun[3].final_state[2]], bins = 100)

convert_results(specific_run="BAv2_run03.jld2")

run!(Simulation(configbatch[1]), name="BAv3")

readdir("results")

reverse(collect(1:10))
stop_at = 0
if stop_at == 0
    stop_at = 2
end
for i in 10:-1:1
    print(i)
end

using JLD
using JLD2

configbatch = Config[]
reverse(configbatch)
for
    agent_count in [100, 1000],
    unfriend in [0.4, 0.8, 1.2],
    addfriends in ["", "neighborsofneighbors", "random"]

    push!(configbatch, Config(
            network = cfg_net(
                agent_count = agent_count,
                m0 = Int(agent_count/10),
                new_follows = 10
            ),
            simulation = cfg_sim(
                ticks = 1000,
                addfriends = addfriends
            ),
            opinion_threshs = cfg_ot(
                backfire = 0.4,
                befriend = 0.2,
                unfriend = unfriend
            ),
            agent_props = cfg_ag(
                own_opinion_weight = 0.95,
                unfriend_rate = 0.05,
                min_friends_count = 5
            )
        )
    )
end

configbatch


text:
using ModelingToolkit
using Base.Test

# Define some variables
@IVar t
@DVar x(t) y(t) z(t)
@Deriv D'~t # Default of first derivative, Derivative(t,1)
@Param σ ρ β
@Const c=0
@Var a

# Define a differential equation
eqs = [D*x ~ σ*(y-x),
       D*y ~ x*(ρ-z)-y,
       D*z ~ x*y - β*z]
de = DiffEqSystem(eqs,[t],[x,y,z],Variable[],[σ,ρ,β])
ModelingToolkit.generate_ode_function(de)
jac_expr = ModelingToolkit.generate_ode_jacobian(de)
jac = ModelingToolkit.calculate_jacobian(de)
f = DiffEqFunction(de)
W = I - jac
simplify_constants.(inv(W))

# Differential equation with automatic extraction of variables on rhs
de2 = DiffEqSystem(eqs, [t])

function test_vars_extraction(de, de2)
    for el in (:ivs, :dvs, :vs, :ps)
        names2 = sort(collect(var.name for var in getfield(de2,el)))
        names = sort(collect(var.name for var in getfield(de,el)))
        @test names2 == names
    end
end
test_vars_extraction(de, de2)

# Conversion to first-order ODEs #17
@Deriv D3'''~t
@Deriv D2''~t
@DVar u(t) u_tt(t) u_t(t) x_t(t)
eqs = [D3*u ~ 2(D2*u) + D*u + D*x + 1
       D2*x ~ D*x + 2]
de = DiffEqSystem(eqs, [t])
de1 = ode_order_lowering(de)
lowered_eqs = [D*u_tt ~ 2u_tt + u_t + x_t + 1
               D*x_t  ~ x_t + 2
               D*u_t  ~ u_tt
               D*u    ~ u_t
               D*x    ~ x_t]
function test_eqs(eqs1, eqs2)
    eq = true
    for i in eachindex(eqs1)
        lhs1, lhs2 = eqs1[i].args[1], eqs2[i].args[1]
        for f in fieldnames(typeof(lhs1))
            eq = eq && isequal(getfield(lhs1, f), getfield(lhs2, f))
        end
        eq = eq && isequal(eqs1[i].args[2], eqs2[i].args[2])
    end
    @test_broken eq
end
test_eqs(de1.eqs, lowered_eqs)

# Internal calculations
eqs = [a ~ y-x,
       D*x ~ σ*a,
       D*y ~ x*(ρ-z)-y,
       D*z ~ x*y - β*z]
de = DiffEqSystem(eqs,[t],[x,y,z],[a],[σ,ρ,β])
ModelingToolkit.generate_ode_function(de)
jac = ModelingToolkit.calculate_jacobian(de)
f = DiffEqFunction(de)

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs,[x,y,z],[σ,ρ,β])
ns2 = NonlinearSystem(eqs)
for el in (:vs, :ps)
    names2 = sort(collect(var.name for var in getfield(ns2,el)))
    names = sort(collect(var.name for var in getfield(ns,el)))
    @test names2 == names
end

ModelingToolkit.generate_nlsys_function(ns)

# Now nonlinear system with only variables
@Var x y z
@Param σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs)
nlsys_func = ModelingToolkit.generate_nlsys_function(ns)
jac = ModelingToolkit.generate_nlsys_jacobian(ns)
f = @eval eval(nlsys_func)

# Intermediate calculations
# Define a nonlinear system
eqs = [a ~ y-x,
       0 ~ σ*a,
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs,[x,y,z],[σ,ρ,β])
nlsys_func = ModelingToolkit.generate_nlsys_function(ns)
jac = ModelingToolkit.calculate_jacobian(ns)
jac = ModelingToolkit.generate_nlsys_jacobian(ns)


text:
using Espresso
using Test

include("utils_test.jl")
include("rewrite_test.jl")
include("simplify_test.jl")
include("exnode_test.jl")
include("exgraph_test.jl")
# include("tracked_test.jl")  -- tracked arrays are discontinued, use Yota.jl instead


text:
@symbol_func function P_L(cur_reactor::AbstractReactor)
  cur_P = P_alpha(cur_reactor)

  cur_P += P_H(cur_reactor)

  cur_P += P_ohmic(cur_reactor)

  cur_P -= P_BR(cur_reactor)

  cur_P
end


text:
# an AbstractIndex is a thing that can be used to look up ordered things by name, but that
# will also accept a position or set of positions or range or other things and pass them
# through cleanly.
# an Index is the usual implementation.
# a SimpleIndex only works if the things are integer indexes, which is weird.
typealias Indices Union(Real, AbstractVector{Real})

abstract AbstractIndex

type Index <: AbstractIndex   # an OrderedDict would be nice here...
    lookup::Dict{ByteString, Indices}      # name => names array position
    names::Vector{ByteString}
end
function Index{T <: ByteString}(x::Vector{T})
    x = make_unique(convert(Vector{ByteString}, x))
    Index(Dict{ByteString, Indices}(tuple(x...), tuple([1:length(x)]...)), x)
end
Index() = Index(Dict{ByteString, Indices}(), ByteString[])
Base.length(x::Index) = length(x.names)
Base.names(x::Index) = copy(x.names)
Base.copy(x::Index) = Index(copy(x.lookup), copy(x.names))
Base.deepcopy(x::Index) = Index(deepcopy(x.lookup), deepcopy(x.names))
Base.isequal(x::Index, y::Index) = isequal(x.lookup, y.lookup) && isequal(x.names, y.names)

# I think this should be Vector{T <: String}
function names!(x::Index, nm::Vector)
    if length(nm) != length(x)
        error("lengths don't match.")
    end
    for i in 1:length(nm)
        delete!(x.lookup, x.names[i])
        x.lookup[nm[i]] = i
    end
    x.names = nm
    return x
end

function rename!(x::Index, nms)
    for (from, to) in nms
        if haskey(x, from)
            if haskey(x, to)
                error("Tried renaming $from to $to, when $to already exists in the Index.")
            end
            x.lookup[to] = col = pop!(x.lookup, from)
            if !isa(col, Array)
                x.names[col] = to
            end
        end
    end
    return x
end

rename!(x::Index, from, to) = rename!(x, zip(from, to))
rename!(x::Index, from::String, to::String) = rename!(x, ((from, to),))
rename!(x::Index, f::Function) = rename!(x, [(x,f(x)) for x in x.names])

rename(x::Index, args...) = rename!(copy(x), args...)

Base.haskey(x::Index, key::String) = haskey(x.lookup, key)
Base.haskey(x::Index, key::Symbol) = haskey(x.lookup, string(key))
Base.haskey(x::Index, key::Real) = 1 <= key <= length(x.names)
Base.keys(x::Index) = names(x)

function Base.push!(x::Index, nm::String)
    x.lookup[nm] = length(x) + 1
    push!(x.names, nm)
    return x
end

function Base.delete!(x::Index, idx::Integer)
    # reset the lookup's beyond the deleted item
    for i in (idx + 1):length(x.names)
        x.lookup[x.names[i]] = i - 1
    end
    gr = get_groups(x)
    delete!(x.lookup, x.names[idx])
    splice!(x.names, idx)
    # fix groups:
    for (k,v) in gr
        newv = [[haskey(x, vv) ? vv : ASCIIString[] for vv in v]...]
        set_group(x, k, newv)
    end
    return x
end

function Base.delete!(x::Index, nm::String)
    if !haskey(x.lookup, nm)
        return x
    end
    idx = x.lookup[nm]
    return delete!(x, idx)
end

Base.getindex(x::Index, idx::String) = x.lookup[idx]
Base.getindex(x::Index, idx::Symbol) = x.lookup[string(idx)]
Base.getindex(x::AbstractIndex, idx::Real) = int(idx)
Base.getindex(x::AbstractIndex, idx::AbstractDataVector{Bool}) = getindex(x, array(idx, false))
Base.getindex{T}(x::AbstractIndex, idx::AbstractDataVector{T}) = getindex(x, dropna(idx))
Base.getindex(x::AbstractIndex, idx::AbstractVector{Bool}) = find(idx)
Base.getindex(x::AbstractIndex, idx::Ranges) = [idx]
Base.getindex{T <: Real}(x::AbstractIndex, idx::AbstractVector{T}) = convert(Vector{Int}, idx)
Base.getindex{T <: String}(x::AbstractIndex, idx::AbstractVector{T}) = [[x.lookup[i] for i in idx]...]
Base.getindex{T <: Symbol}(x::AbstractIndex, idx::AbstractVector{T}) = [[x.lookup[string(i)] for i in idx]...]

type SimpleIndex <: AbstractIndex
    length::Integer
end
SimpleIndex() = SimpleIndex(0)
Base.length(x::SimpleIndex) = x.length
Base.names(x::SimpleIndex) = nothing

# Chris's idea of namespaces adapted by Harlan for column groups
function set_group(idx::Index, newgroup, names)
    if !haskey(idx, newgroup) || isa(idx.lookup[newgroup], Array)
        idx.lookup[newgroup] = [[idx.lookup[nm] for nm in names]...]
    end
end
function set_groups(idx::Index, gr::Dict{ByteString,Vector{ByteString}})
    for (k,v) in gr
        if !haskey(idx, k)
            idx.lookup[k] = [[idx.lookup[nm] for nm in v]...]
        end
    end
end
function get_groups(idx::Index)
    gr = Dict{ByteString,Vector{ByteString}}()
    for (k,v) in idx.lookup
        if isa(v,Array)
            gr[k] = idx.names[v]
        end
    end
    gr
end
function is_group(idx::Index, name::ByteString)
  if haskey(idx, name)
    return isa(idx.lookup[name], Array)
  else
    return false
  end
end

# special pretty-printer for groups, which are just Dicts.
function pretty_show(io::IO, gr::Dict{ByteString,Vector{ByteString}})
    allkeys = keys(gr)
    for k = allkeys
        print(io, "$(k): ")
        print(io, join(gr[k], ", "))
        if k != last(allkeys)
            print(io, "; ")
        end
    end
end


text:
const LSliced = LArray{T,N,D,Syms} where {T,N,D<:AbstractArray{T,N},Syms<:Tuple}
const SLSliced = SLArray{S,T,N,L,Syms} where {S,T,N,L,Syms<:Tuple}
const Sliced = Union{SLSliced,LSliced}

# Allow chained getproperty with partial application of the symbol lookup
struct SlicePartial{Syms1,P}
    __p::P
end

@inline function Base.getproperty(p::SlicePartial{Val{s1}},s2::Symbol) where s1
  getfield(p, :__p)[Val(s1), Val(s2)]
end

@inline function Base.setproperty!(p::SlicePartial{Val{s1}},s2::Symbol,y) where s1
  getfield(p, :__p)[Val(s1), Val(s2)] = y
end

@inline function Base.getproperty(x::Sliced,s::Symbol)
  if s == :__x
    return getfield(x,:__x)
  end
  SlicePartial{Val{s},typeof(x)}(x)
end

#####################################
# Array Interface
# (getindex/getproperty methods are shared with SLArray
#####################################
Base.propertynames(::LArray{T,N,D,S}) where {T,N,D,S<:Tuple{Syms1,Syms2}} where {Syms1,Syms2} = Syms1, Syms2
symnames(::Type{LArray{T,N,D,S}}) where {T,N,D,S<:Tuple{Syms1,Syms2}} where {Syms1,Syms2} = Syms1, Syms2

@inline Base.getindex(x::Sliced,i::Int...) = getfield(x,:__x)[i...]
@inline Base.getindex(x::Sliced,s1::Symbol,s2::Symbol) = getindex(x,Val(s1),Val(s2))
@inline Base.getindex(x::Sliced,i1,s2::Symbol) = getindex(x,i1,Val(s2))
@inline Base.getindex(x::Sliced,s1::Symbol,i2) = getindex(x,Val(s1),i2)

@inline @generated function Base.getindex(x::Sliced,::Val{s1},::Val{s2}) where {s1,s2}
  idx1 = findfirst(y->y==s1,symnames(x)[1])
  idx2 = findfirst(y->y==s2,symnames(x)[2])
  :(getfield(x,:__x)[$idx1, $idx2])
end

@inline @generated function Base.getindex(x::Sliced,i1::Number,::Val{s2}) where s2
  idx2 = findfirst(y->y==s2,symnames(x)[2])
  :(getfield(x,:__x)[i1, $idx2])
end

@inline @generated function Base.getindex(x::Sliced,::Val{s1},i2::Number) where s1
  idx1 = findfirst(y->y==s1,symnames(x)[1])
  :(getfield(x,:__x)[$idx1, i2])
end

@inline Base.setindex!(x::LSliced,y,i...) = getfield(x,:__x)[i...] = y
@inline Base.setindex!(x::LSliced,y,s1::Symbol,s2::Symbol) = setindex!(x,y,Val(s1),Val(s2))
@inline Base.setindex!(x::LSliced,y,s1::Symbol,i2) = setindex!(x,y,Val(s1),i2)
@inline Base.setindex!(x::LSliced,y,i1,s2::Symbol) = setindex!(x,y,i1,Val(s2))

@inline @generated function Base.setindex!(x::LSliced,y,::Val{s1},::Val{s2}) where {s1, s2}
    idx1 = findfirst(y->y==s1,symnames(x)[1])
    idx2 = findfirst(y->y==s2,symnames(x)[2])
    :(x.__x[$idx1, $idx2] = y)
end

@inline @generated function Base.setindex!(x::LSliced,y,::Val{s1},i2::Number) where s1
    idx1 = findfirst(y->y==s1,symnames(x)[1])
    :(x.__x[$idx1, i2] = y)
end

@inline @generated function Base.setindex!(x::LSliced,y,i1::Number,::Val{s2}) where s2
    idx2 = findfirst(y->y==s2,symnames(x)[2])
    :(x.__x[i1, $idx2] = y)
end

#####################################
# Broadcast
#####################################
LAStyle{T,N,Syms}(x::Val{i}) where {T,N,Syms<:Tuple,i} = LAStyle{T,N,Syms}()

"""
    @LSliced Eltype Size Names
    @LSliced Values Names

Creates a `LArray` where the rows and columns are labelled instead of individual cells.
Names are determined from the `Names` tuples and values determined from the `Values` array.
Otherwise, and eltype and size are used to make an LArray with undefined values.

For example:

    a = @LSliced Float64 (4,2) (:a,:b,:c,:d), (:x,:y,:z)
    b = @LSliced [1 2; 3 4; 5 6] (:a,:b,:c), (:x,:y)
"""
macro LSliced(vals,syms)
  vals = esc(vals)
  syms = esc(syms)
  return quote
      LArray{Tuple{$syms...,}}($vals)
  end
end

macro LSliced(type,size,syms)
  type = esc(type)
  size = esc(size)
  syms = esc(syms)
  return quote
      LArray{Tuple{$syms...,}}(Array{$type}(undef,$size...))
  end
end

# """
# dimSymbols(dimSymbols(::LArray{T,N,D,Syms}, dim::Int) where {Syms<:Tuple}

# Returns the labels of the `LArray` associated with dimension `dim`.

# For example:

#     A = @LSliced [1 2; 3 4; 5 6] (:a,:b,:c), (:x, :y)
#     dimSymbols(A,1)  # (:a, :b, :c)
# """
# dimSymbols(::LArray{T,N,D,Syms}, dim::Int) where
# {T,N,D<:AbstractArray{T,N},Syms<:Tuple} = Syms.parameters[dim]

# "returns dimSymbols(,1)"
# rowSymbols(::LArray{T,N,D,Syms}) where
# {T,N,D<:AbstractArray{T,N},Syms<:Tuple} = Syms.parameters[1]

# "returns dimSymbols(,2)"
# colSymbols(::LArray{T,N,D,Syms}) where
# {T,N,D<:AbstractArray{T,N},Syms<:Tuple} = Syms.parameters[2]

symbols(::LArray{T,N,D,Syms}) where
{rows,cols,T,N,D,Syms<:Tuple{rows,cols}} = rows,cols


text:
@testset "math partitioning" begin
    s = raw"A $B \com{C}$ D"
    parts = FP.md_partition(s)
    mathb = parts[2]
    @test mathb.name == :MATH_INLINE
    parts = FP.math_partition(FP.content(mathb))
    @test parts[1].name == :TEXT
    @test parts[2].name == :LX_COMMAND
    @test parts[3].name == :CU_BRACKETS
    @test FP.content(parts[1]) == "B "
    @test FP.content(parts[3]) == "C"
end


text:
struct SharedQSpec{T} <: SpatQSimSpec
    realization::Int
    share_scale::T
    prep_file::String

    function SharedQSpec(realization::Int,
                         share_scale::T,
                         prep_file::String = "prep.h5") where T
        new{T}(realization, share_scale, prep_file)
    end
end

# share_scales() = 0:0.2:1
# share_scale_idx(spec::SharedQSpec) = findfirst(isapprox(spec.share_scale),
#                                                share_scales())

sim_value(spec::SharedQSpec) = spec.share_scale
# survey_qdev_scale(spec::SharedQSpec) = spec.share_scale
function survey_catchability(spec::SharedQSpec, prep::SpatQSimPrep)
    qdevs = transform_log_qdevs(log_qdevs(prep),
                                sim_value(spec) * comm_qdev_scale(spec))
    Catchability(base_catchability().catchability .* qdevs)
end

# function simstudy_dir(spec::SharedQSpec; base_dir = ".")
#     joinpath(base_dir, "sharedq")
# end

# """
#     file_paths(spec::SharedQSpec)

# Returns paths for population state HDF5 array, population total CSV, and catch
# CSV. Will create directories as needed. For example, for the 5th realization
# using the second catchability deviation scaling factor, this function will
# return:

#     repl_05/sharedq_02_popstate.h5
#     repl_05/sharedq_02_popstate.csv
#     repl_05/sharedq_02_catch.csv
# """
# function file_paths(spec::SharedQSpec; base_dir = ".")
#     rlz = realization(spec)

#     base_dir = joinpath(simstudy_dir(spec; base_dir = base_dir),
#                          "repl_" * string(rlz, pad = 2))
#     if !isdir(base_dir)
#         mkpath(base_dir)
#     end

#     file_base = joinpath(base_dir,
#                          "sharedq_" * string(share_scale_idx(spec), pad = 2))
#     file_base .* ["_popstate.h5",
#                   "_popstate.csv",
#                   "_catch.csv"]
# end

# """
#     run_sharedq_sim(n_repl::Int = 100, prep_file = "prep.h5")

# Simulate `n_repl` replicates of a fishery with different magnitudes of spatial
# variation in catchability for the commercial fleet (defined in the function
# `qdev_scales`).
# """
# function run_sharedq_sim(n_repl::Int = 100, prep_file = "prep.h5")
#     share_scale_vec = share_scales()
#     for rlz in 1:n_repl
#         for sc in share_scale_vec
#             spec = SharedQSpec(rlz, sc, prep_file)
#             setup = SpatQSimSetup(spec)
#             result = simulate(setup)
#             save(result)
#         end
#     end
#     nothing
# end


text:
using GeometricAlgebra
using DataStructures
using MacroTools: postwalk, prewalk, @capture
import BenchmarkTools
import Grassmann

include("printing.jl")
include("utils.jl")

function run_benchmark(name, expr::Expr, level)
    indent_level = level * 2
    result = @eval @localbenchmark $expr
    replprint(string(minimum(result)), newline=1, prefix=name * ": "; indent_level)
end

function run_benchmark(name, group, level)
    indent_level = level * 2
    replprint(name; indent_level, newline=1, bold=true, color=color_levels[level])
    for (k, v) ∈ group
        run_benchmark(k, v, level + 1)
    end
end

function run_benchmarks(suite)
    replprint("Benchmarking geometric algebras", newline=1, color=:yellow)

    ga_suite = prepare_suite!(deepcopy(suite))
    grassmann_suite = prepare_suite!(deepcopy(suite); grassmann=true)

    run_benchmark("GeometricAlgebra", ga_suite, 1)
    run_benchmark("Grassmann", grassmann_suite, 1)
end

function prepare_suite!(suite; grassmann=false)
    for (k, v) ∈ suite
        if v isa Expr
            if grassmann
                suite[k] = postwalk(suite[k]) do x
                    if is_blade_symbol(x)
                        Symbol(replace(string(x), r"^v" => "g"))
                    else
                        x
                    end
                end
            end
        else
            prepare_suite!(v; grassmann)
        end
    end
    suite
end

is_blade_symbol(x) = x isa Symbol && startswith(string(x), r"v\d+")

make_suite(base) = make_suite!(OrderedDict(), base)

function make_suite!(suite, base)
    res = DefaultOrderedDict(() -> OrderedDict())
    for (k, v) ∈ base
        if v isa Vector{Expr}
            res[k] = OrderedDict(string.(v) .=> v)
        else
            res[k] = make_suite(v)
        end
    end
    res
end

color_levels = [:red, :cyan]

@basis "+++"
Grassmann.@basis "+++" G g

suite_base = DefaultOrderedDict(() -> OrderedDict())
suite_base["Geometric product"] = [
    :(5v1 * 5v2),
]
suite_base["Addition"] = [
    :(5v1 + 5v2),
]
suite_base["Mixed"]= [
    :((5v1 + 3v3 + 1v12) * 5v2),
]

suite = make_suite(suite_base)

run_benchmarks(suite)


text:
"""
    Hopper
    	model from "Dynamically Stable Legged Locomotion"
		x = (px, py, t, r)
"""
struct Hopper{T}
    n::Int
    m::Int
    d::Int

    mb::T # mass of body
    ml::T # mass of leg
    Jb::T # inertia of body
    Jl::T # inertia of leg

    μ::T  # coefficient of friction
    g::T  # gravity

    qL::Vector{T}
    qU::Vector{T}

    nq
    nu
    nc
    nf
    nb
	ns

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s
end

# Dimensions
nq = 4 # configuration dimension
nu = 2 # control dimension
nc = 1 # number of contact points
nf = 2 # number of faces for friction cone
nb = nc * nf
ns = 1

# Parameters
g = 9.81 # gravity
μ = 1.0  # coefficient of friction
mb = 10.0 # body mass
ml = 1.0  # leg mass
Jb = 2.5 # body inertia
Jl = 0.25 # leg inertia

n = 2 * nq
m = nu + nc + nb + nc + nb + ns
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

# Kinematics
kinematics(::Hopper, q) = [q[1] + q[4] * sin(q[3]), q[2] - q[4] * cos(q[3])]

# Methods
M_func(model::Hopper, q) = Diagonal(@SVector [
								 model.mb + model.ml,
								 model.mb + model.ml,
								 model.Jb + model.Jl,
								 model.ml])

G_func(model::Hopper, q) = @SVector [0.0,
									 (model.mb + model.ml) * model.g,
									 0.0,
									 0.0]

function ϕ_func(::Hopper, q)
    @SVector [q[2] - q[4] * cos(q[3])]
end

N_func(::Hopper, q) = @SMatrix [0.0 1.0 (q[4] * sin(q[3])) (-1.0 * cos(q[3]))]

function P_func(::Hopper, q)
    @SMatrix [1.0 0.0 (q[4] * cos(q[3])) sin(q[3]);
              -1.0 0.0 (-1.0 * q[4] * cos(q[3])) -1.0 * sin(q[3])]
end

B_func(::Hopper, q) = @SMatrix [0.0 0.0 1.0 0.0;
                                -sin(q[3]) cos(q[3]) 0.0 1.0]

function fd(model::Hopper, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

    [q2⁺ - q2⁻;
	((1.0 / h) * (M_func(model, q1) * (SVector{4}(q2⁺) - SVector{4}(q1))
    - M_func(model, q2⁺) * (SVector{4}(q3) - SVector{4}(q2⁺)))
    + transpose(B_func(model, q3)) * SVector{2}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{1}(λ)
    + transpose(P_func(model, q3)) * SVector{2}(b)
    - h * G_func(model, q2⁺))]
end

function maximum_dissipation(model::Hopper, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]
	ψ = u[model.idx_ψ]
	ψ_stack = ψ[1] * ones(model.nb)
	η = u[model.idx_η]
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function friction_cone(model::Hopper, u)
	λ = u[model.idx_λ]
	b = u[model.idx_b]
	return @SVector [model.μ * λ[1] - sum(b)]
end

r = 0.7
qL = -Inf * ones(nq)
qU = Inf * ones(nq)
qL[4] = r / 2.0
qU[4] = r

model = Hopper(n, m, d,
			   mb, ml, Jb, Jl,
			   μ, g,
			   qL, qU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

# Visualization
function visualize!(vis, model::Hopper, q; Δt = 0.1)
    r_foot = 0.05
    r_leg = 0.5 * r_foot

    setobject!(vis["body"], Sphere(Point3f0(0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = RGBA(0, 1, 0, 1.0)))

    setobject!(vis["foot"], Sphere(Point3f0(0),
        convert(Float32, r_foot)),
        MeshPhongMaterial(color = RGBA(1, 0, 0, 1.0)))

    n_leg = 100
    for i = 1:n_leg
        setobject!(vis["leg$i"], Sphere(Point3f0(0),
            convert(Float32, r_leg)),
            MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
    end

    p_leg = [zeros(3) for i = 1:n_leg]
    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        p_body = [q[t][1], 0.0, q[t][2]]
        p_foot = [kinematics(model, q[t])[1], 0.0, kinematics(model, q[t])[2]]

        q_tmp = Array(copy(q[t]))
        r_range = range(0, stop = q[t][4], length = n_leg)
        for i = 1:n_leg
            q_tmp[4] = r_range[i]
            p_leg[i] = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]
        end
        q_tmp[4] = q[t][4]
        p_foot = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]

        z_shift = [0.0; 0.0; r_foot]

        MeshCat.atframe(anim, t) do
            settransform!(vis["body"], Translation(p_body + z_shift))
            settransform!(vis["foot"], Translation(p_foot + z_shift))

            for i = 1:n_leg
                settransform!(vis["leg$i"], Translation(p_leg[i] + z_shift))
            end
        end
    end

    MeshCat.setanimation!(vis, anim)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using SolidStateDetectors
using Plots
# plotlyjs(); # backend for interactive plots (Jupyter)

##
det_config_json = SSD_examples[:BEGe]
simulation = Simulation{Float32}(det_config_json)

##
plot_det = plot(simulation.detector)
# png(plot_det, "plots_md/ssd_geometry.png")

##
using LegendGeSim

##
# import Pkg
# Pkg.add(url="https://github.com/legend-exp/LegendTestData.jl.git")
# Pkg.build("LegendTestData")

##
using LegendTestData
testdata_path = joinpath(LegendTestData.legend_test_data_path(), "data", "ldsim")
# path to my detector metadata
det_metadata = joinpath(testdata_path, "invcoax-metadata.json")

##
# convert LEGEND metadata to SSD configuration (only geometry!)
det_config_ssd = LegendGeSim.ssd_config(det_metadata)
# now we can plug it into SSD
simulation1 = Simulation(SolidStateDetector{Float32}(det_config_ssd))
# it's a bit boring cause we see the same detector as before
plot_leg = plot(simulation1.detector)
# png(plot_leg, "plots_md/metadata_geometry.png")

##
detector = LegendGeSim.simulate_detector(det_metadata, "configs/detector_study_ssd.json")

##
calculate_capacitance(detector)
get_active_volume(detector.point_types)

##
plot(detector.electric_potential)

##
plot_ef = plot(detector.electric_field, φ = 0.0)
plot_electric_fieldlines!(detector, φ = 0.0)
# png(plot_ef, "plots_md/metadata_ef.png")

##
detector = LegendGeSim.simulate_detector(det_metadata, "configs/detector_study_fieldgen.json")