text:
@testset "JLD2" begin

    function test_model_data(model, other)
        @test model.scheduler == other.scheduler
        @test model.rng == other.rng
        @test model.maxid.x == other.maxid.x
    end

    function test_space(space::GridSpace, other)
        @test size(space.s) == size(other.s)
        @test all(length(other.s[pos]) == length(space.s[pos]) for pos in eachindex(space.s))
        @test all(all(x in other.s[pos] for x in space.s[pos]) for pos in eachindex(space.s))
        @test space.metric == other.metric
    end

    function test_space(space::ContinuousSpace, other)
        test_space(space.grid, other.grid)
        @test space.update_vel! == other.update_vel!
        @test space.dims == other.dims
        @test space.spacing == other.spacing
        @test space.extent == other.extent
    end

    test_costmetric(metric, other) = @test false
    test_costmetric(
        metric::Pathfinding.DirectDistance{D},
        other::Pathfinding.DirectDistance{D}
    ) where {D} = @test metric.direction_costs == other.direction_costs

    test_costmetric(
        metric::Pathfinding.MaxDistance{D},
        other::Pathfinding.MaxDistance{D}
    ) where {D} = @test true

    function test_costmetric(
        metric::Pathfinding.PenaltyMap{D},
        other::Pathfinding.PenaltyMap{D}
    ) where {D}
        @test metric.pmap == other.pmap
        test_costmetric(metric.base_metric, other.base_metric)
    end

    function test_astar(astar, other)
        @test typeof(astar) == typeof(other)
        @test astar.agent_paths == other.agent_paths
        @test astar.dims == other.dims
        @test astar.neighborhood == other.neighborhood
        @test astar.admissibility == other.admissibility
        @test astar.walkmap == other.walkmap
        test_costmetric(astar.cost_metric, other.cost_metric)
    end

    function test_pathfinding_model(model, other)
        # agent data
        @test nagents(other) == nagents(model)
        @test all(haskey(other.agents, i) for i in allids(model))
        @test all(model[i].pos == other[i].pos for i in allids(model))
        # model data
        test_model_data(model, other)
        # space data
        @test typeof(model.space) == typeof(other.space)    # to check periodicity
        test_space(model.space, other.space)
        # pathfinder data
        test_astar(model.pathfinder, other.pathfinder)
    end

    @testset "No space" begin
        model, _ = Models.hk()
        AgentsIO.save_checkpoint("test.jld2", model)
        other = AgentsIO.load_checkpoint("test.jld2")

        # agent data
        @test nagents(other) == nagents(model)
        @test all(haskey(other.agents, i) for i in allids(model))
        @test all(model[i].old_opinion == other[i].old_opinion for i in allids(model))
        @test all(model[i].new_opinion == other[i].new_opinion for i in allids(model))
        @test all(model[i].previous_opinion == other[i].previous_opinion for i in allids(model))
        # properties
        @test model.ϵ == other.ϵ
        # model data
        test_model_data(model, other)

        rm("test.jld2")
    end

    @testset "GridSpace" begin
        # predator_prey used since properties is a NamedTuple, and contains an Array
        model, astep, mstep = Models.predator_prey()
        step!(model, astep, mstep, 50)
        AgentsIO.save_checkpoint("test.jld2", model)
        other = AgentsIO.load_checkpoint("test.jld2"; scheduler = Schedulers.by_property(:type))

        # agent data
        @test nagents(other) == nagents(model)
        @test all(haskey(other.agents, i) for i in allids(model))
        @test all(model[i].type == other[i].type for i in allids(model))
        @test all(model[i].energy == other[i].energy for i in allids(model))
        @test all(model[i].reproduction_prob == other[i].reproduction_prob for i in allids(model))
        @test all(model[i].Δenergy == other[i].Δenergy for i in allids(model))
        # properties
        @test model.fully_grown == other.fully_grown
        @test model.countdown == other.countdown
        @test model.regrowth_time == other.regrowth_time
        # model data
        test_model_data(model, other)
        # space data
        @test typeof(model.space) == typeof(other.space)    # to check periodicity
        test_space(model.space, other.space)

        rm("test.jld2")
    end

    @testset "ContinuousSpace" begin
        model, astep, mstep = Models.social_distancing(N = 300)
        step!(model, astep, mstep, 100)
        AgentsIO.save_checkpoint("test.jld2", model)
        other = AgentsIO.load_checkpoint("test.jld2")

        # agent data
        @test nagents(other) == nagents(model)
        @test all(haskey(other.agents, i) for i in allids(model))
        @test all(model[i].pos == other[i].pos for i in allids(model))
        @test all(model[i].vel == other[i].vel for i in allids(model))
        @test all(model[i].mass == other[i].mass for i in allids(model))
        @test all(model[i].days_infected == other[i].days_infected for i in allids(model))
        @test all(model[i].status == other[i].status for i in allids(model))
        @test all(model[i].β == other[i].β for i in allids(model))
        # properties
        @test model.infection_period == other.infection_period
        @test model.reinfection_probability == other.reinfection_probability
        @test model.detection_time == other.detection_time
        @test model.death_rate == other.death_rate
        @test model.interaction_radius == other.interaction_radius
        @test model.dt == other.dt
        # model data
        test_model_data(model, other)
        # space data
        @test typeof(model.space) == typeof(other.space)    # to check periodicity
        test_space(model.space, other.space)

        rm("test.jld2")
    end

    @testset "GraphSpace" begin
        struct ModelData
            i::Int
            f::Float32
            d::Dict{Int,String}
        end

        model = ABM(
            Agent7,
            GraphSpace(complete_graph(10));
            properties = ModelData(3, 4.2f32, Dict(1 => "foo", 2 => "bar")),
            rng = MersenneTwister(42)
        )

        for i in 1:30
            add_agent_pos!(Agent7(i, i % 10 + 1, rand(model.rng) < 0.5, rand(model.rng, Int)), model)
        end

        AgentsIO.save_checkpoint("test.jld2", model)
        other = AgentsIO.load_checkpoint("test.jld2")

        # agent data
        @test nagents(other) == nagents(model)
        @test all(haskey(other.agents, i) for i in allids(model))
        @test all(model[i].pos == other[i].pos for i in allids(model))
        @test all(model[i].f1 == other[i].f1 for i in allids(model))
        @test all(model[i].f2 == other[i].f2 for i in allids(model))
        # properties
        @test model.i == other.i
        @test model.f == other.f
        @test all(haskey(other.d, k) for k in keys(model.d))
        @test all(other.d[k] == v for (k, v) in model.d)
        # model data
        test_model_data(model, other)
        # space data
        @test model.space.graph == other.space.graph
        @test length(model.space.s) == length(other.space.s)
        @test all(length(model.space.s[pos]) == length(other.space.s[pos]) for pos in eachindex(model.space.s))
        @test all(all(x in other.space.s[pos] for x in model.space.s[pos]) for pos in eachindex(model.space.s))

        rm("test.jld2")
    end

    @testset "Grid Pathfinder" begin
        astep!(a, m) = move_along_route!(a, m, m.pathfinder)
        walk = BitArray(fill(true, 10, 10))
        walk[2, 2] = false
        walk[9, 9] = false
        pmap = abs.(rand(Int, 10, 10)) .% 10
        direct = Pathfinding.DirectDistance{2}([0, 10])
        maxd = Pathfinding.MaxDistance{2}()
        hmm = Pathfinding.PenaltyMap(pmap)

        function setup_model(; kwargs...)
            space = GridSpace((10, 10); periodic = false)
            pathfinder = Pathfinding.AStar(space; kwargs...)
            model = ABM(
                Agent1,
                space;
                properties = (pathfinder = pathfinder,),
                rng = MersenneTwister(42)
            )
            add_agent!((1, 1), model)
            plan_route!(model[1], (10, 10), model.pathfinder)
            step!(model, astep!, dummystep)

            AgentsIO.save_checkpoint("test.jld2", model)
            other = AgentsIO.load_checkpoint("test.jld2")
            return model, other
        end

        test_pathfinding_model(setup_model()...)
        test_pathfinding_model(setup_model(; diagonal_movement = true)...)
        test_pathfinding_model(setup_model(; admissibility = 0.5)...)
        test_pathfinding_model(setup_model(; walkmap = walk)...)
        test_pathfinding_model(setup_model(; cost_metric = direct)...)
        test_pathfinding_model(setup_model(; cost_metric = maxd)...)
        test_pathfinding_model(setup_model(; cost_metric = hmm)...)
        test_pathfinding_model(setup_model(; cost_metric = Pathfinding.PenaltyMap(pmap, direct))...)
        test_pathfinding_model(setup_model(; cost_metric = Pathfinding.PenaltyMap(pmap, maxd))...)
        test_pathfinding_model(setup_model(; cost_metric = Pathfinding.PenaltyMap(pmap, hmm))...)

        rm("test.jld2")
    end

    @testset "Continuous Pathfinder" begin
        astep!(a, m) = move_along_route!(a, m, m.pathfinder, 0.89, 0.56)
        walk = BitArray(fill(true, 10, 10))
        walk[2, 2] = false
        walk[9, 9] = false
        pmap = abs.(rand(Int, 10, 10)) .% 10
        direct = Pathfinding.DirectDistance{2}([0, 10])
        maxd = Pathfinding.MaxDistance{2}()
        hmm = Pathfinding.PenaltyMap(pmap)

        function setup_model(; kwargs...)
            space = ContinuousSpace((10.0, 10.0); periodic = false)
            pathfinder = Pathfinding.AStar(space; kwargs...)
            model = ABM(
                Agent6,
                deepcopy(space);
                properties = (pathfinder = pathfinder,),
                rng = MersenneTwister(42)
            )
            add_agent!((1.3, 1.5), model, (0.0, 0.0), 0.0)
            plan_route!(model[1], (9.7, 4.8), model.pathfinder)
            step!(model, astep!, dummystep)

            AgentsIO.save_checkpoint("test.jld2", model)
            other = AgentsIO.load_checkpoint("test.jld2")
            return model, other
        end

        test_pathfinding_model(setup_model(walkmap = trues(10, 10))...)
        test_pathfinding_model(setup_model(walkmap = trues(10, 10), admissibility = 0.5)...)
        test_pathfinding_model(setup_model(walkmap = walk)...)
        test_pathfinding_model(setup_model(walkmap = trues(10, 10), cost_metric = direct)...)
        test_pathfinding_model(setup_model(walkmap = trues(10, 10), cost_metric = maxd)...)
        test_pathfinding_model(setup_model(walkmap = trues(10, 10), cost_metric = hmm)...)
        test_pathfinding_model(setup_model(cost_metric = Pathfinding.PenaltyMap(pmap, direct))...)
        test_pathfinding_model(setup_model(cost_metric = Pathfinding.PenaltyMap(pmap, maxd))...)
        test_pathfinding_model(setup_model(cost_metric = Pathfinding.PenaltyMap(pmap, hmm))...)

        rm("test.jld2")
    end

    @testset "Multi-agent" begin
        model, _ = Models.daisyworld()
        AgentsIO.save_checkpoint("test.jld2", model)
        other = @test_nowarn AgentsIO.load_checkpoint("test.jld2"; scheduler = Models.daisysched, warn = false)

        # agent data
        @test nagents(other) == nagents(model)
        @test all(haskey(other.agents, i) for i in allids(model))
        @test all(model[i].pos == other[i].pos for i in allids(model))
        @test all(model[i].temperature == other[i].temperature for i in allids(model) if model[i] isa Models.Land)
        @test all(model[i].breed == other[i].breed for i in allids(model) if model[i] isa Models.Daisy)
        @test all(model[i].age == other[i].age for i in allids(model) if model[i] isa Models.Daisy)
        @test all(model[i].albedo == other[i].albedo for i in allids(model) if model[i] isa Models.Daisy)
        # properties
        @test model.max_age == other.max_age
        @test model.surface_albedo == other.surface_albedo
        @test model.solar_luminosity == other.solar_luminosity
        @test model.solar_change == other.solar_change
        @test model.scenario == other.scenario
        @test model.tick == other.tick
        # model data
        test_model_data(model, other)
        # space data
        test_space(model.space, other.space)

        rm("test.jld2")
    end

    @testset "OSMSpace" begin
        @agent Zombie OSMAgent begin
            infected::Bool
        end
        model = ABM(Zombie, OpenStreetMapSpace(OSM.test_map()); rng = MersenneTwister(42))

        for id in 1:100
            start = random_position(model)
            finish = OSM.random_road_position(model)
            human = Zombie(id, start, false)
            add_agent_pos!(human, model)
            plan_route!(human, finish, model)
        end

        start = OSM.road((51.530876112711745, 9.945125635913511), model)
        finish = OSM.intersection((51.5328328, 9.9351811), model)
        zombie = add_agent!(start, model, true)
        plan_route!(zombie, finish, model)

        AgentsIO.save_checkpoint("test.jld2", model)
        @test_throws AssertionError AgentsIO.load_checkpoint("test.jld2")
        other = AgentsIO.load_checkpoint("test.jld2"; map = OSM.test_map())

        # agent data
        @test nagents(other) == nagents(model)
        @test all(haskey(other.agents, i) for i in allids(model))
        @test all(OSM.latlon(model[i].pos, model) == OSM.latlon(other[i].pos, other) for i in allids(model))
        @test all(model[i].infected == other[i].infected for i in allids(model))
        # model data
        test_model_data(model, other)
        @test sort(collect(keys(model.space.routes))) == sort(collect(keys(other.space.routes)))
        @test all(model.space.routes[i].route == other.space.routes[i].route for i in keys(model.space.routes))
        @test all(model.space.routes[i].start == other.space.routes[i].start for i in keys(model.space.routes))
        @test all(model.space.routes[i].dest == other.space.routes[i].dest for i in keys(model.space.routes))
        @test all(model.space.routes[i].return_route == other.space.routes[i].return_route for i in keys(model.space.routes))
        @test all(model.space.routes[i].has_to_return == other.space.routes[i].has_to_return for i in keys(model.space.routes))

        rm("test.jld2")
    end
end


text:
export analytical_jacobian_strength!,
       analytical_jacobian_strength,
       analytical_jacobian_position!,
       analytical_jacobian_position,
       analytical_jacobian_pressure

function analytical_jacobian_strength!(dpdS, dpdSstar, Css, Cts, Ctsblob, wtarget, target::Vector{ComplexF64}, source::T, freestream,
	                                   idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t;
	                                   iscauchystored::Bool = false, issourcefixed::Bool = false) where T <: Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}
	Ny = size(target, 1)
	Nx = size(source, 1)

	δ = source[1].δ

	U = freestream.U
	cst = inv(2*π)

	# @assert size(dpdS) == (Ny, Nx) && size(dpdSstar) == (Ny, Nx) && size(Css) == (Nx, Nx) && size(Cts) == (Ny, Nx) && size(wtarget) == (Ny,)

	fill!(dpdS, complex(0.0)); fill!(dpdSstar, complex(0.0))

	if iscauchystored == false
	  fill!(Css, complex(0.0));
	  fill!(Cts, complex(0.0));
	  fill!(Ctsblob, complex(0.0));

	  # Account for the presence of the free-stream
	  fill!(wtarget, conj(U))

	  # Construct the Cauchy matrix to store 1/(zJ-zK) in Css
	  @inbounds for K=1:Nx
		  zK = source[K].z
		  # exploit anti-symmetry of the velocity kernel
		  # Diagonal entries are equal to 0
		  for J=K+1:Nx
			  zJ = source[J].z
			  tmp = abs2(zJ - zK) + δ^2
			  Css[J,K] = conj(zJ - zK)/tmp
			  Css[K,J] = -Css[J,K]
		  end
	  end

	  # Construct the Cauchy matrix to store (z̄-z̄J)/(|z-zJ|^2 + δ^2) and the total induced velocity in Cts and wtarget
	  @inbounds for J=1:Nx
		  zJ = source[J].z
		  SJ = -im*source[J].S
		  for i=1:Ny
			  zi = target[i]
			  tmp = abs2(zi - zJ) + δ^2
			  Cts[i,J] = inv(zi - zJ)
			  Ctsblob[i,J] = conj(zi - zJ)/tmp
			  wtarget[i] += cst*SJ*Ctsblob[i,J]
		  end
	  end
	end

	# Evaluate ∂(-0.5v^2)/∂SL, ∂(-0.5v^2)/∂S̄L
	for L=1:Nx
	   zL = source[L].z
	   SL = -im*source[L].S
	   for i=1:Ny
	       wt = wtarget[i]
	       dpdS[i,L] = -0.5*conj(wt)*cst*Ctsblob[i,L]
	       # dpdSstar[i,L] = conj(dpdS[i,L])
	   end
	end

	# Evaluate ∂(-∂ϕ/∂t)/∂SL, ∂(-∂ϕ/∂t)/∂S̄L
	@inbounds for L=1:Nx
		SL = -im*source[L].S
		for K=1:Nx
			SK = -im*source[K].S
			for i=1:Ny
				term1 = 0.5*cst^2*conj(SK)*Cts[i,L]*conj(Css[L,K])
				dpdS[i,L] += term1
				# dpdSstar[i,L] += conj(term1)
			end
		end
		# add the freestream contribution
		if abs(U) > 100*eps()
			for i=1:Ny
				term3 = 0.5*cst*Cts[i,L]*U
				dpdS[i,L] += term3
				# dpdSstar[i,L] += conj(term3)
			end
		end

	    for J=1:Nx
	        if issourcefixed == true && typeof(source[J]) <: PotentialFlow.Blobs.Blob{ComplexF64, Float64}
	            continue
	        else
	        	SJ = -im*source[J].S
	        end

	  		for i=1:Ny
	           term2 = 0.5*cst^2*conj(SJ)*conj(Cts[i,J])*Css[J,L]
	           dpdS[i,L] += term2
	           # dpdSstar[i,L] += conj(term2)
	       end
	   end
	end
	dpdSstar .= conj(dpdS)
	nothing
end

"""
Returns the Jacobian of the pressure computed from the unsteady Bernoulli equation with respect to the strength and conjugate strength of the singularities.
Note that we multiply the strength `point.S` of a singularity by -i to move from the convention Γ+iQ (used in PotentialFlow.jl) to S = Q-iΓ.
"""
function analytical_jacobian_strength!(dpdS, dpdSstar, Css, Cts, wtarget, target::Vector{ComplexF64}, source::T, freestream,
	                                   idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t;
                                       iscauchystored::Bool = false) where T <: Vector{PotentialFlow.Points.Point{Float64, Float64}}
	Ny = size(target, 1)
    Nx = size(source, 1)

    U = freestream.U
    cst = inv(2*π)

    # @assert size(dpdS) == (Ny, Nx) && size(dpdSstar) == (Ny, Nx) && size(Css) == (Nx, Nx) && size(Cts) == (Ny, Nx) && size(wtarget) == (Ny,)

    fill!(dpdS, complex(0.0)); fill!(dpdSstar, complex(0.0))

    if iscauchystored == false

        fill!(Css, complex(0.0)); fill!(Cts, complex(0.0))

        # Account for the presence of the free-stream
        fill!(wtarget, conj(U))

        # Construct the Cauchy matrix to store 1/(zJ-zK) in Css
        @inbounds for K=1:Nx
            zK = source[K].z
            # exploit anti-symmetry of the velocity kernel
            for J=K+1:Nx
                zJ = source[J].z
                Css[J,K] = inv(zJ - zK)
                Css[K,J] = -Css[J,K]
            end
        end

        # Construct the Cauchy matrix to store 1/(z-zJ) and the total induced velocity in Cts and wtarget
        for J=1:Nx
            zJ = source[J].z
            SJ = -im*source[J].S
            for i=1:Ny
                zi = target[i]
                Cts[i,J] = inv(zi - zJ)
                wtarget[i] += Cts[i,J]*cst*SJ
            end
        end
    end

    # Evaluate ∂(-0.5v^2)/∂SL, ∂(-0.5v^2)/∂S̄L
    for L in idx #1:Nx
        zL = source[L].z
        SL = -im*source[L].S
        for i=1:Ny
            wt = wtarget[i]
            dpdS[i,L] = -0.5*conj(wt)*cst*Cts[i,L]
            # dpdSstar[i,L] = conj(dpdS[i,L])
        end
    end

    # Evaluate ∂(-∂ϕ/∂t)/∂SL, ∂(-∂ϕ/∂t)/∂S̄L
    @inbounds for L in idx #1:Nx
		SL = -im*source[L].S
		for K=1:Nx
			SK = -im*source[K].S
			for i=1:Ny
				term1 = 0.5*cst^2*conj(SK)*Cts[i,L]*conj(Css[L,K])
				dpdS[i,L] += term1
				# dpdSstar[i,L] += conj(term1)
			end
		end
		# add the freestream contribution
		if abs(U)>100*eps()
			for i=1:Ny
				term3 = 0.5*cst*Cts[i,L]*U
				dpdS[i,L] += term3
				# dpdSstar[i,L] += conj(term3)
			end
		end

	    for J=1:Nx
	        SJ = -im*source[J].S
       		for i=1:Ny
                term2 = 0.5*cst^2*conj(SJ)*conj(Cts[i,J])*Css[J,L]
                dpdS[i,L] += term2
                # dpdSstar[i,L] += conj(term2)
            end
        end
    end
	dpdSstar .= conj(dpdS)
    nothing
end

analytical_jacobian_strength!(dpdS, dpdSstar, Css, Cts, Ctsblob, wtarget, target::Vector{ComplexF64}, source, freestream, t;
	                          iscauchystored::Bool = false) =
analytical_jacobian_strength!(dpdS, dpdSstar, Css, Cts, Ctsblob, wtarget, target, source, freestream, 1:length(source), t;
							  iscauchystored = iscauchystored)

function analytical_jacobian_strength(target::Vector{ComplexF64}, source, freestream,
	                                  idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t;
                                      iscauchystored::Bool = false)
	Ny = size(target, 1)
	Nx = size(source, 1)
	Css = zeros(ComplexF64, Nx, Nx)
	Cts = zeros(ComplexF64, Ny, Nx)
	wtarget = zeros(ComplexF64, Ny)

	dpdS = zeros(ComplexF64, Ny, Nx)
	dpdSstar = zeros(ComplexF64, Ny, Nx)
	if typeof(source)<:Vector{PotentialFlow.Points.Point{Float64, Float64}}
		analytical_jacobian_strength!(dpdS, dpdSstar, Css, Cts, wtarget, target, source, freestream, idx, t;
	                                  iscauchystored = iscauchystored)
	elseif typeof(source)<:Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}
		Ctsblob = zeros(ComplexF64, Ny, Nx)
		analytical_jacobian_strength!(dpdS, dpdSstar, Css, Cts, Ctsblob, wtarget, target, source, freestream, idx, t;
									  iscauchystored = iscauchystored)
	end

	return dpdS, dpdSstar
end

analytical_jacobian_strength(target::Vector{ComplexF64}, source, freestream, t;
                                      iscauchystored::Bool = false) =
analytical_jacobian_strength(target, source, freestream, 1:length(source), t;
                                      iscauchystored = iscauchystored)

# Version for point singularities
"""
Returns the Jacobian of the pressure computed from the unsteady Bernoulli equation with respect to the position and conjugate position of the singularities.
Note that we multiply the strength `point.S` of a singularity by -i to move from the convention Γ+iQ (used in PotentialFlow.jl) to Q-iΓ.
"""
function analytical_jacobian_position!(dpdz, dpdzstar, Css, Cts, wtarget, target::Vector{ComplexF64}, source::T, freestream,
	                                   idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t;
                                       iscauchystored::Bool = false) where T <: Vector{PotentialFlow.Points.Point{Float64, Float64}}
    Nx = size(source, 1)
	Ny = size(target, 1)

    U = freestream.U
    cst = inv(2*π)

    # @assert size(dpdz) == (Ny, Nx) && size(dpdzstar) == (Ny, Nx) && size(Css) == (Nx, Nx) && size(Cts) == (Ny, Nx) && size(wtarget) == (Ny,)

    fill!(dpdz, complex(0.0)); fill!(dpdzstar, complex(0.0))

    if iscauchystored == false
        fill!(Css, complex(0.0)); fill!(Cts, complex(0.0))

        # Account for the presence of the free-stream
        fill!(wtarget, conj(U))

        # Construct the Cauchy matrix to store 1/(zJ-zK) in Css
        @inbounds for K=1:Nx
            zK = source[K].z
            # exploit anti-symmetry of the velocity kernel
            for J=K+1:Nx
                zJ = source[J].z
                Css[J,K] = inv(zJ - zK)
                Css[K,J] = -Css[J,K]
            end
        end

        # Construct the Cauchy matrix to store 1/(z-zJ) and the total induced velocity in Cts and wtarget
        @inbounds for J=1:Nx
            zJ = source[J].z
            SJ = -im*source[J].S
            for i=1:Ny
                zi = target[i]
                Cts[i,J] = inv(zi - zJ)
                wtarget[i] += Cts[i,J]*cst*SJ
            end
        end
    end
	# @show norm(Cts)
	# @show norm(real(Css)), norm(imag(Css)), norm(Css), norm(real(wtarget))

    # Evaluate ∂(-0.5v^2)/∂zL, ∂(-0.5v^2)/∂z̄L
    @inbounds for L in idx #1:Nx
        zL = source[L].z
        SL = -im*source[L].S
        for i=1:Ny
            wt = wtarget[i] #+ conj(U)
            dpdz[i,L] = -0.5*conj(wt)*SL*cst*Cts[i,L]^2
            # dpdzstar[i,L] = conj(dpdz[i,L])
        end
    end

    # Evaluate ∂(-∂ϕ/∂t)/∂zL, ∂(-∂ϕ/∂t)/∂z̄L
	@inbounds for L in idx #1:Nx
		SL = -im*source[L].S
		# @show SL
		for K=1:Nx
			SK = -im*source[K].S
			for i=1:Ny
				term1 = 0.5*cst^2*SL*conj(SK)*(Cts[i,L])^2*conj(Css[L,K])
				term1 += 0.5*cst^2*conj(SL)*SK*conj(Cts[i,L])*(-Css[L,K]^2)
				dpdz[i,L] += term1
				# dpdzstar[i,L] += conj(term1)
			end
		end
		if abs(U) > 100*eps()
			for i=1:Ny
				term3 = 0.5*cst*SL*Cts[i,L]^2*U
				dpdz[i,L] += term3
				# dpdzstar[i,L] += conj(term3)
			end
		end

		for J=1:Nx
            SJ = -im*source[J].S
			for i=1:Ny
                term2 = 0.5*cst^2*conj(SJ)*SL*conj(Cts[i,J])*(Css[J,L])^2
                dpdz[i,L] += term2
                # dpdzstar[i,L] += conj(term2)
            end
        end
	end
	dpdzstar .= conj(dpdz)
    nothing
end

# Version for regularized point singularities (aka blobs)
"""
Returns the Jacobian of the pressure computed from the unsteady Bernoulli equation with respect to the position and conjugate position of the singularities.
Note that we multiply the strength `point.S` of a singularity by -i to move from the convention Γ+iQ (used in PotentialFlow.jl) to Q-iΓ.
"""
function analytical_jacobian_position!(dpdz, dpdzstar, Css, ∂Css, Cts, Ctsblob, ∂Ctsblob, wtarget, target::Vector{ComplexF64}, source::T, freestream,
									   idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t;
                                       iscauchystored::Bool = false) where T <: Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}
	Ny = size(target, 1)
	Nx = size(source, 1)

	δ = source[1].δ

	U = freestream.U
	cst = inv(2*π)

	# @assert size(dpdz) == (Ny, Nx) && size(dpdzstar) == (Ny, Nx) && size(Css) == (Nx, Nx) && size(Cts) == (Ny, Nx) && size(wtarget) == (Ny,)

	fill!(dpdz, complex(0.0)); fill!(dpdzstar, complex(0.0))

	if iscauchystored == false
	   fill!(Css, complex(0.0));
	   fill!(∂Css, 0.0);
	   fill!(Cts, complex(0.0));
	   fill!(Ctsblob, complex(0.0));
	   fill!(∂Ctsblob, 0.0)

	   # Account for the presence of the free-stream
	   fill!(wtarget, conj(U))

	   # Construct the Cauchy matrix to store 1/(zJ-zK) in Css
	   @inbounds for K=1:Nx
	       zK = source[K].z
	       # exploit anti-symmetry of the velocity kernel
		   # Diagonal entries are equal to 0
	       for J=K+1:Nx
	           zJ = source[J].z
			   tmp = abs2(zJ - zK) + δ^2
	           Css[J,K] = conj(zJ - zK)/tmp
			   Css[K,J] = -Css[J,K]
			   ∂Css[J,K] = -(δ/tmp)^2
			   ∂Css[K,J] =  ∂Css[J,K]
	       end
	   end

	   # Construct the Cauchy matrix to store (z̄-z̄J)/(|z-zJ|^2 + δ^2) and the total induced velocity in Cts and wtarget
	   @inbounds for J=1:Nx
	       zJ = source[J].z
	       SJ = -im*source[J].S
	       for i=1:Ny
	           zi = target[i]
			   tmp = abs2(zi - zJ) + δ^2
			   Cts[i,J] = inv(zi - zJ)
	           Ctsblob[i,J] = conj(zi - zJ)/tmp
			   ∂Ctsblob[i,J] = -(δ/tmp)^2
	           wtarget[i] += cst*SJ*Ctsblob[i,J]
	       end
	   end
	end
	# @show 0.5*(wtarget .- conj(U))*U + 0.5*conj(wtarget .- conj(U))*conj(U)

	# Evaluate ∂(-0.5v^2)/∂zL, ∂(-0.5v^2)/∂z̄L
	@inbounds for L in idx #1:Nx
	   zL = source[L].z
	   SL = -im*source[L].S
	   for i=1:Ny
	       wt = wtarget[i]
		   dpdz[i,L] = -0.5*cst*conj(SL)*∂Ctsblob[i,L]*wt
	       dpdz[i,L] += -0.5*conj(wt)*cst*SL*Ctsblob[i,L]^2
	       # dpdzstar[i,L] = conj(dpdz[i,L])
	   end
	end

	# Evaluate ∂(-∂ϕ/∂t)/∂zL, ∂(-∂ϕ/∂t)/∂z̄L
	@inbounds for L in idx #1:Nx
		SL = -im*source[L].S
		for K=1:Nx
			SK = -im*source[K].S
			for i=1:Ny
				termJL = 0.5*cst^2*SL*conj(SK)*Cts[i,L]^2*conj(Css[L,K])
				termJL += 0.5*cst^2*SL*conj(SK)*Cts[i,L]*(-∂Css[L,K])
				termJL += 0.5*cst^2*conj(SL)*SK*conj(Cts[i,L])*(-Css[L,K]^2)
				dpdz[i,L] += termJL
				# dpdzstar[i,L] += conj(termJL)
			end
		end
		if abs(U) > 100*eps()
			for i=1:Ny
				termU = 0.5*cst*SL*Cts[i,L]^2*U
				dpdz[i,L] += termU

				# dpdzstar[i,L] += conj(termU)
			end
		end

		for J=1:Nx
		SJ = -im*source[J].S
			for i=1:Ny
	           termKL = 0.5*cst^2*SJ*conj(SL)*Cts[i,J]*∂Css[J,L]
			   termKL += 0.5*cst^2*conj(SJ)*SL*conj(Cts[i,J])*Css[J,L]^2
	           dpdz[i,L] += termKL

	           # dpdzstar[i,L] += conj(termKL)
		   end
	   end
	end
	dpdzstar .= conj(dpdz)
	nothing
end

analytical_jacobian_position!(dpdz, dpdzstar, Css, Cts, wtarget, target::Vector{ComplexF64}, source::T, freestream,
	                                   idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t;
                                       iscauchystored::Bool = false) where T =
analytical_jacobian_position!(dpdz, dpdzstar, Css, Cts, wtarget, target, source, freestream,
	                                  1:length(source), t;
                                      iscauchystored = iscauchystored)

function analytical_jacobian_position(target::Vector{ComplexF64}, source, freestream,
									  idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t;
                                      iscauchystored::Bool = false)
    Ny = size(target, 1)
    Nx = size(source, 1)
    Css = zeros(ComplexF64, Nx, Nx)
    Cts = zeros(ComplexF64, Ny, Nx)
    wtarget = zeros(ComplexF64, Ny)

    dpdz = zeros(ComplexF64, Ny, Nx)
    dpdzstar = zeros(ComplexF64, Ny, Nx)

	if typeof(source)<:Vector{PotentialFlow.Points.Point{Float64, Float64}}
	    analytical_jacobian_position!(dpdz, dpdzstar, Css, Cts, wtarget, target, source, freestream, idx, t;
	                                  iscauchystored = iscauchystored)
	elseif typeof(source)<:Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}
		∂Css = zeros(Nx, Nx)
		Ctsblob = zeros(ComplexF64, Ny, Nx)
		∂Ctsblob = zeros(Ny, Nx)
		analytical_jacobian_position!(dpdz, dpdzstar, Css, ∂Css, Cts, Ctsblob, ∂Ctsblob, wtarget, target, source, freestream, idx, t;
									  iscauchystored = iscauchystored)
	end

    return dpdz, dpdzstar
end

analytical_jacobian_position(target::Vector{ComplexF64}, source, freestream, t;
                                      iscauchystored::Bool = false) =
analytical_jacobian_position(target, source, freestream, 1:length(source), t;
                                      iscauchystored = iscauchystored)

analytical_jacobian_pressure(target, source, freestream, t) = analytical_jacobian_pressure(target, source, freestream, 1:length(source), t)

# Version with allocations for point vortices
function analytical_jacobian_pressure(target, source::T, freestream, idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t) where T <: Vector{PotentialFlow.Points.Point{Float64, Float64}}
	Nv = size(source, 1)
	Nx = 3*Nv
	Ny = size(target, 1)

	J = zeros(Ny, Nx)

	wtarget = zeros(ComplexF64, Ny)

	dpd = zeros(ComplexF64, Ny, Nv)
	dpdstar = zeros(ComplexF64, Ny, Nv)

	Css = zeros(ComplexF64, Nv, Nv)
	Cts = zeros(ComplexF64, Ny, Nv)

	analytical_jacobian_pressure!(J, wtarget, dpd, dpdstar, Css, Cts, target, source, freestream, idx, t)
	return J
end

# Version with allocations for regularized vortices
function analytical_jacobian_pressure(target, source::T, freestream, idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t) where T <: Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}
	Nv = size(source, 1)
	Nx = 3*Nv
	Ny = size(target, 1)

	J = zeros(Ny, Nx)

	wtarget = zeros(ComplexF64, Ny)

	dpd = zeros(ComplexF64, Ny, Nv)
	dpdstar = zeros(ComplexF64, Ny, Nv)

	Css = zeros(ComplexF64, Nv, Nv)
	Cts = zeros(ComplexF64, Ny, Nv)

	∂Css = zeros(Nv, Nv)
	Ctsblob = zeros(ComplexF64, Ny, Nv)
	∂Ctsblob = zeros(Ny, Nv)

	analytical_jacobian_pressure!(J, wtarget, dpd, dpdstar, Css, Cts, ∂Css, Ctsblob, ∂Ctsblob, target, source, freestream, idx, t)
	return J
end

# In-place version for point vortices
function analytical_jacobian_pressure!(J, wtarget, dpd, dpdstar, Css, Cts, target, source::T, freestream, idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t) where T <: Vector{PotentialFlow.Points.Point{Float64, Float64}}
	Nv = size(source, 1)
	Nx = 3*Nv
	Ny = size(target, 1)

	@assert size(J) == (Ny, Nx)
	@assert size(wtarget) == (Ny,)
	@assert size(dpd) == (Ny, Nv)
	@assert size(dpdstar) == (Ny, Nv)
	@assert size(Css) == (Nv, Nv)
	@assert size(Cts) == (Ny, Nv)

	fill!(J, 0.0)
	fill!(wtarget, 0.0*im)
	fill!(dpd, 0.0*im)
	fill!(dpdstar, 0.0*im)
	fill!(Css, 0.0*im)
	fill!(Cts, 0.0*im)

	cst = inv(2*π)

	src = vcat(source...)

	# Account for the presence of the free-stream
	U = freestream.U
	fill!(wtarget, conj(U))

	# Construct the Cauchy matrix to store 1/(zJ-zK) in Css
	@inbounds for K=1:Nv
		zK = src[K].z
		# exploit anti-symmetry of the velocity kernel
		# Diagonal entries are equal to 0
		for J=K+1:Nv
			zJ = src[J].z
			Css[J,K] = inv(zJ - zK)
			Css[K,J] = -Css[J,K]
		end
	end

	# Construct the Cauchy matrix to store (z̄-z̄J)/(|z-zJ|^2 + δ^2) and the total induced velocity in Cts and wtarget
	@inbounds for J=1:Nv
		zJ = src[J].z
		SJ = -im*src[J].S
		for i=1:Ny
			zi = target[i]
			Cts[i,J] = inv(zi - zJ)
			wtarget[i] += cst*SJ*Cts[i,J]
		end
	end

	analytical_jacobian_position!(dpd, dpdstar, Css, Cts, wtarget, target, src, freestream, idx, t;
								  iscauchystored = true)

	# Fill dpdpx and dpdy
	# J[:, 1:3:3*(Nv-1)+1] .= 2*real.(dpd[:,1:Nv])
	# J[:, 2:3:3*(Nv-1)+2] .= -2*imag.(dpd[:,1:Nv])
	@inbounds for L in idx
		J[:, 3*(L-1)+1] .= 2*real.(view(dpd,:,L))
		J[:, 3*(L-1)+2] .= -2*imag.(view(dpd,:,L))
	end

	analytical_jacobian_strength!(dpd, dpdstar, Css, Cts, wtarget, target, src, freestream, idx, t;
								  iscauchystored = true)

	@inbounds for L in idx
  		J[:, 3*(L-1)+3] .=  2*imag.(view(dpd,:,L))
	end
	# Vortices
	# J[:, 3:3:3*(Nv-1)+3] .= 2imag.(dpd[:,1:Nv])

	return J
end

# In-place version for regularized vortices
function analytical_jacobian_pressure!(J, wtarget, dpd, dpdstar, Css, Cts, ∂Css, Ctsblob, ∂Ctsblob, target, source::T, freestream, idx::Union{Int64, Vector{Int64}, UnitRange{Int64}}, t) where T <: Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}
    Nv = size(source, 1)
    Nx = 3*Nv
    Ny = size(target, 1)

	@assert size(J) == (Ny, Nx)
	@assert size(wtarget) == (Ny,)
	@assert size(dpd) == (Ny, Nv)
	@assert size(dpdstar) == (Ny, Nv)
	@assert size(Css) == (Nv, Nv)
	@assert size(Cts) == (Ny, Nv)
	@assert size(∂Css) == (Nv, Nv)
	@assert typeof(∂Css) <: Matrix{Float64}
	@assert size(Ctsblob) == (Ny, Nv)
	@assert typeof(Ctsblob) <: Matrix{ComplexF64}
	@assert size(∂Ctsblob) == (Ny, Nv)
	@assert typeof(∂Ctsblob) <: Matrix{Float64}

	fill!(J, 0.0)
	fill!(wtarget, 0.0*im)
	fill!(dpd, 0.0*im)
	fill!(dpdstar, 0.0*im)
	fill!(Css, 0.0*im)
	fill!(Cts, 0.0*im)
	fill!(∂Css, 0.0)
	fill!(Ctsblob, 0.0*im)
	fill!(∂Ctsblob, 0.0)

    cst = inv(2*π)

    src = vcat(source...)
	δ = src[1].δ

	# Account for the presence of the free-stream
	U = freestream.U
	fill!(wtarget, conj(freestream.U))

	# Construct the Cauchy matrix to store 1/(zJ-zK) in Css
	@inbounds for K=1:Nv
		zK = src[K].z
		# exploit anti-symmetry of the velocity kernel
		# Diagonal entries are equal to 0
		for J=K+1:Nv
			zJ = src[J].z
			tmp = abs2(zJ - zK) + δ^2
			Css[J,K] = conj(zJ - zK)/tmp
			Css[K,J] = -Css[J,K]
			∂Css[J,K] = -(δ/tmp)^2
			∂Css[K,J] =  ∂Css[J,K]
		end
	end

	# Construct the Cauchy matrix to store (z̄-z̄J)/(|z-zJ|^2 + δ^2) and the total induced velocity in Cts and wtarget
	@inbounds for J=1:Nv
		zJ = src[J].z
		SJ = -im*src[J].S
		for i=1:Ny
			zi = target[i]
			tmp = abs2(zi - zJ) + δ^2
			Cts[i,J] = inv(zi - zJ)
			Ctsblob[i,J] = conj(zi - zJ)/tmp
			∂Ctsblob[i,J] = -(δ/tmp)^2
			wtarget[i] += cst*SJ*Ctsblob[i,J]
		end
	end

	analytical_jacobian_position!(dpd, dpdstar, Css, ∂Css, Cts, Ctsblob, ∂Ctsblob, wtarget, target, src, freestream, idx, t;
								  iscauchystored = true)

	# Fill dpdpx and dpdy
    # J[:, 1:3:3*(Nv-1)+1] .= 2*real.(dpd[:,1:Nv])
    # J[:, 2:3:3*(Nv-1)+2] .= -2*imag.(dpd[:,1:Nv])
	@inbounds for L in idx
		J[:, 3*(L-1)+1] .= 2*real.(view(dpd,:,L))
		J[:, 3*(L-1)+2] .= -2*imag.(view(dpd,:,L))
	end

	analytical_jacobian_strength!(dpd, dpdstar, Css, Cts, Ctsblob, wtarget, target, src, freestream, idx, t;
								  iscauchystored = true)

	@inbounds for L in idx
		J[:, 3*(L-1)+3] .=  2*imag.(view(dpd,:,L))
	end

    # Vortices
    # J[:, 3:3:3*(Nv-1)+3] .= 2imag.(dpd[:,1:Nv])
    return J
end

# function analytical_jacobian_pressure(target, source::T, freestream, t; issourcefixed::Bool=false) where T <: Union{Vector{PotentialFlow.Points.Point{Float64, Float64}}, Vector{PotentialFlow.Points.Point{ComplexF64, Float64}}, Vector{PotentialFlow.Points.Point{S, Float64} where S<:Number}}
#
#
#     Nv = size(source[1], 1)
#     Ns = size(source[2], 1)
#     Nsing =  Nv + Ns
#     Nx = 3*Nv + Ns
#     Ny = size(target, 1)
#
#     J = zeros(Ny, Nx)
#     cst = inv(2*π)
#
#     wtarget = fill(conj(freestream.U), Ny)
#
#     dpd = zeros(ComplexF64, Ny, Nsing)
#     dpdstar = zeros(ComplexF64, Ny, Nsing)
#
#     Css = zeros(ComplexF64, Nsing, Nsing)
#     Cts = zeros(ComplexF64, Ny, Nsing)
#
#     src = vcat(source...)
#
#     # Construct the Cauchy matrix to store 1/(zJ-zK) in Css
#     @inbounds for K=1:Nsing
#         zK = src[K].z
#         # exploit anti-symmetry of the velocity kernel
#         for J=K+1:Nsing
#             zJ = src[J].z
#             Css[J,K] = inv(zJ - zK)
#             Css[K,J] = -Css[J,K]
#         end
#     end
#
#     # Construct the Cauchy matrix to store 1/(z-zJ) and the total induced velocity in Cts and wtarget
#     @avx for J=1:Nsing
#         zJ = src[J].z
#         SJ = -im*src[J].S
#         for i=1:Ny
#             zi = target[i]
#             Cts[i,J] = inv(zi - zJ)
#             wtarget[i] += Cts[i,J]*cst*SJ
#         end
#     end
#
#     analytical_jacobian_position!(dpd, dpdstar, Css, Cts, wtarget, target, src, freestream, t; iscauchystored = true, issourcefixed = issourcefixed)
#
# 	# Fill dpdpx and dpdy
#     J[:, 1:3:3*(Nv-1)+1] .= 2*real.(dpd[:,1:Nv])
#     J[:, 2:3:3*(Nv-1)+2] .= -2*imag.(dpd[:,1:Nv])
#
# 	analytical_jacobian_strength!(dpd, dpdstar, Css, Cts, wtarget, target, src, freestream, t; iscauchystored = true, issourcefixed = issourcefixed)
#
#     # Vortices
#     J[:, 3:3:3*(Nv-1)+3] .= 2imag.(dpd[:,1:Nv])
#     # # Sources
#     J[:, 3*Nv+1:3*Nv+Ns] .= 2*real.(dpd[:,Nv+1:Nv+Ns])
#     return J
# end

ext: jl


text:
function nlsolve2!(F::OnceDiffed, x, approach::Tuple{<:Newton, <:Dogleg}; maxiter=200, f_abstol=1e-8, f_reltol=1e-12)
    scheme, linesearch = approach

    function normed_residual(Jx, Fx, x)
        Fx, Jx = F(x, Fx, Jx)
        f = (norm(Fx)^2)/2
        if Jx isa Nothing
            if Fx isa Nothing
                throw(ErrorException("Error happened in trust region nlsolve"))
            end
            return f, Jx'*Fx
        else
            return f, Jx'*Fx, Jx'*Jx
        end
    end
    td = TwiceDiffed(normed_residual)
    res = minimize!(td, x, approach)
    return res[1], res[3], res[4]
end


text:
export DefConst

struct DefConst{T} <: AbstractPhysicsConst{T}
    quantity::String
    value::T
    unit::String
end

value(x::DefConst) = x.value
eltype(::Type{DefConst{T}}) where T = T


text:
VERSION >= v"0.4.0-dev+6641" && __precompile__()
module ELF
    include("constants.jl")
    using StructIO
    using FileIO
    using ObjFileBase
    import Base: start, next, done, endof, length, getindex
    import Base: read, readuntil, write, readbytes, seek, seekstart, position
    import Base: show, showcompact
    import Base: sizeof
    import ObjFileBase: readmeta, debugsections, deref, sectionoffset, sectionaddress,
        sectionsize, Section, endianness, replace_sections_from_memory, strtab_lookup,
        getSectionLoadAddress, sectionname, load_strtab, handle, symname, isundef,
        symbolvalue, symbolnum, intptr
    import StructIO: unpack

    abstract ELFFile

    #
    # Represents the actual ELF file
    #
    immutable ELFHandle{T<:IO} <: ObjectHandle
        io::T
        start::Int
        file::ELFFile
    end
    __init__() = push!(ObjFileBase.ObjHandles, ELFHandle)
    Base.eof(handle::ELFHandle) = eof(handle.io)
    ObjFileBase.handle(handle::ELFHandle) = handle

    abstract ELFHeader
    abstract ELFSectionHeader <: Section{ELFHandle}
    abstract ELFSymtabEntry
    abstract ELFRel
    abstract ELFRela
    abstract ELFProgramHeader

    sheader(::Type{ELFHeader}) = ELFSectionHeader
    pheader(::Type{ELFHeader}) = ELFProgramHeader
    symtype(::Type{ELFSectionHeader}) = ELFSymtabEntry

    module ELF32
        import ELF
        using StructIO

        @struct immutable Header <: ELF.ELFHeader
            e_type::UInt16
            e_machine::UInt16
            e_version::UInt32
            e_entry::UInt32
            e_phoff::UInt32
            e_shoff::UInt32
            e_flags::UInt32
            e_ehsize::UInt16
            e_phentsize::UInt16
            e_phnum::UInt16
            e_shentsize::UInt16
            e_shnum::UInt16
            e_shstrndx::UInt16
        end

        @struct immutable SectionHeader <: ELF.ELFSectionHeader
            sh_name::UInt32
            sh_type::UInt32
            sh_flags::UInt32
            sh_addr::UInt32
            sh_offset::UInt32
            sh_size::UInt32
            sh_link::UInt32
            sh_info::UInt32
            sh_addralign::UInt32
            sh_entsize::UInt32
        end

        @struct immutable SymtabEntry <: ELF.ELFSymtabEntry
            st_name::UInt32
            st_value::UInt32
            st_size::UInt32
            st_info::UInt8
            st_other::UInt8
            st_shndx::UInt16
        end

        @struct immutable Rel <: ELF.ELFRel
            r_offset::UInt32
            r_info::UInt32
        end

        @struct immutable Rela <: ELF.ELFRel
            r_offset::UInt32
            r_info::UInt32
            r_addend::Int32
        end

        @struct immutable ProgramHeader <: ELF.ELFProgramHeader
            p_type::UInt32
            p_offset::UInt32
            p_vaddr::UInt32
            p_paddr::UInt32
            p_filesz::UInt32
            p_memsz::UInt32
            p_flags::UInt32
            p_align::UInt32
        end

        type File <: ELF.ELFFile
            endianness::Symbol
            ei_version::UInt8
            ei_osabi::UInt8
            header::Header
        end
        ELF.sheader(::Type{File}) = SectionHeader
        ELF.pheader(::Type{File}) = ProgramHeader
        ELF.symtype(::Type{SectionHeader}) = SymtabEntry
    end

    module ELF64
        import ELF
        using StructIO

        @struct immutable Header <: ELF.ELFHeader
            e_type::UInt16
            e_machine::UInt16
            e_version::UInt32
            e_entry::UInt64
            e_phoff::UInt64
            e_shoff::UInt64
            e_flags::UInt32
            e_ehsize::UInt16
            e_phentsize::UInt16
            e_phnum::UInt16
            e_shentsize::UInt16
            e_shnum::UInt16
            e_shstrndx::UInt16
        end

        @struct immutable SectionHeader <: ELF.ELFSectionHeader
            sh_name::UInt32
            sh_type::UInt32
            sh_flags::UInt64
            sh_addr::UInt64
            sh_offset::UInt64
            sh_size::UInt64
            sh_link::UInt32
            sh_info::UInt32
            sh_addralign::UInt64
            sh_entsize::UInt64
        end

        @struct immutable SymtabEntry <: ELF.ELFSymtabEntry
            st_name::UInt32
            st_info::UInt8
            st_other::UInt8
            st_shndx::UInt16
            st_value::UInt32
            st_size::UInt64
        end

        @struct immutable Rel <: ELF.ELFRel
            r_offset::UInt64
            r_info::UInt64
        end

        @struct immutable Rela <: ELF.ELFRel
            r_offset::UInt64
            r_info::UInt64
            r_addend::Int64
        end

        @struct immutable ProgramHeader <: ELF.ELFProgramHeader
            p_type::UInt32
            p_flags::UInt32
            p_offset::UInt64
            p_vaddr::UInt64
            p_paddr::UInt64
            p_filesz::UInt64
            p_memsz::UInt64
            p_align::UInt64
        end

        type File <: ELF.ELFFile
            endianness::Symbol
            ei_version::UInt8
            ei_osabi::UInt8
            header::Header
        end
        ELF.sheader(::Type{File}) = SectionHeader
        ELF.pheader(::Type{File}) = ProgramHeader
        ELF.symtype(::Type{SectionHeader}) = SymtabEntry
    end

    sectionsize(sh::ELFSectionHeader) = sh.sh_size
    sectionoffset(sh::ELFSectionHeader) = sh.sh_offset
    sectionaddress(sh::ELFSectionHeader) = sh.sh_addr

    # Definitions for ELF Handle
    ELFHandle{T<:IO}(io::T, file::ELFFile) = ELFHandle{T}(io,position(io),file)
    function show(io::IO,h::ELFHandle)
        print(io,"ELF Handle (")
        if typeof(h.file) == ELF32.File
            print(io,"32-bit")
        elseif typeof(h.file) == ELF64.File
            print(io,"64-bit")
        else
            error("Unrecognized ELF type")
        end
        print(io,")")
    end

    intptr(h::ELFHandle) = isa(h.file, ELF32.File) ? UInt32 : UInt64

    for f in (:readuntil,:write)
        @eval $(f){T<:IO}(io::ELFHandle{T},args...) = $(f)(io.io,args...)
    end
    readbytes{T<:IO}(io::ELFHandle{T},num::Integer) = readbytes(io.io,num)
    seek{T<:IO}(io::ELFHandle{T},pos::Integer) = seek(io.io,io.start+pos)
    seekstart(io::ELFHandle) = seek(io.io,io.start)
    position{T<:IO}(io::ELFHandle{T}) = position(io.io)-io.start
    unpack{T,ioT<:IO}(h::ELFHandle{ioT},::Type{T}) = unpack(h.io,T)
    endianness(h::ELFHandle) = h.file.endianness

    function endianness(ei_data::UInt8)
        if ei_data == ELFDATA2MSB
            :BigEndian
        elseif ei_data == ELFDATA2LSB
            :LittleEndian
        else
            error("Invalid Data Specification")
        end
    end

    function readmeta(io::IO,::Type{ELFHandle})
        start = position(io)
        mag0 = read(io,UInt8)
        mag1 = read(io,UInt8)
        mag2 = read(io,UInt8)
        mag3 = read(io,UInt8)
        if((mag0 != UInt32('\177')) || (mag1 != UInt32('E')) || (mag2 != UInt32('L')) || (mag3 != UInt32('F')))
            throw(ObjFileBase.MagicMismatch("Magic Number does not match"))
        end
        class = read(io,UInt8)
        data = read(io,UInt8)
        version = read(io,UInt8)
        osabi = read(io,UInt8)
        abiversion = read(io,UInt8)
        skip(io,7)
        if class == ELFCLASS32
            header = unpack(io,ELF32.Header,endianness(data))
            file = ELF32.File(endianness(data),version,osabi,header)
        elseif class == ELFCLASS64
            header = unpack(io,ELF64.Header,endianness(data))
            file = ELF64.File(endianness(data),version,osabi,header)
        else
            error("Invalid File Class")
        end
        seek(io,start)
        ELFHandle(io,file)
    end
    FileIO.load(s::Stream{format"ELF"}) = readmeta(stream(s))
    function FileIO.save(f::File{format"ELF"},oh::ELFHandle)
        # First we need to calculate the total size of the object file. To this,
        # iterate through all sections and fine which one (if placed last)
        # yields the largest object file
        header = oh.file.header
        size = max(
            maximum(map(sec->(sectionoffset(sec)+sectionsize(sec)),Sections(oh))),
            (header.e_shoff+header.e_shentsize*header.e_shnum))
        open(f,"w") do s
            seekstart(oh)
            write(s,readbytes(oh,size))
        end
    end
    FileIO.save(s::AbstractString,oh::ELFHandle) =
        save(File{format"ELF"}(s),oh)

    function read(io::IO,::Type{ELFProgramHeader},f::ELFFile)
        s = sizeof(pheader(typeof(f)))
        if s > f.header.e_phentsize
            error("Missing data for program header")
        end
        ret = unpack(io,pheader(typeof(f)),f.endianness)
        skip(io,f.header.e_phentsize-s)
        ret
    end

    function read(io::IO,::Type{ELFSectionHeader},f::ELFFile)
        s = sizeof(sheader(typeof(f)))
        if s > f.header.e_shentsize
            error("Missing data for section header")
        end
        ret = unpack(io,sheader(typeof(f)),f.endianness)
        skip(io,f.header.e_shentsize-s)
        ret
    end

    function strtab_lookup(io,strtable::ELFSectionHeader,index)
        seek(io,strtable.sh_offset+index)
        strip(readuntil(io,'\0'),'\0')
    end

    function read(io::IO,x::Array{UInt8,1},file::ELFFile,header::ELFProgramHeader)
        seek(io,header.p_offset)
        read(io,x)
    end

    function read(io::IO,x::Array{UInt8,1},file::ELFFile,header::ELFSectionHeader)
        seek(io,header.sh_offset)
        read(io,x)
    end

    function read(io::IO,file::ELFFile,header::ELFSectionHeader)
        x = Array(UInt8,header.sh_size)
        read(io,x,file,header)
        x
    end

    function read(io::IO,file::ELFFile,header::ELFProgramHeader)
        x = Array(UInt8,header.p_filesz)
        read(io,x,file,header)
        x
    end

    # Program Header iteration
    immutable ProgramHeaders
        handle::ELFHandle
    end
    handle(phs::ProgramHeaders) = phs.handle
    endof(phs::ProgramHeaders) = handle(phs).file.header.e_phnum
    length(phs::ProgramHeaders) = endof(phs)
    start(phs::ProgramHeaders) = 1
    done(phs::ProgramHeaders,n) = n > length(phs)
    next(phs::ProgramHeaders,n) = (phs[n],n+1)
    function getindex(phs::ProgramHeaders, n)
        @assert 0 < n <= length(phs)
        file = handle(phs).file
        seek(handle(phs),file.header.e_phoff + (n-1)*file.header.e_phentsize)
        read(handle(phs),ELFProgramHeader,file)
    end

    function show(io::IO, header::ELFProgramHeader; strtab = nothing, sections = nothing)
        printentry(io,"Type",P_TYPE[header.p_type])
        printentry(io,"Offset","0x",hex(header.p_offset))
        printentry(io,"Virtual Address","0x",hex(header.p_vaddr))
        printentry(io,"Physical Address","0x",hex(header.p_paddr))
        printentry(io,"Size in File","0x",hex(header.p_filesz))
        printentry(io,"Size in Memory","0x",hex(header.p_memsz))
        printentry(io,"Flags","0x",hex(header.p_flags))
        printentry(io,"Align","0x",hex(header.p_align))
    end

    # Access to sections
    function sectionname(header::ELFSectionHeader; strtab = nothing, errstrtab = true)
        if strtab == nothing
            errstrtab && error("No Strtab given")
            return string("strtab@",header.sh_name)
        end
        return strtab_lookup(strtab,header.sh_name)
    end

    sizeof(header::ELFSectionHeader) = header.sh_size

    function secttype(sh_type::UInt32)
        if haskey(SHT_TYPES, sh_type)
            return SHT_TYPES[sh_type]
        end
        return string("Unknown (0x",hex(sh_type),")")
    end

    function show(io::IO, header::ELFSectionHeader; strtab = nothing, sections = nothing)
        printentry(io,"Name",sectionname(header;strtab=strtab,errstrtab=false))
        printentry(io,"Type",secttype(header.sh_type))
        printentry(io,"Size","0x",hex(header.sh_size))
        printentry(io,"Offset","0x",hex(header.sh_offset))
        printentry(io,"Load Address","0x",hex(header.sh_addr))
        if header.sh_link != 0
            target = ""
            if strtab !== nothing && sections !== nothing
                target = string(" -> ",sectionname(sections[header.sh_link+1].header; strtab = strtab))
            end
            printentry(io,"Link Section",header.sh_link,target)
        end
        if header.sh_info != 0
            if header.sh_type == SHT_REL || header.sh_type == SHT_RELA
                target = ""
                if strtab !== nothing && sections !== nothing
                    target = string(" -> ",sectionname(sections[header.sh_info+1].header; strtab = strtab))
                end
                printentry(io,"Info Section",header.sh_info,target)
            elseif header.sh_type == SHT_SYMTAB || header.sh_type == SHT_DYNSYM
                printentry(io,"Last Local Idx",header.sh_info)
            end
        end
        flags = String[]
        for (k,v) in SHF_FLAGS
            ((k&header.sh_flags) != 0) && push!(flags, v)
        end
        !isempty(flags) && printentry(io,"Flags",join(flags,","))
        printentry(io,"Align","0x",hex(header.sh_addralign))
        if header.sh_entsize != 0
            printentry(io,"Entry Size","0x",hex(header.sh_entsize))
        end
    end

    function filetype(e_type)
        if haskey(ET_TYPES, e_type)
            return ET_TYPES[e_type]
        end
        return string("Unknown (0x",hex(e_type),")")
    end
    function machinetype(e_machine)
        if haskey(EM_MACHINES, e_machine)
            return EM_MACHINES[e_machine]
        end
        return string("Unknown (0x",hex(e_machine),")")
    end
    function show(io::IO, header::ELFHeader)
        printentry(io,"Type",filetype(header.e_type))
        printentry(io,"Machine",machinetype(header.e_machine))
        # Skip e_version (not particularly useful)
        printentry(io,"Entrypoint","0x",hex(header.e_entry))
        printentry(io,"PH Offset","0x",hex(header.e_phoff))
        printentry(io,"SH Offset","0x",hex(header.e_shoff))
        # Skip flags
        printentry(io,"Header Size","0x",hex(header.e_ehsize))
        printentry(io,"PH Entry Size","0x",hex(header.e_phentsize))
        printentry(io,"PH Entry Count",dec(header.e_phnum))
        printentry(io,"SH Entry Size","0x",hex(header.e_shentsize))
        printentry(io,"SH Entry Count",dec(header.e_shnum))
        printentry(io,"Strtab Index",dec(header.e_shstrndx))
    end

    ObjFileBase.isrelocatable(handle::ELFHandle) =
        handle.file.header.e_type == ELF.ET_REL

    ObjFileBase.isexecutable(handle::ELFHandle) =
        handle.file.header.e_type == ELF.ET_EXEC

    immutable SectionRef{T<:ELFHandle, hdr} <: ObjFileBase.SectionRef{ELFHandle}
        handle::T
        header::hdr
    end
    @Base.pure ObjFileBase.SectionRef{T<:ELFHandle}(::Type{T}) = SectionRef{T}
    handle(sec::SectionRef) = sec.handle
    sectionname(sec::SectionRef; strtab=load_strtab(sec.handle), errstrtab = true) = sectionname(sec.header; strtab = strtab, errstrtab = true)
    show(io::IO, sr::SectionRef; strtab = load_strtab(sr.handle), sections = nothing) = show(io,sr.header; strtab = strtab, sections = sections)
    sizeof(s::SectionRef) = sizeof(s.header)
    deref(s::SectionRef) = s.header
    seek(s::SectionRef, offs) = seek(s.handle, sectionoffset(s) + offs)
    Base.seekstart(s::SectionRef) = seek(handle(s), sectionoffset(s))
    Base.read(s::SectionRef) = (seek(s,0); read(s.handle.io, sectionsize(s)))

    ObjFileBase.isBSS(x::ELFSectionHeader) = x.sh_type == SHT_NOBITS

    immutable StrTab{T<:SectionRef} <: ObjFileBase.StrTab
        strtab::T
    end
    ObjFileBase.StrTab(strtab::SectionRef) = StrTab(strtab)
    strtab_lookup(s::StrTab,index) = strtab_lookup(s.strtab.handle,s.strtab.header,index)

    immutable Sections{T<:ELFHandle} <: ObjFileBase.Sections{ELFHandle}
        handle::T
    end
    ObjFileBase.handle(sections::Sections) = sections.handle
    ObjFileBase.Sections(handle::ELFHandle) = Sections(handle)
    ObjFileBase.mangle_sname(h::ELFHandle, name) = string(".", name)
    endof(s::Sections) = s.handle.file.header.e_shnum
    length(s::Sections) = endof(s)
    function getindex{T}(s::Sections{T}, n)
        @assert 0 < n <= length(s)
        file = s.handle.file
        seek(s.handle,file.header.e_shoff + (n-1)*file.header.e_shentsize)
        SectionRef(s.handle, read(s.handle,ELFSectionHeader,file))
    end

    load_strtab(h::ELFHandle) = StrTab(Sections(h)[h.file.header.e_shstrndx+1])
    load_strtab{T}(s::Sections{T}) = StrTab(s[s.handle.file.header.e_shstrndx+1])
    load_strtab{T}(s::SectionRef{T}) = StrTab(s)
    const strtab = load_strtab

    start(s::Sections) = 1
    done(s::Sections,n) = n > length(s)
    next(s::Sections,n) = (s[n],n+1)

    function show(io::IO,s::Sections)
        println(io,"ELF Section Table")
        for section in s
            show(io,section; strtab=strtab(s), sections = s)
            println(io)
        end
    end

    link_sec(sec::SectionRef) = Sections(sec.handle)[sec.header.sh_link+1]
    info_sec(sec::SectionRef) = Sections(sec.handle)[sec.header.sh_info+1]

    # # Symbols
    immutable Symbols{T<: SectionRef} <: ObjFileBase.Symbols{ELFHandle}
        symtab::T
    end
    ObjFileBase.Symbols(sec::SectionRef) = Symbols(sec)
    handle(s::Symbols) = handle(s.symtab)
    Symbols(h::ELFHandle) =
        Symbols(first(filter(x->sectionname(x) in (".dynsym",".symtab"),Sections(h))))
    ObjFileBase.Symbols(h::ELFHandle) = Symbols(h)
    ObjFileBase.StrTab(symtab::Symbols) = StrTab(symtab)
    StrTab(symtab::Symbols) = StrTab(link_sec(symtab.symtab))

    immutable SymbolRef <: ObjFileBase.SymbolRef{ELFHandle}
        syms::Symbols
        num::UInt32
        offset::Int
        entry::ELFSymtabEntry
    end
    handle(sym::SymbolRef) = handle(sym.syms)
    symname(sym::SymbolRef; strtab = StrTab(sym.syms), kwargs...) = symname(sym.entry; strtab=strtab, kwargs...)
    deref(ref::SymbolRef) = ref.entry
    symbolnum(ref::SymbolRef) = ref.num

    function symname(sym::ELFSymtabEntry; strtab = nothing, errstrtab = true)
        if strtab == nothing
            errstrtab && error("No Strtab given")
            return string("strtab@",sym.st_name)
        end
        return strtab_lookup(strtab,sym.st_name)
    end

    st_bind(st_info) = st_info>>4
    st_type(st_info::Integer) = st_info & 0xf
    st_type(x::SymbolRef) = st_type(deref(x))
    st_type(x::ELFSymtabEntry) = st_type(x.st_info)
    isglobal(x) = (st_bind(x.st_info) & STB_GLOBAL) != 0
    islocal(x) = !isglobal(x)
    isweak(x) = (st_bind(x.st_info) & STB_WEAK) != 0
    isdebug(x) = false
    isundef(x::ELFSymtabEntry) = x.st_shndx == SHN_UNDEF
    isundef(x::SymbolRef) = isundef(deref(x))

    function compute_first_ph_vaddr(h)
        for ph in ProgramHeaders(h)
            ph.p_type != PT_LOAD && continue
            return ph.p_vaddr
        end
        return 0
    end

    function symbolvalue(sym::Union{SymbolRef, ELFSymtabEntry}, sections,
            first_ph_vaddr=compute_first_ph_vaddr(handle(sections)))
        value = deref(sym).st_value
        shndx = deref(sym).st_shndx
        if shndx != ELF.SHN_UNDEF && shndx < ELF.SHN_LORESERVE
            sec = sections[shndx+1]
            # What to do here depends on the object kind. Shared Libraries and
            # executable's st_value's are virtual addresses
            if deref(sec).sh_addr != 0 && handle(sec).file.header.e_type == ET_REL
                value += deref(sec).sh_addr
            end
            if handle(sec).file.header.e_type == ET_EXEC ||
                 handle(sec).file.header.e_type == ET_DYN
                value -= first_ph_vaddr
            end
        end
        value
    end

    # Symbol printing stuff
    function showcompact(io::IO, x::SymbolRef; shstrtab = load_strtab(handle(x)), strtab = StrTab(x.syms), sections = Sections(handle(x)))
        print(io,'[')
        printfield(io,dec(symbolnum(x)),5)
        print(io,"] ")
        showcompact(io, x.entry; shstrtab = shstrtab, strtab = strtab, sections = sections)
    end
    show(io::IO, x::SymbolRef) = showcompact(io, x)

    # Try to follow the same format as llvm-objdump
    function showcompact(io::IO,x::ELFSymtabEntry; shstrtab = nothing, strtab = nothing, sections = nothing)
        # Value
        print(io,string("0x",hex(x.st_value,2*sizeof(x.st_value))))
        print(io," ")

        # Size
        print(io,string("0x",hex(x.st_size,2*sizeof(x.st_size))))
        print(io," ")

        # Symbol flags
        print(io, isglobal(x) ? "g" : islocal(x) ? "l" : "-")
        print(io, isweak(x) ? "w" : "-")
        print(io, "-"^3) # Unsupported
        print(io, isdebug(x) ? "d" : "-")
        STT = st_type(x)
        # Symbol type
        print(io, STT == STT_FILE ? "F" : STT == STT_FUNC ? "f" : STT == STT_OBJECT ? "O" : "-")

        print(io, " ")
        if x.st_shndx == SHN_UNDEF
            printfield(io,"*UND*",20; align = :left)
        elseif x.st_shndx == SHN_COMMON
            printfield(io,"*COM*",20; align = :left)
        elseif x.st_shndx == SHN_ABS
            printfield(io,"*ABS*",20; align = :left)
        elseif sections !== nothing
            printfield(io, sectionname(sections[x.st_shndx+1];
                strtab = shstrtab, errstrtab=true), 20; align = :left)
        else
            printfield(io, "Section #$(x.st_shndx)", 20; align = :left)
        end
        print(io, " ")

        print(io,symname(x; strtab = strtab, errstrtab = false))
    end

    function show(io::IO, s::Symbols)
        h = s.symtab.handle
        shstrtab = strtab(h)
        symstrtab = strtab(Sections(h)[s.symtab.header.sh_link+1])
        for sym in s
            showcompact(io, sym; shstrtab = shstrtab, strtab = symstrtab)
            println(io)
        end
    end

    SymtabEntrySize(s::Symbols) = sizeof(symtype(typeof(s.symtab.header)))
    endof(s::Symbols) = div(s.symtab.header.sh_size,SymtabEntrySize(s))
    function getindex(s::Symbols,n)
        if n < 1 || n > endof(s)
            throw(BoundsError())
        end
        h = s.symtab.handle
        offset = s.symtab.header.sh_offset + (n-1)*SymtabEntrySize(s)
        seek(h,offset)
        SymbolRef(s,n,offset,unpack(h, symtype(typeof(s.symtab.header))))
    end

    start(s::Symbols) = 1
    done(s::Symbols,n) = n > endof(s)
    next(s::Symbols,n) = (x=s[n];(x,n+1))
    length(s::Symbols) = endof(s)

    # Access to relocations
    immutable Relocations{T <: ELFRel, S <: SectionRef}
        sec::S
    end
    function Relocations(sec::SectionRef)
        is64 = isa(sec.handle.file,ELF64.File)
        isRela = sec.header.sh_type == SHT_RELA
        Relocations{is64 ? (isRela ? ELF64.Rela : ELF64.Rel) : (isRela ? ELF32.Rela : ELF32.Rel),typeof(sec)}(sec)
    end

    immutable RelocationRef{T <: ELFRel} <: ObjFileBase.RelocationRef{ELFHandle}
        h::ELFHandle
        reloc::T
    end

    deref(x::RelocationRef) = x.reloc

    entrysize{T}(s::Relocations{T}) = sizeof(T)
    endof{T}(s::Relocations{T}) = div(s.sec.header.sh_size,entrysize(s))
    length(r::Relocations) = endof(r)
    function getindex{T}(s::Relocations{T},n)
        if n < 1 || n > length(s)
            throw(BoundsError())
        end
        offset = sectionoffset(s.sec) + (n-1)*entrysize(s)
        seek(s.sec.handle,offset)
        RelocationRef{T}(s.sec.handle,unpack(s.sec.handle, T))
    end

    start(s::Relocations) = 1
    done(s::Relocations,n) = n > length(s)
    next(s::Relocations,n) = (x=s[n];(x,n+1))

    # DWARF support
    #=
    function read(io::IO,file::ELFFile,h::ELFSectionHeader,::Type{DWARF.ARTable})
        seek(io,h.sh_offset)
        ret = DWARF.ARTable(Array(DWARF.ARTableSet,0))
        while position(io) < h.sh_offset + h.sh_size
            push!(ret.sets,read(io,DWARF.ARTableSet,f.endianness))
        end
        ret
    end

    function read(io::IO,file::ELFFile,h::ELFSectionHeader,::Type{DWARF.PUBTable})
        seek(io,h.sh_offset)
        ret = DWARF.PUBTable(Array(DWARF.PUBTableSet,0))
        while position(io) < h.sh_offset + h.sh_size
            push!(ret.sets,read(io,DWARF.PUBTableSet,f.endianness))
        end
        ret
    end

    function read(io::IO,f::ELFFile,h::ELFSectionHeader,::Type{DWARF.AbbrevTableSet})
        seek(io,h.sh_offset)
        read(io,AbbrevTableSet,f.endianness)
    end

    function read(io::IO,f::ELFFile,h::ELFSectionHeader,s::DWARF.PUBTableSet,::Type{DWARF.DWARFCUHeader})
        seek(io,h.sh_offset+s.header.debug_info_offset)
        read(io,DWARF.DWARFCUHeader,f.endianness)
    end

    function read(io::IO,f::ELFFile,debug_info::ELFSectionHeader,debug_abbrev::ELFSectionHeader,
        s::DWARF.PUBTableSet,e::DWARF.PUBTableEntry,header::DWARF.DWARFCUHeader,::Type{DWARF.DIE})
        ats = read(io,f,debug_abbrev,header,DWARF.AbbrevTableSet)
        seek(io,debug_info.sh_offset+s.header.debug_info_offset+e.offset)
        read(io,header,ats,DWARF.DIE)
    end

    function read(io::IO,f::ELFFile,h::ELFSectionHeader,s::DWARF.DWARFCUHeader,::Type{DWARF.AbbrevTableSet})
        seek(io,h.sh_offset+s.debug_abbrev_offset)
        read(io,AbbrevTableSet,f.endianness)
    end

    function debugsections(io::IO,f::ELFFile)
        snames = names(io,f,f.sheaders)
        sections = Dict{String,ELFSectionHeader}()
        for i in 1:length(snames)
            # Remove leading "."
            ind = findfirst(DEBUG_SECTIONS,snames[i][2:end])
            if ind != 0
                sections[DEBUG_SECTIONS[ind]] = f.sheaders[ind]
            end
        end
        sections
    end

    function read(io::IO,f::ELFFile,debug_info::ELFSectionHeader,debug_abbrev::ELFSectionHeader,
        s::DWARF.PUBTableSet,e::DWARF.PUBTableEntry,header::DWARF.DWARFCUHeader,::Type{DWARF.DIETree})
        ats = read(io,f,debug_abbrev,header,DWARF.AbbrevTableSet)
        seek(io,debug_info.sh_offset+s.header.debug_info_offset+e.offset)
        ret = DIETree(Array(DWARF.DIETreeNode,0))
        read(io,header,ats,ret,DWARF.DIETreeNode,f.endianness)
        ret
    end
    =#

    immutable dl_phdr_info
        dlpi_addr::UInt64
        dlpi_name::Ptr{UInt8}
        dlpi_phdr::Ptr{Void}
        dlpi_phnum::UInt16
    end

    function callback(info::Ptr{dl_phdr_info},size::Csize_t, data::Ptr{Void})
        push!(unsafe_pointer_to_objref(data),unsafe_load(info))
        convert(Cint,0)
    end

    function loaded_libraries()
        x = Array(dl_phdr_info,0)
        ccall(:dl_iterate_phdr, Cint, (Ptr{Void}, Any), cfunction(callback, Cint, (Ptr{dl_phdr_info},Csize_t,Ptr{Void})), x)
        x
    end

    ## DWARF Support
    function debugsections(h::ELFHandle)
        sects = collect(Sections(h))
        strt = strtab(h)
        snames = map(s->sectionname(s.header;strtab=strt),sects)
        sections = Dict{String,SectionRef}()
        for i in 1:length(snames)
            # remove leading "."
            ind = findfirst(ObjFileBase.DEBUG_SECTIONS,string(snames[i])[2:end])
            if ind != 0
                sections[ObjFileBase.DEBUG_SECTIONS[ind]] = sects[i]
            end
        end
        ObjFileBase.DebugSections(h,sections)
    end

    # JIT Utils
    function replace_sections_from_memory(h::ELFHandle, new_buffer)
        for sec in Sections(h)
            if ObjFileBase.is_jit_section(sec)
                seek(new_buffer,sectionoffset(sec))
                write(new_buffer,pointer_to_array(
                    reinterpret(Ptr{UInt8},sectionaddress(sec)),
                    sectionsize(sec),false))
            end
        end
        seekstart(new_buffer)
        new_buffer
    end

    # Other things
    include("relocate.jl")
end


text:
using JuMP

export Layer, NeuralNet

JuMPLinearType = Union{JuMP.VariableRef,JuMP.AffExpr}
JuMPReal = Union{Real,JuMPLinearType}

include("net_components/core_ops.jl")

"""
$(TYPEDEF)

Supertype for all types storing the parameters of each layer. Inherit from this
to specify your own custom type of layer. Each implementation is expected to:

1. Implement a callable specifying the output when any input of type `JuMPReal` is provided.
"""
abstract type Layer end

"""
An array of `Layers` is interpreted as that array of layer being applied
to the input sequentially, starting from the leftmost layer. (In functional programming
terms, this can be thought of as a sort of `fold`).
"""
chain(x::Array{<:JuMPReal}, ps::Array{<:Layer,1}) = length(ps) == 0 ? x : chain(ps[1](x), ps[2:end])

(ps::Array{<:Layer,1})(x::Array{<:JuMPReal}) = chain(x, ps)

function check_size(input::AbstractArray, expected_size::NTuple{N,Int})::Nothing where {N}
    input_size = size(input)
    @assert input_size == expected_size "Input size $input_size did not match expected size $expected_size."
end

include("net_components/layers.jl")

"""
$(TYPEDEF)

Supertype for all types storing the parameters of a neural net. Inherit from this
to specify your own custom architecture. Each implementation
is expected to:

1. Implement a callable specifying the output when any input of type `JuMPReal` is provided
2. Have a `UUID` field for the name of the neural network.
"""
abstract type NeuralNet end

include("net_components/nets.jl")


text:
module BenchCollect

using BenchmarkTools
using Transducers

suite = BenchmarkGroup()

function collect_unordered(xf, xs; basesize = 1, kwargs...)
    T = eltype(xs)
    return collect(
        T,
        channel_unordered(
            xf,
            xs;
            eltype = T,
            size = basesize,
            basesize = basesize,
            kwargs...,
        ),
    )
end

xf = Map() do x
    for _ in 1:1000
        x = sin(x)
    end
    return x
end

xs = randn(2^15)

let
    ys = collect(xf, xs) |> sort!
    ys_unordered_1 = collect_unordered(xf, xs; basesize = 1) |> sort!
    ys_unordered_32 = collect_unordered(xf, xs; basesize = 32) |> sort!
    ys_unordered_1024 = collect_unordered(xf, xs; basesize = 1024) |> sort!
    @assert ys == ys_unordered_1
    @assert ys == ys_unordered_32
    @assert ys == ys_unordered_1024
end

suite["seq"] = @benchmarkable collect($xf, $xs)
suite_assoc = suite["assoc"] = BenchmarkGroup()
for basesize in [1, 2^5, 2^10]
    suite_assoc["basesize=$basesize"] = @benchmarkable tcollect(
        $xf,
        $xs;
        basesize = $basesize,
    )
end
suite_unordered = suite["unordered"] = BenchmarkGroup()
for basesize in [1, 2^5, 2^10]
    suite_unordered["basesize=$basesize"] =
        @benchmarkable collect_unordered($xf, $xs; basesize = $basesize)
end

end  # module
BenchCollect.suite


text:
"""

E-step of the EM-algorithm.

	filtersmoother(model,data,w)

Calibrates a w[x,y,x',y'] array of weights such that the EM-step is obtained by maximizing
sum_{x,y,x',y'} w[x,y,x',y']\\*log(m[x,y,x',y']) which in Julia is simply `sum(w.*m)`.
The weights are obtained by running the filter/smoother (aka forward/backward) algorithm.
The `EMalgorithm` function wraps `filtersmoother` along with a generic numerical optimization
step, but more efficient (even closed-form) M-step optimization may be available in particular
models, such as Hidden Markov Models: if this the case call `filtersmoother` to build the weights and
implement your own optimization routine.

"""
function estep(model::DynamicDiscreteModel,data::Array{Int,1},w::Array{Float64,4})
	T=length(data)
	dx,dy=size(model.mu)

	filter=Array(Float64,dx,T)
	filter[:,1]=model.mu[:,data[1]]/sum(model.mu[:,data[1]])
	for t=2:T
		rho=0.0
		for jx=1:dx
			filter[jx,t]=0
			for ix=1:dx
				filter[jx,t]+=model.m[ix,data[t-1],jx,data[t]]*filter[ix,t-1]
			end
			rho+=filter[jx,t]
		end
		for jx=1:dx
			filter[jx,t]/=rho
		end
	end

	smoother=Array(Float64,dx,T)
	smoother[:,T]=1
	for t=T-1:-1:1
		rho=0.0
		for ix=1:dx
			smoother[ix,t]=0
			for jx=1:dx
				smoother[ix,t]+=model.m[ix,data[t],jx,data[t+1]]*smoother[jx,t+1]
			end
			rho+=smoother[ix,t]
		end
		for ix=1:dx
			smoother[ix,t]/=rho
		end
	end

	conditional=Array(Float64,dx,dx)
	for t=1:T-1
		tempsum=0.0
		for jx=1:dx
			for ix=1:dx
				conditional[ix,jx]=model.m[ix,data[t],jx,data[t+1]]
				conditional[ix,jx]*=filter[ix,t]
				conditional[ix,jx]*=smoother[jx,t+1]
				tempsum+=conditional[ix,jx]
			end
		end
		for jx=1:dx
			for ix=1:dx
				#note that there is no risk of division by zero
				#reason: this is called only on _observed_ data which must thus have non-zero probability
				w[ix,data[t],jx,data[t+1]]+=conditional[ix,jx]/tempsum
			end
		end
	end

end

# Wrapper to handle several time-series of observations.
function estep(model::DynamicDiscreteModel,data::Array{Array,1},w::Array{Float64,4})
	for i=1:length(data)
		estep(model,data[i],w)
	end
end


text:
using Documenter
using OrbitalTrajectories

makedocs(
    sitename = "OrbitalTrajectories",
    format = Documenter.HTML(),
    modules = [OrbitalTrajectories],
    pages = [
        "OrbitalTrajectories.jl" => "index.md",
        "Examples" => [
            "tutorials/DynamicModels.md",
            "tutorials/SunMarsCR3BP.md",
            "tutorials/STM.md",
            "tutorials/QSOFamilies.md",
            "tutorials/Koopman.md"
        ],
        "Package API" => "API.md"
    ],
    strict = true
)

REPO = get(ENV, "GITHUB_REPOSITORY", nothing)
@info "REPO: $(REPO)"
if !isnothing(REPO)
    deploydocs(
        repo = "github.com/$(REPO)",
        forcepush = true
    )
end


text:
# ---
# title: 48. Rotate Image
# id: problem48
# author: Indigo
# date: 2022-03-03
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/rotate-image/description/>
# hidden: true
# ---
#
# You are given an _n_ x _n_ 2D `matrix` representing an image, rotate the image
# by 90 degrees (clockwise).
#
# You have to rotate the image [**in-place**](https://en.wikipedia.org/wiki/In-
# place_algorithm), which means you have to modify the input 2D matrix directly.
# **DO NOT** allocate another 2D matrix and do the rotation.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/28/mat1.jpg)
#
#
#
#     Input: matrix = [[1,2,3],[4,5,6],[7,8,9]]
#     Output: [[7,4,1],[8,5,2],[9,6,3]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/28/mat2.jpg)
#
#
#
#     Input: matrix = [[5,1,9,11],[2,4,8,10],[13,3,6,7],[15,14,12,16]]
#     Output: [[15,13,2,5],[14,3,4,1],[12,6,8,9],[16,7,10,11]]
#
#
# **Example 3:**
#
#
#
#     Input: matrix = [[1]]
#     Output: [[1]]
#
#
# **Example 4:**
#
#
#
#     Input: matrix = [[1,2],[3,4]]
#     Output: [[3,1],[4,2]]
#
#
#
#
# **Constraints:**
#
#   * `matrix.length == n`
#   * `matrix[i].length == n`
#   * `1 <= n <= 20`
#   * `-1000 <= matrix[i][j] <= 1000`
#
#
## @lc code=start
using LeetCode

function rotate_48(matrix::Matrix{Int})
    n = size(matrix, 1)
    for i in 1:div(n, 2), j in i:(n - i)
        matrix[i, j], matrix[j, n + 1 - i], matrix[n + 1 - i, n + 1 - j], matrix[n + 1 - j, i] = matrix[
            n + 1 - j, i
        ],
        matrix[i, j], matrix[j, n + 1 - i],
        matrix[n + 1 - i, n + 1 - j]
    end
    return matrix
end
## @lc code=end


text:
module LambdaRuntime

using HTTP
using JSON

function post_error(path::AbstractString, ex::Exception)
    headers = [
        "Content-type" => "application/json",
        "Lambda-Runtime-Function-Error-Type" => "Unhandled",
    ]
    body = JSON.json(Dict(
        "errorType" => string(typeof(ex)),
        "errorMessage" => string(ex)
    ))
    local url
    try
        url = "$RUNTIME_URL/$path"
        @async HTTP.post(url, headers, JSON.json(body))
        @info "Notified lambda runtime about the error" url ex
    catch failure
        @error "Unable to notify lambda runtime about the error" url ex failure
    end
end

function init(mod::Module)
    try
        global AWS_LAMBDA_RUNTIME_API = get(ENV, "AWS_LAMBDA_RUNTIME_API", "")
        global LAMBDA_TASK_ROOT = get(ENV, "LAMBDA_TASK_ROOT", "")
        global HANDLER = Symbol(get(ENV, "_HANDLER", "handle_event"))
        global RUNTIME_URL = "http://$AWS_LAMBDA_RUNTIME_API/2018-06-01"
        @info "Settings" AWS_LAMBDA_RUNTIME_API LAMBDA_TASK_ROOT HANDLER
    catch ex
        @error "Initializtion error" ex
        post_error("runtime/init/error", ex)
        @info "Shutting down container"
        exit(1)
    end
end

function main_loop(mod::Module)
    @info "Start processing requests"
    while true
        local state, request, request_id
        try
            state = :start
            request = HTTP.get("$RUNTIME_URL/runtime/invocation/next")
            request_id = Dict(request.headers)["Lambda-Runtime-Aws-Request-Id"]
            @info "Received event" request_id
            state = :received

            response = @eval $(mod).$(HANDLER)(String($request.body), $request.headers)
            @info "Got response from handler" response
            state = :handled

            @async HTTP.post("$RUNTIME_URL/runtime/invocation/$request_id/response", [], response)
            @info "notified lambda runtime"
            state = :finished
        catch ex
            if state == :start
                @error "Unable to receive request" ex
            elseif state == :received
                @error "Unable to handle event" ex request
                post_error("runtime/invocation/$request_id/error", ex)
            elseif state == :handled
                @error "Unable to notify lambda runtime (invocation response)" ex request
            else
                @error "Unknown state" ex state
            end
            sleep(0.2)
        end
    end
    @info "Shutdown complete"
    exit(0)
end

end


text:
include(joinpath(pwd(), "src/direct_policy_optimization/dpo.jl"))
include(joinpath(@__DIR__, "car_obstacles.jl"))

# Additive noise model
model = Car(3, 2, 3)

function fd(model::Car, x⁺, x, u, w, h, t)
    midpoint_implicit(model, x⁺, x, u, w, h) - w
end

# Nominal solution
X̄, Ū = unpack(Z̄, prob)
prob_nom = prob.prob

# DPO
N = 2 * model.n
D = 2 * model.d

β = 1.0
δ = 1.0e-3

# initial samples
x1_sample = resample(x1, Diagonal([1.0; 1.0; 0.1]), 1.0e-1)

# mean problem
prob_mean = trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				h = h,
				dynamics = false
				)

# sample problems
prob_sample = [trajectory_optimization(
				model,
				EmptyObjective(),
				T,
				h = h,
				xl = state_bounds(model, T, x1 = x1_sample[i])[1],
				xu = state_bounds(model, T, x1 = x1_sample[i])[2],
				ul = ul,
				uu = uu,
				dynamics = false,
				con = con_obstacles
				) for i = 1:N]

# sample objective
Q = [(t < T ? Diagonal([10.0; 10.0; 1.0])
	: Diagonal(100.0 * ones(model.n))) for t = 1:T]
R = [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1]

obj_sample = sample_objective(Q, R)
policy = linear_feedback(model.n, model.m)
dist = disturbances([Diagonal(δ * [1.0; 1.0; 0.1]) for t = 1:T-1])
sample = sample_params(β, T)

prob_dpo = dpo_problem(
	prob_nom, prob_mean, prob_sample,
	obj_sample,
	policy,
	dist,
	sample)

# TVLQR policy
K = tvlqr(model, X̄, Ū, Q, R, h)

z0_dpo = zeros(prob_dpo.num_var)
z0_dpo[prob_dpo.prob.idx.nom] = pack(X̄, Ū, prob_nom)
z0_dpo[prob_dpo.prob.idx.mean] = pack(X̄, Ū, prob_nom)
for i = 1:N
	z0_dpo[prob_dpo.prob.idx.sample[i]] = pack(X̄, Ū, prob_nom)
end
for j = 1:(N + D)
	z0_dpo[prob_dpo.prob.idx.slack[j]] = vcat(X̄[2:end]...)
end
for t = 1:T-1
	z0_dpo[prob_dpo.prob.idx.policy[prob_dpo.prob.idx.θ[t]]] = vec(copy(K[t]))
end

# Solve
Z = solve(prob_dpo, copy(z0_dpo),
	tol = 1.0e-3, c_tol = 1.0e-3,
	max_iter = 1000)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using ITensors
using LinearAlgebra
using NDTensors

"""
Tensor data structure for representing tensors and keeping track of hyper indices.
"""

export QXTensor, Index, hyperindices, contract_tensors, tensor_data

"""Datastructure representing tensors"""
mutable struct QXTensor
    rank::Int64
    indices::Array{<:Index}
    hyper_indices::Array{Array{Int64, 1}, 1}
    storage::AbstractArray
end

"""Custom show for QXTensors"""
function Base.show(io::IO, ::MIME"text/plain", t::QXTensor)
    print(io, "QXTensor, rank: $(t.rank), " *
              "dims: $(Tuple(dim.(t.indices))), " *
              "storage: $(typeof(t.storage)), " *
              "hyper_indices: $(t.hyper_indices)")
end

"""Custom size function"""
Base.size(a::QXTensor) = Tuple(dim.(a.indices))
"""Implement inds for QXTensor"""
ITensors.inds(a::QXTensor) = a.indices
"""Implement store for QXTensor"""
ITensors.store(a::QXTensor) = a.storage

"""
    QXTensor(a::T) where T <: Number

QXTensor constructor which creates a new instance of QXTensor corresponding to a scalar
"""
function QXTensor(a::Number)
    QXTensor(0, Array{Index, 1}(), Array{Array{Index, 1}, 1}(), fill(a, ()))
end

"""
    QXTensor(indices::Vector{<:Index},
             hyper_indices::Union{Nothing, Vector{<:Vector{Int64}}},
             storage::Union{Nothing, <: AbstractArray}=nothing;
             diagonal_check::Bool=true)

QXTensor constructor creates a new instance of QXTensor with the given indices
and hyper indices. If no storage data structure is given then a MockTensor of that shape
is added as the storage. If diagonal_check is true, it will automaticallly check which indices are hyper indices
and record in the hyper_indices field. If hyper_indices are given, then these are used.
"""
function QXTensor(indices::Vector{<:Index},
                  hyper_indices::Union{Nothing, Vector{<:Vector{Int64}}}=nothing,
                  storage::Union{Nothing, <: AbstractArray}=nothing;
                  diagonal_check::Bool=true)
    if storage === nothing
        storage = MockTensor(Tuple(dim.(indices)))
    end

    # if hyper indices not given, diagonal check enabled and not mock tensor,
    # then attempt to detect hyperindices directly from the tensor
    if !(storage isa MockTensor)
        if diagonal_check && hyper_indices === nothing
            hyper_indices = find_hyper_edges(storage)
        end
        # if it is not already in reduced form then we reduce it
        if length(hyper_indices) > 0 && ndims(storage) == length(indices)
            storage = reduce_tensor(storage, hyper_indices)
        end
    end

    if hyper_indices === nothing
        hyper_indices = Vector{Vector{Int64}}()
    end

    QXTensor(length(indices), indices, hyper_indices, storage)
end

QXTensor(i::Vector{<:Index}, t::Union{Nothing, <: AbstractArray}; kwargs...) =
         QXTensor(i, nothing, t; kwargs...)
QXTensor(t::Union{Nothing, <: AbstractArray}; kwargs...) =
         QXTensor([Index(x) for x in size(t)], nothing, t; kwargs...)

"""
    tensor_data(tensor::QXTensor; consider_hyperindices::Bool=true)

Get the data associated with the given tensor. If the consider_hyperindices flag is true
then the rank is reduced to merge related indices. For example for a 5 rank tensor where
the 2nd and 4th indices form a group of hyper indices, with this option set to true would
return a rank 4 tensor where the 2nd index has been merged with the 4th. With `consider_hyperindices` set to false a rank 5 tensor
is returned.
"""
function tensor_data(tensor::QXTensor; consider_hyperindices::Bool=true)
    data = store(tensor)
    hi = hyperindices(tensor)
    if consider_hyperindices || length(hi) == 0
        # data is stored in reduced form
        return data
    else
        hi_ranks = indices2ranks(tensor, hyperindices(tensor))
        return collect(expand_tensor(data, hi_ranks))
    end
end

"""
    indices2ranks(tensor::QXTensor, hi::Vector{<:Vector{<:Index}})

Function convert groups of indices to groups of index positions (ranks).
"""
function indices2ranks(tensor::QXTensor, hi::Vector{<:Vector{<:Index}})
    # create an array of the ranks from the groups of hyper indices
    hi_ranks = Array{Int64, 1}[]
    all_indices = inds(tensor)
    for group in hi
        push!(hi_ranks, sort(map(x -> findfirst(y -> y == x, all_indices) ,group)))
    end
    hi_ranks
end

"""
    hyperindices(t::QXTensor)

Function to get the hyper indices as an array of Indices. If the
all_indices flag is true, then all indices are returned, if false
then just the groups of 2 or more are returned.
"""
function hyperindices(t::QXTensor; all_indices=false)
    indices = copy(t.indices)
    hyper_indices = map(x -> Index[x], indices)
    for group in t.hyper_indices
        group = sort(group)
        append!(hyper_indices[group[1]], indices[group[2:end]])
        empty!.(hyper_indices[group[2:end]])
    end
    filter(x -> length(x) > (all_indices ? 0 : 1), hyper_indices)
end

"""
    disable_hyperindices!(t::QXTensor)

Function to disable use of hyper indices with this tensor by removing the hyper
indices and reshaping storage
"""
function disable_hyperindices!(t::QXTensor)
    if length(t.hyper_indices) > 0
        t.storage = collect(expand_tensor(t.storage, t.hyper_indices))
        empty!(t.hyper_indices)
    end
    nothing
end


text:
#############################
### Google routing module ###
#############################

"""
Dictionary for Google Distances API requests:

**Keys**
* `:url` : url for google API, only JSON files outputs are accepted
* `:mode` : transportation mode used in simulation, in the current library scope only driving is accepted
* `:avoid` : map features to avoid (to mantain compatibility with OSM routes ferries should be avoided)
* `:units` : unit system for displaing distances (changing to *imperial* needs deeper changes in both OSMsim and OpenStreetMapX modules)

"""
const googleAPI_parameters = Dict{Symbol,String}(
:url => "https://maps.googleapis.com/maps/api/directions/json?",
:mode  => "driving",
:avoid => "ferries",
:units => "metric",
)

"""
    node_to_string(node_id::Int,map_data::MapData)

Convert node coordinates (stored in ENU system in the `nodes` field of `map_data`)
identified by `node_id` to string with LLA system coordinates

"""
function node_to_string(node_id::Int,map_data::MapData)
    coords = LLA(map_data.nodes[node_id],map_data.bounds)
    return string(coords.lat,",",coords.lon)
end

"""
Gets Google Distances API request url with three points (origin, destination, waypoint between)

**Arguments**
* `origin` : unique node id
* `destination` : unique node id
* `waypoint` : unique node id
* `map_data` : `OpenStreetMapX.MapData` object
* `googleapi_key`: Goole API key
* `googleapi_parameters` : dictionary with assumptions about Google Distances API request

"""
function get_googleapi_url(origin::Int,destination::Int, waypoint::Int,
                            map_data::OpenStreetMapX.MapData, googleapi_key::String;
                            googleapi_parameters::Dict{Symbol,String} = OpenStreetMapX.googleAPI_parameters)
    origin = OpenStreetMapX.node_to_string(origin, map_data)
    destination = OpenStreetMapX.node_to_string(destination, map_data)
    waypoint = OpenStreetMapX.node_to_string(waypoint, map_data)
    return googleapi_parameters[:url]*"origin="*origin*"&destination="*destination*"&waypoints="*waypoint*
    "&avoid="*googleapi_parameters[:avoid]*"&units="*googleapi_parameters[:units]*
    "&mode="*googleapi_parameters[:mode]*"&key="*googleapi_key
end

"""
Gets Google Distances API request url with two points (origin and destination)

**Arguments**
* `origin` : unique node id
* `destination` : unique node id
* `map_data` : `OpenStreetMapX.MapData;` object
* `googleapi_key`: Google API key
* `googleapi_parameters` : dictionary with assumptions about Google Distances API request

"""
function get_googleapi_url(origin::Int,destination::Int,
                            map_data::OpenStreetMapX.MapData, googleapi_key::String;
                            googleapi_parameters::Dict{Symbol,String} = OpenStreetMapX.googleAPI_parameters)
    origin = OpenStreetMapX.node_to_string(origin, map_data)
    destination = OpenStreetMapX.node_to_string(destination, map_data)
    return googleapi_parameters[:url]*"origin="*origin*"&destination="*destination*
    "&avoid="*googleapi_parameters[:avoid]*"&units="*googleapi_parameters[:units]*
    "&mode="*googleapi_parameters[:mode]*"&key="*googleapi_key
end

"""
Get JSON file from Google Distances API request and extract results

**Arguments**
* `url` : string with proper url

"""
function parse_google_url(url::String)
    status, routes = nothing, nothing
    res_json = JSON.parse(String(HTTP.get(url).body))
    status, routes = res_json["status"], res_json["routes"]
    return status, routes
end

"""
Extract route from Google API results

**Arguments**
* `routes` : dictionary with informations about the route

"""
function extract_google_route(routes::Dict)
    res = Array{Tuple{Float64,Float64},1}[]
    legs = routes["legs"]
    for leg in legs
        steps = leg["steps"]
        for step in steps
            push!(res,OpenStreetMapX.decode(step["polyline"]["points"]))
        end
    end
    return vcat(res...)
end

"""
Match Google route with vertices of map network

**Arguments**
* `route` : array with LLA coordinates of crucial route points
* `map_data` : `OpenStreetMapX.MapData;` object

"""
function google_route_to_network(route::Array{Tuple{Float64,Float64},1},map_data::MapData)
    route = [ENU(LLA(coords[1], coords[2]),map_data.bounds) for coords in route]
    res = [nearest_node(map_data, route[1])]
    index = 2
    for i = 2:length(route)
        node = OpenStreetMapX.nearest_node(map_data, route[i])
        if node != res[index-1]
            push!(res,node)
            index += 1
        end
    end
    return res
end

"""
    get_google_route(origin::Int, destination::Int, waypoint::Int,
                     map_data:MapData, googleapi_key::String;
                     googleapi_parameters::Dict{Symbol,String} = googleAPI_parameters)

Get route from to based on Google Distances API with three points
(`origin`, `destination` and `waypoint` between)
on map `map_data` using Google API key `googleapi_key` with optional
Google Distances API request parameters `googleapi_parameters`.

"""
function get_google_route(origin::Int,destination::Int,waypoint::Int,
                            map_data::OpenStreetMapX.OpenStreetMapX.MapData, googleapi_key::String;
                            googleapi_parameters::Dict{Symbol,String} = googleAPI_parameters)
    url = OpenStreetMapX.get_googleapi_url(origin, destination, waypoint,map_data,googleapi_key,googleapi_parameters = googleapi_parameters)
    status, routes = OpenStreetMapX.parse_google_url(url)
    if status == "OK"
        route = OpenStreetMapX.extract_google_route(routes[1])
        return OpenStreetMapX.google_route_to_network(route,map_data), "google"
    elseif status =="ZERO_RESULTS"
        return Int[], "google"
    elseif status == "OVER_QUERY_LIMIT"
        sleep(0.5)
        return OpenStreetMapX.get_google_route(origin,destination,waypoint,map_data,googleapi_key,googleapi_parameters = googleapi_parameters)
    else
        #get route based on OSM routing
        @warn "Google Distances API cannot get a proper results - route will be calculated with OSMSim Routing module"
		if rand() < 0.5
			route_nodes, distance, route_time = OpenStreetMapX.shortest_route(map_data, origin, waypoint, destination)
			return route_nodes, "shortest"
		else
			route_nodes, distance, route_time = OpenStreetMapX.fastest_route(map_data, origin, waypoint, destination)
			return route_nodes, "fastest"
		end
    end
end

"""
    get_google_route(origin::Int, destination::Int,
                     map_data:MapData, googleapi_key::String;
                     googleapi_parameters::Dict{Symbol,String} = googleAPI_parameters)

Get route from to based on Google Distances API with two points (`origin` and `destination`)
on map `map_data` using Google API key `googleapi_key` with optional
Google Distances API request parameters `googleapi_parameters`.

"""
function get_google_route(origin::Int,destination::Int,
                            map_data::MapData,googleapi_key::String;
                            googleapi_parameters::Dict{Symbol,String} = googleAPI_parameters)
    url = OpenStreetMapX.get_googleapi_url(origin, destination,map_data,googleapi_key,googleapi_parameters = googleapi_parameters)
    status, routes = OpenStreetMapX.parse_google_url(url)
    if status == "OK"
        route = OpenStreetMapX.extract_google_route(routes[1])
        return OpenStreetMapX.google_route_to_network(route,map_data), "google"
    elseif status =="ZERO_RESULTS"
        return Int[],"google"
    elseif status == "OVER_QUERY_LIMIT"
        sleep(0.5)
        return OpenStreetMapX.get_google_route(origin,destination,map_data,googleapi_key,googleapi_parameters = googleapi_parameters)
    else
        #get route based on OSM routing
        @warn "Google Distances API cannot get a proper results - route will be calculated with OSMSim Routing module"
		if rand() < 0.5
			route_nodes, distance, route_time = OpenStreetMapX.shortest_route(map_data, origin, destination)
			return route_nodes, "shortest"
		else
			route_nodes, distance, route_time = OpenStreetMapX.fastest_route(map_data, origin, destination)
			return route_nodes, "fastest"
		end
    end
end


text:
"""
    Parities

Implementation of parity queries using Fast Hadamard Walsh Transform and Gosper
iteration.

Releasing privacy prserving marginals for contingency table.
The measurments in this case are not the marginals, but the fourier coefficients of the contingency table,
that are sufficient for recovering the marginals.

Details can be found in:
"B. Barak, K. Chaudhuri, C. Dwork, S. Kale,F. McSherry, and K. Talwar.
Privacy, accuracy, and consistency too: a holistic solution to contingency table release. In PODS, 2007."
"""
struct Parities <: Queries
    dimension::Int
    order::Int
    idx::Array{Int, 1}
end

"""
    hadamard_basis_vector(index, dimension)

Returns hadamard basis vector given by `index`.
"""
function hadamard_basis_vector(index::Int, dimension::Int)
    hadamard = zeros(Float64, 1 << dimension)
    hadamard[1] = 1.0
    for i = 0:dimension-1
        sign = (index & (1 << i)) > 0 ? -1.0 : 1.0
        @simd for j = 1:(1 << i)
            @inbounds hadamard[j + (1 << i)] = sign * hadamard[j]
        end
    end
    hadamard
end

"""
return the downward closure of all the incides of the k-way marginals, of dimension d.

# Example:
    k = 2, d = 4
    k-way indices= [0011, 0101, 0110, 1001, 1010, 1100]
    downward_closure(0011) = [0000, 0001, 0010, 0011]
"""
function complete_way_marginals_indices(k::Int64, d::Int64)

    idx = Set(zeros(Int64, 0))

    for s = gosper(d, k)
        for alpha in BinaryItr(d, s)
            temp = alpha + 1
            push!(idx, temp)
        end
    end

    sort(collect(idx))
end

function Parities(dimension, order)
    Parities(dimension, order, complete_way_marginals_indices(order, dimension))
end

function get(queries::Parities, i::Int)
    HistogramQuery(hadamard_basis_vector(queries.idx[i]-1, queries.dimension))
end

"""
evaluate the measurements of the linear queries defined by
queries.idx (the indices needed for calculating the marginals of the contingency table)
"""
function evaluate(queries::Parities, h::Histogram)
    2^queries.dimension * fwht_natural(h.weights)[queries.idx]
end

function fourierCoefficients(queries::Parities, h::Histogram)
    2^queries.dimension * fwht_natural(h.weights)
end

"""
Calculate the marginal C_beta : R^(2^d) -> R^(2^(norm1(beta))), for beta in {0,1}^d,
of the fourier basis vector in the location alpha
"""
function fourierMarginal(dimension::Int64, beta::Int64, alpha::Int64)

    @assert beta & alpha == alpha

    beta_norm1, beta_indices = norm_1_with_indices(beta)
    alpha_norm1, alpha_indices = norm_1_with_indices(alpha)

    marginal = zeros(Int64, 2^beta_norm1)

    inverse_beta = setdiff(collect(1:dimension), beta_indices)

    for (gamma_index, gamma) in enumerate(BinaryItr(dimension, beta))
        marginal[gamma_index] = 2^length(inverse_beta)*((-1)^(norm_1(alpha&gamma)))
    end

    marginal
end

"""
Calculate the marginal C_beta : R^(2^d) -> R^(2^(norm1(beta))), for beta in {0,1}^d,
of the contingencyTable (represent as histogram).
"""
function calculateMarginal(queries::Parities, h::Histogram, beta::Int64)

    beta_norm1, beta_indices = norm_1_with_indices(beta)

    marginal = zeros(Int64, 2^beta_norm1)

    coefficients = fourierCoefficients(queries, h)

    for alpha in BinaryItr(queries.dimension, beta)
        marginal += fourierMarginal(queries.dimension, beta, alpha)*coefficients[alpha+1]
    end

    marginal/(2^queries.dimension)
end

function Marginals(queries::Parities, h::Histogram)

    marginals = []

    for beta in queries.idx
        push!(marginals, calculateMarginal(queries, h, beta))
    end

    marginals
end

function initialize(queries::Parities, data::Tabular, ps::MWParameters)
    initialize(queries, Histogram(data), ps)
end

struct FactorParity <: FactorHistogramQuery
    attributes::Array{Int, 1}
end

"""
    FactorParities

Implementation of parity queries over factored histograms.
"""
struct FactorParities <: FactorHistogramQueries
    queries::Array{FactorParity, 1}
end

attributes(q::FactorParity) = q.attributes
default_value(q::FactorParity) = 0.0
get(qs::FactorParities, i::QueryIndex) = qs.queries[i]

function restrict(q::FactorParity, attributes::Array{Int, 1})
    idx = 0
    d = length(attributes)
    for a in q.attributes
        i = findfirst(x -> x==a, attributes)
        idx += 2^(d-i)
    end
    HistogramQuery(hadamard_basis_vector(idx, d))
end

function FactorParities(dimension, order)
    qs = FactorParity[]
    for k = 1:order
        append!(qs, [FactorParity(s) for s = subsets(collect(1:dimension), k)])
    end
    FactorParities(qs)
end

function parity(x::Array{Float64, 1}, attributes::Array{Int64, 1})
    (-1.0)^sum(x[attributes])
end

function evaluate(q::FactorParity, table::Tabular)
    s = 0.0
    n = size(table.data)[2]
    a = attributes(q)
    @simd for i = 1:n
        @inbounds s += parity(table.data[:, i], a)
    end
    s/n
end

function evaluate(qs::FactorParities, table::Tabular)
    evals = zeros(Float64, length(qs.queries))
    @simd for i in 1:length(qs.queries)
        @inbounds evals[i] = evaluate(qs.queries[i], table)
    end
    evals
end

function evaluate(qs::FactorParities, d::FactorHistogram)
    [evaluate(q, d) for q in qs.queries]
end


text:
#=
Visual illustration of when the law of large numbers fails

@author : Spencer Lyon <spencer.lyon@nyu.edu>

References
----------
Based off the original python file cauchy_samples.py
=#
using PyPlot
using Distributions

srand(12)  # reproducible results
n = 200
dist = Cauchy()
data = rand(dist, n)

function plot_draws()
    fig, ax = subplots()
    ax[:plot](1:n, data, "bo", alpha=0.5)
    ax[:vlines](1:n, 0, data, lw=0.2)
    ax[:set_title]("$n observations from the Cauchy distribution")
    Void
end

function plot_means()
    # == Compute sample mean at each n == #
    sample_mean = Array(Float64, n)
    for i=1:n
        sample_mean[i] = mean(data[1:i])
    end

    # == Plot == #
    fig, ax = subplots()
    ax[:plot](1:n, sample_mean, "r-", lw=3, alpha=0.6, label=L"$\bar{X}_n$")
    ax[:plot](1:n, zeros(n), "k--", lw=0.5)
    ax[:legend]()
    Void
end

plot_draws()
plot_means()