text:
module Wine
    using Gallium
    using Gallium: load, RemotePtr

    const REQ_new_process                  =  0
    const REQ_get_new_process_info         =  1
    const REQ_new_thread                   =  2
    const REQ_get_startup_info             =  3
    const REQ_init_process_done            =  4
    const REQ_init_thread                  =  5
    const REQ_terminate_process            =  6
    const REQ_terminate_thread             =  7
    const REQ_get_process_info             =  8
    const REQ_set_process_info             =  9
    const REQ_get_thread_info              = 10
    const REQ_get_thread_times             = 11
    const REQ_set_thread_info              = 12
    const REQ_get_dll_info                 = 13
    const REQ_suspend_thread               = 14
    const REQ_resume_thread                = 15
    const REQ_load_dll                     = 16
    const REQ_unload_dll                   = 17
    const REQ_queue_apc                    = 18
    const REQ_get_apc_result               = 19
    const REQ_close_handle                 = 20
    const REQ_set_handle_info              = 21
    const REQ_dup_handle                   = 22
    const REQ_open_process                 = 23
    const REQ_open_thread                  = 24
    const REQ_select                       = 25
    const REQ_create_event                 = 26
    const REQ_event_op                     = 27
    const REQ_query_event                  = 28
    const REQ_open_event                   = 29
    const REQ_create_keyed_event           = 30
    const REQ_open_keyed_event             = 31
    const REQ_create_mutex                 = 32
    const REQ_release_mutex                = 33
    const REQ_open_mutex                   = 34
    const REQ_create_semaphore             = 35
    const REQ_release_semaphore            = 36
    const REQ_query_semaphore              = 37
    const REQ_open_semaphore               = 38
    const REQ_create_file                  = 39
    const REQ_open_file_object             = 40
    const REQ_alloc_file_handle            = 41
    const REQ_get_handle_unix_name         = 42
    const REQ_get_handle_fd                = 43
    const REQ_flush                        = 44
    const REQ_lock_file                    = 45
    const REQ_unlock_file                  = 46
    const REQ_create_socket                = 47
    const REQ_accept_socket                = 48
    const REQ_accept_into_socket           = 49
    const REQ_set_socket_event             = 50
    const REQ_get_socket_event             = 51
    const REQ_get_socket_info              = 52
    const REQ_enable_socket_event          = 53
    const REQ_set_socket_deferred          = 54
    const REQ_alloc_console                = 55
    const REQ_free_console                 = 56
    const REQ_get_console_renderer_events  = 57
    const REQ_open_console                 = 58
    const REQ_get_console_wait_event       = 59
    const REQ_get_console_mode             = 60
    const REQ_set_console_mode             = 61
    const REQ_set_console_input_info       = 62
    const REQ_get_console_input_info       = 63
    const REQ_append_console_input_history = 64
    const REQ_get_console_input_history    = 65
    const REQ_create_console_output        = 66
    const REQ_set_console_output_info      = 67
    const REQ_get_console_output_info      = 68
    const REQ_write_console_input          = 69
    const REQ_read_console_input           = 70
    const REQ_write_console_output         = 71
    const REQ_fill_console_output          = 72
    const REQ_read_console_output          = 73
    const REQ_move_console_output          = 74
    const REQ_send_console_signal          = 75
    const REQ_read_directory_changes       = 76
    const REQ_read_change                  = 77
    const REQ_create_mapping               = 78
    const REQ_open_mapping                 = 79
    const REQ_get_mapping_info             = 80
    const REQ_get_mapping_committed_range  = 81
    const REQ_add_mapping_committed_range  = 82
    const REQ_create_snapshot              = 83
    const REQ_next_process                 = 84
    const REQ_next_thread                  = 85
    const REQ_wait_debug_event             = 86
    const REQ_queue_exception_event        = 87
    const REQ_get_exception_status         = 88
    const REQ_continue_debug_event         = 89
    const REQ_debug_process                = 90
    const REQ_debug_break                  = 91
    const REQ_set_debugger_kill_on_exit    = 92
    const REQ_read_process_memory          = 93
    const REQ_write_process_memory         = 94
    const REQ_create_key                   = 95
    const REQ_open_key                     = 96
    const REQ_delete_key                   = 97
    const REQ_flush_key                    = 98
    const REQ_enum_key                     = 99
    const REQ_set_key_value                = 100
    const REQ_get_key_value                = 101
    const REQ_enum_key_value               = 102
    const REQ_delete_key_value             = 103
    const REQ_load_registry                = 104
    const REQ_unload_registry              = 105
    const REQ_save_registry                = 106
    const REQ_set_registry_notification    = 107
    const REQ_create_timer                 = 108
    const REQ_open_timer                   = 109
    const REQ_set_timer                    = 110
    const REQ_cancel_timer                 = 111
    const REQ_get_timer_info               = 112
    const REQ_get_thread_context           = 113
    const REQ_set_thread_context           = 114
    const REQ_get_selector_entry           = 115
    const REQ_add_atom                     = 116
    const REQ_delete_atom                  = 117
    const REQ_find_atom                    = 118
    const REQ_get_atom_information         = 119
    const REQ_set_atom_information         = 120
    const REQ_empty_atom_table             = 121
    const REQ_init_atom_table              = 122
    const REQ_get_msg_queue                = 123
    const REQ_set_queue_fd                 = 124
    const REQ_set_queue_mask               = 125
    const REQ_get_queue_status             = 126
    const REQ_get_process_idle_event       = 128
    const REQ_send_message                 = 129
    const REQ_post_quit_message            = 130
    const REQ_send_hardware_message        = 131
    const REQ_get_message                  = 132
    const REQ_reply_message                = 133
    const REQ_accept_hardware_message      = 134
    const REQ_get_message_reply            = 135
    const REQ_set_win_timer                = 136
    const REQ_kill_win_timer               = 137
    const REQ_is_window_hung               = 138
    const REQ_get_serial_info              = 139
    const REQ_set_serial_info              = 140
    const REQ_register_async               = 141
    const REQ_cancel_async                 = 142
    const REQ_read                         = 143
    const REQ_write                        = 144
    const REQ_ioctl                        = 145
    const REQ_set_irp_result               = 146
    const REQ_get_irp_result               = 147
    const REQ_create_named_pipe            = 148
    const REQ_get_named_pipe_info          = 149
    const REQ_set_named_pipe_info          = 150
    const REQ_create_window                = 151
    const REQ_destroy_window               = 152
    const REQ_get_desktop_window           = 153
    const REQ_set_window_owner             = 154
    const REQ_get_window_info              = 155
    const REQ_set_window_info              = 156
    const REQ_set_parent                   = 157
    const REQ_get_window_parents           = 158
    const REQ_get_window_children          = 159
    const REQ_get_window_children_from_point = 160
    const REQ_get_window_tree                = 161
    const REQ_set_window_pos                 = 162
    const REQ_get_window_rectangles          = 163
    const REQ_get_window_text                = 164
    const REQ_set_window_text                = 165
    const REQ_get_windows_offset             = 166
    const REQ_get_visible_region             = 167
    const REQ_get_surface_region             = 168
    const REQ_get_window_region              = 169
    const REQ_set_window_region              = 170
    const REQ_get_update_region              = 171
    const REQ_update_window_zorder           = 172
    const REQ_redraw_window                  = 173
    const REQ_set_window_property            = 174
    const REQ_remove_window_property         = 175
    const REQ_get_window_property            = 176
    const REQ_get_window_properties          = 177
    const REQ_create_winstation              = 178
    const REQ_open_winstation                = 179
    const REQ_close_winstation               = 180
    const REQ_get_process_winstation         = 181
    const REQ_set_process_winstation         = 182
    const REQ_enum_winstation                = 183
    const REQ_create_desktop                 = 184
    const REQ_open_desktop                   = 185
    const REQ_open_input_desktop             = 186
    const REQ_close_desktop                  = 187
    const REQ_get_thread_desktop             = 188
    const REQ_set_thread_desktop             = 189
    const REQ_enum_desktop                   = 190
    const REQ_set_user_object_info           = 191
    const REQ_register_hotkey                = 192
    const REQ_unregister_hotkey              = 193
    const REQ_attach_thread_input            = 194
    const REQ_get_thread_input               = 195
    const REQ_get_last_input_time            = 196
    const REQ_get_key_state                  = 197
    const REQ_set_key_state                  = 198
    const REQ_set_foreground_window          = 199
    const REQ_set_focus_window               = 200
    const REQ_set_active_window              = 201
    const REQ_set_capture_window             = 202
    const REQ_set_caret_window               = 203
    const REQ_set_caret_info                 = 204
    const REQ_set_hook                       = 205
    const REQ_remove_hook                    = 206
    const REQ_start_hook_chain               = 207
    const REQ_finish_hook_chain              = 208
    const REQ_get_hook_info                  = 209
    const REQ_create_class                   = 210
    const REQ_destroy_class                  = 211
    const REQ_set_class_info                 = 212
    const REQ_set_clipboard_info             = 213
    const REQ_empty_clipboard                = 214
    const REQ_open_token                     = 215
    const REQ_set_global_windows             = 216
    const REQ_adjust_token_privileges        = 217
    const REQ_get_token_privileges           = 218
    const REQ_check_token_privileges         = 219
    const REQ_duplicate_token                = 220
    const REQ_access_check                   = 221
    const REQ_get_token_sid                  = 222
    const REQ_get_token_groups               = 223
    const REQ_get_token_default_dacl         = 224
    const REQ_set_token_default_dac          = 225
    const REQ_set_security_object            = 226
    const REQ_get_security_object            = 227
    const REQ_get_system_handles             = 228
    const REQ_create_mailslot                = 229
    const REQ_set_mailslot_info              = 230
    const REQ_create_directory               = 231
    const REQ_open_directory                 = 232
    const REQ_get_directory_entry            = 233
    const REQ_create_symlink                 = 234
    const REQ_open_symlink                   = 235
    const REQ_query_symlink                  = 236
    const REQ_get_object_info                = 237
    const REQ_get_object_type                = 238
    const REQ_unlink_object                  = 239
    const REQ_get_token_impersonation_level  = 240
    const REQ_allocate_locally_unique_id     = 241
    const REQ_create_device_manager          = 242
    const REQ_create_device                  = 243
    const REQ_delete_device                  = 244
    const REQ_get_next_device_request        = 245
    const REQ_make_process_system            = 246
    const REQ_get_token_statistics           = 247
    const REQ_create_completion              = 248
    const REQ_open_completion                = 249
    const REQ_add_completion                 = 250
    const REQ_remove_completion              = 251
    const REQ_query_completion               = 252
    const REQ_set_completion_info            = 253
    const REQ_add_fd_completion              = 254
    const REQ_set_fd_disp_info               = 255
    const REQ_set_fd_name_info               = 256
    const REQ_get_window_layered_info        = 257
    const REQ_set_window_layered_info        = 258
    const REQ_alloc_user_handle              = 259
    const REQ_free_user_handle               = 260
    const REQ_set_cursor                     = 261
    const REQ_update_rawinput_devices        = 262
    const REQ_get_suspend_context            = 263
    const REQ_set_suspend_context            = 264
    const REQ_create_job                     = 265
    const REQ_open_job                       = 266
    const REQ_assign_job                     = 267
    const REQ_process_in_job                 = 268
    const REQ_set_job_limits                 = 269
    const REQ_set_job_completion_port        = 270
    const REQ_terminate_job                  = 271

    immutable request_header
        req::Cint
        request_size::Cint
        reply_size::Cint
    end

    immutable init_thread_request
        header::request_header
        unix_pid::Cint
        unix_tid::Cint
        debug_level::Cint
        teb::RemotePtr{Void}
        entry::RemotePtr{Void}
        reply_fd::Cint
        wait_fd::Cint
        cpu::Cint
    end

    function load_request(vm, addr)
        req = load(vm, RemotePtr{Cint}(addr))
        @show req
        if req == REQ_init_thread
            return load(vm, RemotePtr{init_thread_request}(addr))
        else
            error("unknown request")
        end
    end

end # module


text:
#
#	Process
#
function create_grids(
	N::Int,
	l::AbstractFloat)

	# compute derived parameters:
	Δt = 2*l/N
	h  = π*N/2/l
	Δv = 2*h/N

	# compute grids:
	t  = -l .+ collect(0:N-1).*Δt
	v  = -h .+ collect(0:N-1).*Δv

	return t, v
end

function phase_shift!(
	x::Vector{Complex{T}} where T<:AbstractFloat)

	for i in 1:length(x)
		x[i] *= (-1)^(i-1)
	end

	x
end

function ift!(
	G::Array{Complex{T}} where T<:AbstractFloat,
	F::Array{Complex{T}} where T<:AbstractFloat,
	N::Integer,
	l::AbstractFloat)

	G .= F
	phase_shift!(ifft!(phase_shift!(G)))
	G .*= sqrt(pi/2)*(N/l)
end

#
#	Field
#

# Phase shift
function phase_shift!(Y::Cscalar)

	for i=1:Npoints(Y)
		Y.values[i] *= (-1)^sum(tovec(i,Y.cl.N))
	end
end

# ift: F ↦ G
function ift!(G::Cscalar, F::Cscalar)

	G.values .= F.values

	dims = Tuple(n for n in G.cl.N)

	phase_shift!(G)
	ifft!(reshape(G.values,dims))
	phase_shift!(G)

	L = G.cl.points[end,1:end]
	N = G.cl.N
	d = dim(G)
	G.values .*= (pi/2)^(d/2)*prod(N)/prod(L)

end

# inverse Hankel?

ext: jl


text:
# global variable tracking the number of bytes written in the current execution
# request
const stdio_bytes = Ref(0)


text:
const UrbanState = Records.Frame{Records.Entity{AutomotiveDrivingModels.VehicleState,AutomotiveDrivingModels.VehicleDef,Int64}}

const UrbanObs = Vector{Float64}

#### Action type

struct UrbanAction
    acc::Float64
end

function Base.copyto!(a::UrbanAction, b::UrbanAction)
    a.acc = b.acc
end

function Base.hash(a::UrbanAction, h::UInt64 = zero(UInt64))
    return hash(a.acc, h)
end

function Base.:(==)(a::UrbanAction, b::UrbanAction)
    return a.acc == b.acc
end

#### POMDP type

@with_kw mutable struct UrbanPOMDP <: POMDP{UrbanState, UrbanAction, UrbanObs}
    env::UrbanEnv = UrbanEnv()
    obs_dist::ObstacleDistribution = ObstacleDistribution(UrbanEnv())
    lidar::Bool = false
    sensor::AbstractSensor = GaussianSensor()
    models::Dict{Int64, DriverModel} = Dict{Int64, DriverModel}(1=>EgoDriver(UrbanAction(0.)))
    ego_type::VehicleDef = VehicleDef()
    car_type::VehicleDef = VehicleDef()
    ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)
    max_cars::Int64 = 10
    max_peds::Int64 = 10
    max_obstacles::Int64 = 3
    max_acc::Float64 = 2.0
    ego_start::Float64 = env.params.stop_line - ego_type.length/2
    ego_goal::LaneTag = LaneTag(2,2)
    off_grid::VecSE2 = VecSE2(UrbanEnv().params.x_min+VehicleDef().length/2, UrbanEnv().params.y_min+VehicleDef().width/2, 0)
    car_models::Dict{SVector{2, LaneTag}, DriverModel} = get_car_models(env, get_ttc_model)
    n_features::Int64 = 4
    ΔT::Float64 = 0.5 # decision frequency
    car_birth::Float64 = 0.3
    ped_birth::Float64 = 0.3
    a_noise::Float64 = 1.0
    v_noise::Float64 = 1.
    pos_obs_noise::Float64 = 0.5
    vel_obs_noise::Float64 = 0.5
    collision_cost::Float64 = -1.
    action_cost::Float64 = 0.
    goal_reward::Float64 = 1.
    γ::Float64 = 0.95 # discount factor
end

## HELPERS

function POMDPs.discount(pomdp::UrbanPOMDP)
    return pomdp.γ
end

POMDPs.actions(pomdp::UrbanPOMDP) = [UrbanAction(-2.0), UrbanAction(0.0), UrbanAction(2.0), UrbanAction(4.0)]

function POMDPs.actionindex(pomdp::UrbanPOMDP, action::UrbanAction)
    if action.acc == -4.0
        return 1
    elseif action.acc == -2.0
        return 2
    elseif action.acc == 0.
        return 3
    else
        return 4
    end
end

# only works for single lane intersection!
function get_car_models(env::UrbanEnv, get_model::Function)
    d = Dict{SVector{2, LaneTag}, DriverModel}()

    r1 = SVector(LaneTag(1,1), LaneTag(2,1))
    d[r1] = get_model(env, r1)

    r2 = SVector(LaneTag(1,1), LaneTag(5,1))
    d[r2] = get_model(env, r2)

    r3 = SVector(LaneTag(3,1), LaneTag(4,1))
    d[r3] = get_model(env, r3)

    r4 = SVector(LaneTag(3,1), LaneTag(5,1))
    d[r4] = get_model(env, r4)
    return d
end


text:
module JuliaKara_noGUI

include("ActorsWorld.jl"); using .ActorsWorld

export
    World,
    place_kara,
    place_tree,
    place_leaf,
    place_mushroom,
    move,
    turnLeft,
    turnRight,
    putLeaf,
    removeLeaf,
    treeLeft,
    treeFront,
    treeRight,
    mushroomFront,
    onLeaf,
    load_world,
    save_world,
    get_kara,
    reset!,
    World_State,
    world_state_save

const ACTOR_DEFINITIONS = Dict(
    :kara => Actor_Definition(
        moveable=true,
        turnable=true
    ),
    :tree => Actor_Definition(),
    :mushroom => Actor_Definition(
        moveable=true
    ),
    :leaf => Actor_Definition(
        grabable=true,
        layer=0
    )
)

"""
    place_kara(wo::World,x::Int,y::Int,direction::Symbol=:NORTH)

Places kara in world `wo` at location `x`, `y` in direction `direction`.
`direction` is either of :NORTH, :EAST, :SOUTH, :WEST.
Returns a referenc to the created object.
"""
function place_kara(wo::World,x::Int,y::Int,direction::Symbol=ActorsWorld.DIRECTIONS[1])
    actor_create!(
        wo,
        ACTOR_DEFINITIONS[:kara],
        Location(x,y),
        Orientation(direction)
    )
end

"""
    place_tree(wo::World,x::Int,y::Int)

Places a tree in world `wo` at location `x`, `y`.
Returns a referenc to the created object.
"""
function place_tree(wo::World,x::Int,y::Int)
    actor_create!(
        wo,
        ACTOR_DEFINITIONS[:tree],
        Location(x,y),
        Orientation(ActorsWorld.DIRECTIONS[1])
    )
end

"""
    place_leaf(wo::World,x::Int,y::Int)

Places a leaf in world `wo` at location `x`, `y`.
Returns a referenc to the created object.
"""
function place_leaf(wo::World,x::Int,y::Int)
    actor_create!(
        wo,
        ACTOR_DEFINITIONS[:leaf],
        Location(x,y),
        Orientation(ActorsWorld.DIRECTIONS[1])
    )
end

"""
    place_mushroom(wo::World,x::Int,y::Int)

Places a mushroom in world `wo` at location `x`, `y`.
Returns a referenc to the created object.
"""
function place_mushroom(wo::World,x::Int,y::Int)
    actor_create!(
        wo,
        ACTOR_DEFINITIONS[:mushroom],
        Location(x,y),
        Orientation(ActorsWorld.DIRECTIONS[1])
    )
end

"""
    get_kara(wo::World)

Returns all actors of type kara in the world `wo`.
In case there is only one kara no vector is returned.
"""
function get_kara(wo::World)
    k_found = Actor[]
    for ac in wo.actors
        if ac.actor_definition == ACTOR_DEFINITIONS[:kara]
            push!(k_found,ac)
        end
    end
    # Return kara explicitly if only one was found
    # This is not a good style but in terms of usability
    # its better.
    return length(k_found) == 1 ? k_found[1] : k_found
end

"""
    move(wo::World,ac::Actor)

Moves the actor `ac` a step forward in the world `wo`.
"""
function move(wo::World,ac::Actor)
    actor_move!(wo,ac,ac.orientation.value)
end

"""
    turnLeft(wo::World,ac::Actor)

Turns the actor `ac` counter clockwise.
"""
function turnLeft(wo::World,ac::Actor)
    actor_rotate!(ac,false)
end

"""
    turnRight(wo::World,ac::Actor)

Turns the actor `ac` clockwise.
"""
function turnRight(wo::World,ac::Actor)
    actor_rotate!(ac,true)
end

"""
    removeLeaf(wo::World,ac::Actor)

Removes an actor of type leaf from the location `ac` is at.
"""
function removeLeaf(wo::World,ac::Actor)
    actor_pickup!(wo,ac)
end

"""
    putLeaf(wo::World,ac::Actor)

Places an actor of type leaf the location `ac` is at.
"""
function putLeaf(wo::World,ac::Actor)
    actor_putdown!(wo,ac,ACTOR_DEFINITIONS[:leaf])
end

"""
    treeLeft(wo::World,ac::Actor)

Checks if there is an actor of type tree left of actor `ac`.
"""
function treeLeft(wo::World,ac::Actor)
    is_actor_definition_left(wo,ac,ACTOR_DEFINITIONS[:tree])
end

"""
    treeRight(wo::World,ac::Actor)

Checks if there is an actor of type tree right of actor `ac`.
"""
function treeRight(wo::World,ac::Actor)
    is_actor_definition_right(wo,ac,ACTOR_DEFINITIONS[:tree])
end

"""
    treeFront(wo::World,ac::Actor)

Checks if there is an actor of type tree in front of actor `ac`.
"""
function treeFront(wo::World,ac::Actor)
    is_actor_definition_front(wo,ac,ACTOR_DEFINITIONS[:tree])
end

"""
    mushroomFront(wo::World,ac::Actor)

Checks if there is an actor of type mushroom in front of actor `ac`.
"""
function mushroomFront(wo::World,ac::Actor)
    is_actor_definition_front(wo,ac,ACTOR_DEFINITIONS[:mushroom])
end

"""
    onLeaf(wo::World,ac::Actor)

Checks if there is an actor of type leaf below of actor `ac`.
"""
function onLeaf(wo::World,ac::Actor)
    is_actor_definition_here(wo,ac,ACTOR_DEFINITIONS[:leaf])
end

include("JuliaKara_interface_xml.jl")
save_world = xml_save_world
load_world = xml_load_world

end


text:
importall Bukdu
importall Bukdu.Octo
using Base.Test

vector = Assoc(Vector([(:name, String)]))

@test String == vector[:name]
@test_throws KeyError vector[:age]

push!(vector, (:age, Int))

@test [:name, :age] == keys(vector)
@test [String, Int] == values(vector)
@test haskey(vector, :age)
@test !haskey(vector, :locker)

for (k,v) in vector
end

vector = Assoc(Vector([(:join, 0), (:where, 1), (:where, 2)]))
@test [:join, :where, :where] == keys(vector)
@test [0, 1, 2] == values(vector)
@test [:where, :where] == [k for (k,v) in vector if :where==k]
@test [1, 2] == [v for (k,v) in vector if :where==k]
@test [] == [v for (k,v) in vector if :locker==k]
@test haskey(vector, :where)

l = Assoc([(:author,"bar")])
r = Assoc([(:title,"title")])
@test r == Assoc([(:title,"title")])
@test Assoc([(:title,"title")]) == setdiff(r, l)

merge!(l, r)
@test Assoc([(:author,"bar"), (:title,"title")]) == l

@test !isempty(l)
empty!(l)
@test isempty(l)


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

function getversion(headerfile)
    a = b = c = 0
    for l in readlines(headerfile)
        if startswith(l, "#define JLCXX_VERSION_MAJOR")
            a = parse(Int, split(l)[end])
        end
        if startswith(l, "#define JLCXX_VERSION_MINOR")
            b = parse(Int, split(l)[end])
        end
        if startswith(l, "#define JLCXX_VERSION_PATCH")
            c = parse(Int, split(l)[end])
            break
        end
    end
    return VersionNumber(a, b, c)
end

basepath = dirname(@__DIR__)

name = "libcxxwrap_julia"
version = getversion(joinpath(basepath, "include", "jlcxx", "jlcxx_config.hpp"))

julia_versions = [v"1.6.0", v"1.7.0", v"1.8.0"]

# Collection of sources required to complete build
sources = [
    DirectorySource(basepath; target = "libcxxwrap-julia")
]

# Bash recipe for building across all platforms
script = raw"""
mkdir build
cd build

cmake \
    -DJulia_PREFIX=$prefix \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_FIND_ROOT_PATH=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    ../libcxxwrap-julia
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
install_license $WORKSPACE/srcdir/libcxxwrap-julia*/LICENSE.md
"""

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

function libjulia_platforms(julia_version)
    platforms = supported_platforms(; experimental = julia_version ≥ v"1.7")

    filter!(p -> libc(p) != "musl" && arch(p) != "i686" && !contains(arch(p), "arm") && arch(p) != "aarch64" && arch(p) != "powerpc64le" && !Sys.isfreebsd(p), platforms)

    for p in platforms
        p["julia_version"] = string(julia_version)
    end

    platforms = expand_cxxstring_abis(platforms)

    filter!(p -> cxxstring_abi(p) != "cxx03", platforms)

    return platforms
end

platforms = vcat(libjulia_platforms.(julia_versions)...)

# The products that we will ensure are always built
products = [
    LibraryProduct("libcxxwrap_julia", :libcxxwrap_julia; dlopen_flags = [:RTLD_GLOBAL]),
    LibraryProduct("libcxxwrap_julia_stl", :libcxxwrap_julia_stl; dlopen_flags = [:RTLD_GLOBAL]),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("libjulia_jll"),
]

GITHUB_REF_NAME = haskey(ENV, "GITHUB_REF_NAME") ? ENV["GITHUB_REF_NAME"] : ""
deployingargs = deepcopy(ARGS)
if !isempty(GITHUB_REF_NAME) && GITHUB_REF_NAME != "main"
    push!(deployingargs, "--deploy=local")
end

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

if GITHUB_REF_NAME == "main"
    const testreg = "https://github.com/barche/CxxWrapTestRegistry.git"

    const repo = "barche/libcxxwrap_julia_jll.jl"

    Pkg.Registry.add(RegistrySpec(url = testreg))
    Pkg.develop(PackageSpec(url = "https://github.com/$(repo).git"))

    libcxxwrap_jll_project_toml = TOML.parsefile(joinpath(Pkg.devdir(), "libcxxwrap_julia_jll", "Project.toml"))
    build_version = parse(VersionNumber, libcxxwrap_jll_project_toml["version"])
    if version == VersionNumber(build_version.major, build_version.minor, build_version.patch)
        build_version = VersionNumber(build_version.major, build_version.minor, build_version.patch, (), build_version.build .+ 1)
    else
        build_version = version
    end

    mktemp() do jsonfile, _
        push!(deployingargs, "--meta-json=$jsonfile")
        build_tarballs(deployingargs, name, version, sources, script, platforms, products, dependencies;
            preferred_gcc_version = v"9", julia_compat = "1.6")

        json = String(read(jsonfile))
        buff = IOBuffer(strip(json))
        objs = []
        while !eof(buff)
            push!(objs, BinaryBuilder.JSON.parse(buff))
        end
        json_obj = objs[1]

        BinaryBuilder.cleanup_merged_object!(json_obj)
        json_obj["dependencies"] = Dependency[dep for dep in json_obj["dependencies"] if !isa(dep, BuildDependency)]

        tag = "$(name)-v$(build_version)"
        upload_prefix = "https://github.com/$(repo)/releases/download/$(tag)"

        download_dir = "products"
        BinaryBuilder.rebuild_jll_package(json_obj; download_dir, upload_prefix, build_version)
        BinaryBuilder.push_jll_package(name, build_version; deploy_repo = repo)
        BinaryBuilder.upload_to_github_releases(repo, tag, download_dir; verbose = true)
    end

    import libcxxwrap_julia_jll
    using LocalRegistry
    register(libcxxwrap_julia_jll, registry = testreg)
end

max_stars_repo_name: mrdaly/ColorCode

max_issues_repo_name: mrdaly/ColorCode

max_forks_repo_name: mrdaly/ColorCode


text:
c=0
for n=1:100, k=1:n
  global c
  est = try
    binomial(n,k)
  catch
    2e6
  end
  if est > 1e6
    c+=1
  end
end
println(c)

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff


text:
using Pkg; Pkg.activate(".");
using StringDistances
using ConceptnetNumberbatch
using BenchmarkTools
using Random
using Languages
using Serialization

fid = open("./_conceptnet_/numberbatch-en-17.06.txt.bin")
cptnet = deserialize(fid)
close(fid)
words = (key for key in keys(cptnet[Languages.English()]) if isascii(key))

target="sstring"
for dist in [Jaro(), Levenshtein(), DamerauLevenshtein(), Cosine(), QGram(2), QGram(3)]
    _, idx = findmin(map(x->evaluate(dist, target, x), words))
    println("---------------")
    @time _, idx = findmin(map(x->evaluate(dist, target, x), words))
    println("[$dist], best match: $(collect(words)[idx])")
end


text:
# ---
# title: Nested layouts
# cover: assets/nested_layouts.png
# description: AlgebraOfGraphics plots within a Makie figure.
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# All AlgebraOfGraphics plots can be inserted in any figure position, where the rest
# of the figure is managed by vanilla Makie.
# For example

resolution = (800, 600)
fig = Figure(; resolution)
ax = Axis(fig[1, 1], title="Some plot")

df = (
    x=rand(500),
    y=rand(500),
    i=rand(["a", "b", "c"], 500),
    j=rand(["d", "e", "f"], 500),
    k=rand(Bool, 500),
    l=rand(Bool, 500)
)
plt = data(df) * mapping(:x, :y, col=:i, row=:j, color=:k, marker=:l)

subfig = fig[1, 2:3]
ag = draw!(subfig, plt)
for ae in ag
    ae.axis.xticklabelrotation[] = π/2
end
legend!(fig[end+1, 2], ag, orientation=:horizontal, tellheight=true)
fig

# The above also works in more nested situations.

resolution = (800, 600)
f = Figure(; resolution)
ax = Axis(f[1, 1], title="Some plot")
subfig = f[1, 2]
ax2 = Axis(subfig[1, 1])

df = (
    x=rand(500),
    y=rand(500),
    c=rand(["a", "b", "c"], 500),
)
plt = data(df) * mapping(:x, :y, color=:c)

draw!(subfig[2, 1], plt)
f

# It is also possible to let Makie control the axis and plot directly on top of it.

resolution = (800, 300)
f = Figure(; resolution)
ax1 = Axis(f[1, 1])
ax2 = Axis(f[1, 2])

df = (x=rand(100), y=rand(100), c=rand(["a", "b", "c"], 100))
plt = data(df) * mapping(:x, :y, color=:c)

scatter!(ax1, rand(10), rand(10), color=:black)
grid = draw!(ax2, plt)
legend!(f[1, 3], grid)

f

# save cover image #src
mkpath("assets") #src
save("assets/nested_layouts.png", fig) #src


text:
type PMLCoef
     vz::SparseMatrixCSC{Float64, Int64}
     vx::SparseMatrixCSC{Float64, Int64}
     pz::SparseMatrixCSC{Float64, Int64}
     px::SparseMatrixCSC{Float64, Int64}
end

function InitCpml(nz::Int64, nx::Int64, ext::Int64, dz::Float64, dx::Float64, vmax::Float64, iflag::Int64)
      if ext == 10
         R = 0.01
      elseif ext == 20
         R = 0.001
      elseif ext == 30
         R = 0.0001
      else
         error("unsupported damping layers")
      end
      ll = 2;
      if iflag == 1
         Nz = nz + 2*ext
         Nx = nx + 2*ext
         vz = spzeros(Nz, Nx)
         vx = spzeros(Nz, Nx)
         pz = spzeros(Nz, Nx)
         px = spzeros(Nz, Nx)
         #  upper and left
         for i = 1 : ext
             vz[i,:] = -1.5 * vmax/(ext*dz) * log(R) *  ((ext+1-i)/ext)^ll
             vx[:,i] = -1.5 * vmax/(ext*dx) * log(R) *  ((ext+1-i)/ext)^ll
             pz[i,:] = -1.5 * vmax/(ext*dz) * log(R) *  ((ext+1-i)/ext)^ll
             px[:,i] = -1.5 * vmax/(ext*dx) * log(R) *  ((ext+1-i)/ext)^ll
         end
         #  lower
         for iz = nz+ext+1 : Nz
             vz[iz,:] = -1.5 * vmax/(ext*dz) * log(R) * ((iz-nz-ext)/ext)^ll
             pz[iz,:] = -1.5 * vmax/(ext*dz) * log(R) * ((iz-nz-ext)/ext)^ll
         end
         #  right
         for ix = nx+ext+1 : Nx
             vx[:,ix] = -1.5 * vmax/(ext*dx) * log(R) * ((ix-nx-ext)/ext)^ll
             px[:,ix] = -1.5 * vmax/(ext*dx) * log(R) * ((ix-nx-ext)/ext)^ll
         end
      elseif iflag == 2
         Nz = nz +   ext
         Nx = nx + 2*ext
         vz = spzeros(Nz, Nx)
         vx = spzeros(Nz, Nx)
         pz = spzeros(Nz, Nx)
         px = spzeros(Nz, Nx)
         #  lower
         for iz = nz+1 : Nz
             vz[iz,:] = -1.5 * vmax/(ext*dz) * log(R) * ((iz-nz)/ext)^ll
             pz[iz,:] = -1.5 * vmax/(ext*dz) * log(R) * ((iz-nz)/ext)^ll
         end
         #  left
         for ix = 1 : ext
             vx[:,ix] = -1.5 * vmax/(ext*dx) * log(R) * ((ext+1-ix)/ext)^ll
             px[:,ix] = -1.5 * vmax/(ext*dx) * log(R) * ((ext+1-ix)/ext)^ll
         end
         #  right
         for ix = nx+ext+1 : Nx
             vx[:,ix] = -1.5 * vmax/(ext*dx) * log(R) * ((ix-nx-ext)/ext)^ll
             px[:,ix] = -1.5 * vmax/(ext*dx) * log(R) * ((ix-nx-ext)/ext)^ll
         end
      end
      Cpml = PMLCoef(vz, vx, pz, px)
      return Cpml
end

# =============================partial derivative matrix========================
type FidMtx
     nz :: Int64
     nx :: Int64
     ext:: Int64
     iflag :: Int64
     dz :: Float64
     dx :: Float64
     dt :: Float64
     MvzBvz :: SparseMatrixCSC{Float64,Int64}
     MvzBp  :: SparseMatrixCSC{Float64,Int64}
     MvxBvx :: SparseMatrixCSC{Float64,Int64}
     MvxBp  :: SparseMatrixCSC{Float64,Int64}
     MpzBpz :: SparseMatrixCSC{Float64,Int64}
     MpzBvz :: SparseMatrixCSC{Float64,Int64}
     MpxBpx :: SparseMatrixCSC{Float64,Int64}
     MpxBvx :: SparseMatrixCSC{Float64,Int64}
end

function modExpand(par::Array{Float64,2}, ext::Int64, iflag::Int64)
    temp1 = repmat(par[:,  1], 1, ext)
    temp2 = repmat(par[:,end], 1, ext)
    par   = hcat(temp1, par, temp2)
    if iflag == 1
       temp1 = repmat(par[1,  :]', ext, 1)
       temp2 = repmat(par[end,:]', ext, 1)
       par   = vcat(temp1, par, temp2)
    elseif iflag == 2
       temp2 = repmat(par[end,:]', ext, 1)
       par   = vcat(par, temp2)
    end
    return par
end

function Mvx(vx_pml::SparseMatrixCSC{Float64,Int64}, dx::Float64, dt::Float64, ext::Int64)
    (m,n) = size(vx_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dx; c2 = a2/dx;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt + vx_pml[iz,ix]/2
            b = 1/dt - vx_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = 1 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(n, n)
    tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
    for ix = 2: n-2
        tmp[ix,ix  ] = -c1; tmp[ix, ix+1] = c1;
        tmp[ix,ix-1] = -c2; tmp[ix, ix+2] = c2;
    end
    tmp[n-1,n-1] = -c1; tmp[n-1, n  ] = c1; tmp[n-1,n-2] = -c2;
    tmp[n  ,n  ] = -c1; tmp[n  , n-1] =-c2;
    MvxBvx = spdiagm(C3)
    MvxBp  = spdiagm(denum) * kron(tmp, speye(m))
    return MvxBvx, MvxBp
end

function Mvz(vz_pml::SparseMatrixCSC{Float64,Int64}, dz::Float64, dt::Float64, ext::Int64)
    (m,n) = size(vz_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dz; c2 = a2/dz;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt + vz_pml[iz,ix]/2
            b = 1/dt - vz_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = 1 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(m, m)
    tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
    for iz = 2: m-2
        tmp[iz,iz  ] = -c1; tmp[iz, iz+1] = c1;
        tmp[iz,iz-1] = -c2; tmp[iz, iz+2] = c2;
    end
    tmp[m-1,m-1] = -c1; tmp[m-1, m  ] = c1; tmp[m-1, m-2] = -c2;
    tmp[m  ,m  ] = -c1; tmp[m  , m-1] =-c2;
    MvzBvz = spdiagm(C3)
    MvzBp  = spdiagm(denum) * kron(speye(n), tmp)
    return MvzBvz, MvzBp
end

function Mpx(px_pml::SparseMatrixCSC{Float64,Int64}, v::Array{Float64,2}, dx::Float64, dt::Float64, ext::Int64, iflag::Int64)
    v = modExpand(v, ext, iflag)
    (m, n) = size(px_pml)
    a1 = 9/8  ; a2 = -1/24;
    c1 = a1/dx; c2 = a2/dx;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt + px_pml[iz,ix]/2
            b = 1/dt - px_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = (v[iz,ix])^2 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(n, n)
    tmp[1,1] =  c1; tmp[1,2] = c2;
    tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
    for ix = 3: n-1
        tmp[ix,ix-1] = -c1; tmp[ix, ix  ] = c1;
        tmp[ix,ix-2] = -c2; tmp[ix, ix+1] = c2;
    end
    tmp[n,n-1] = -c1; tmp[n,n] = c1; tmp[n,n-2] = -c2;
    MpxBpx = spdiagm(C3)
    MpxBvx = spdiagm(denum) * kron(tmp, speye(m))
    return MpxBpx, MpxBvx
end

function Mpz(pz_pml::SparseMatrixCSC{Float64,Int64}, v::Array{Float64,2}, dz::Float64, dt::Float64, ext::Int64, iflag::Int64)
    v = modExpand(v, ext, iflag)
    (m, n) = size(pz_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dz; c2 = a2/dz;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt + pz_pml[iz,ix]/2
            b = 1/dt - pz_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = (v[iz,ix])^2 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(m, m)
    tmp[1,1] =  c1; tmp[1,2] = c2;
    tmp[2,1] = -c1; tmp[2,2] = c1; tmp[2,3] = c2;
    for iz = 3: m-1
        tmp[iz,iz-1] = -c1; tmp[iz, iz  ] = c1;
        tmp[iz,iz-2] = -c2; tmp[iz, iz+1] = c2;
    end
    tmp[m,m-1] = -c1; tmp[m,m] = c1; tmp[m, m-2] = -c2
    MpzBpz = spdiagm(C3)
    MpzBvz = spdiagm(denum) * kron(speye(n), tmp)
    return MpzBpz, MpzBvz
end

function DispStable!(vmax::Float64, vmin::Float64, f0::Float64, dz::Float64, dx::Float64, dt::Float64)
    h = vmin / 5 / f0
    if dz >= h || dx >=h
       warn("large grid size may produce frequency dispersion, h should less than $h")
    end
    h = dx
    if dz < dx
       h = dz
    end
    tt = 6*h / (7*sqrt(3)*vmax)
    if dt >= tt
       warn("large time step may cause unstable, dt should less than $tt")
    end
    return nothing
end

function InitFidMtx(nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dz::Float64, dx::Float64, dt::Float64, vmax::Float64, vmin::Float64, f0::Float64, v::Array{Float64,2})
    DispStable!(vmax, vmin, f0, dz, dx, dt)
    Cpml = InitCpml(nz, nx, ext, dz, dx, vmax, iflag)
    (MvzBvz, MvzBp)  = Mvz(Cpml.vz,    dz, dt, ext       )
    (MvxBvx, MvxBp)  = Mvx(Cpml.vx,    dx, dt, ext       )
    (MpzBpz, MpzBvz) = Mpz(Cpml.pz, v, dz, dt, ext, iflag)
    (MpxBpx, MpxBvx) = Mpx(Cpml.px, v, dx, dt, ext, iflag)
    fidMtx = FidMtx(nz, nx, ext, iflag, dz, dx, dt, MvzBvz, MvzBp, MvxBvx, MvxBp, MpzBpz, MpzBvz, MpxBpx, MpxBvx)
    return fidMtx
end

# =========================forward modeling ====================================
type SnapShot
     isz:: Int64
     isx:: Int64
     nz :: Int64
     nx :: Int64
     ext:: Int64
     iflag :: Int64
     dt :: Float64
     it :: Int64
     vz :: Array{Float64, 1}
     vx :: Array{Float64, 1}
     pz :: Array{Float64, 1}
     px :: Array{Float64, 1}
end

function InitSnapShot(isz::Int64, isx::Int64, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dt::Float64, it::Int64)
    if iflag == 1
       Nz = nz + 2*ext
       Nx = nx + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       Nx = nx + 2*ext
    end
    spt = SnapShot(isz, isx, nz, nx, ext, iflag, dt, it, zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx))
    return spt
end

function Base.copy(snapShot1::SnapShot, snapShot2::SnapShot)
    if snapShot2.nz != snapShot1.nz || snapShot2.nx != snapShot1.nx || snapShot2.ext != snapShot1.ext || snapShot2.dt != snapShot1.dt || snapShot2.iflag != snapShot1.iflag
       error("the two snapShot are different")
    end
    snapShot1.isz = snapShot2.isz
    snapShot1.isx = snapShot2.isx
    snapShot1.it  = snapShot2.it
    snapShot1.vz[:] = snapShot2.vz[:]
    snapShot1.vx[:] = snapShot2.vx[:]
    snapShot1.pz[:] = snapShot2.pz[:]
    snapShot1.px[:] = snapShot2.px[:]
    return nothing
end

# =========================forward modeling ====================================
type Source
     iz :: Int64
     ix :: Int64
     dt :: Float64
     nt :: Int64
     p  :: Array{Float64,1}
end

function Ricker(f0::Float64, dt::Float64)
	  nw = 2.2/f0/dt
	  nw = 2*floor(Int,nw/2)+1
	  nc = floor(Int,nw/2)
	  w  = zeros(nw)
	  k  = collect(1:nw)
    k  = vec(k)
	  alpha = (nc-k+1)*f0*dt*pi
	  beta = alpha.^2
	  w = (1.-beta.*2).*exp(-beta)
	  return w
end

function InitSource(isz::Int64, isx::Int64, f0::Float64, dt::Float64)
    p = Ricker(f0, dt)
    nt= length(p)
    src = Source(isz, isx, dt, nt, p)
    return src
end

function AddSource!(spt::SnapShot, src::Source)
    dt = spt.dt
    if spt.dt != src.dt
       error("time sampling dismatch")
    end
    nz = spt.nz
    nx = spt.nx
    ext= spt.ext
    iflag = spt.iflag
    if iflag == 1
       zupper = ext
       Nz  = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz  = nz +   ext
    end
    Nx  = nx + 2*ext
    it  = spt.it
    if 0.0 <= (it-1)*dt <= (src.nt-1)*dt
       indt = it
       indz = src.iz + zupper
       indx = src.ix + ext
       pos = (indx-1) * Nz + indz
      #  spt.vz[pos] = spt.vz[pos] + src.p[indt]
       spt.pz[pos] = spt.pz[pos] + (src.p[indt])/2
       spt.px[pos] = spt.px[pos] + (src.p[indt])/2
    end
    return nothing
end

type ShotGather
     isz :: Int64
     isx :: Int64
     irz :: Array{Int64, 1}
     irx :: Array{Int64, 1}
     ot  :: Float64
     nt  :: Int64
     dt  :: Float64
     d   :: Array{Float64, 2}
end

function InitShotGather(isz::Int64, isx::Int64, pos::Array{Int64,2}, ot::Float64, nt::Int64, dt::Float64)
    nrec = size(pos, 1)
    irecz= pos[:, 1]
    irecx= pos[:, 2]
    d    = zeros(nt, nrec)
    shot = ShotGather(isz, isx, irecz, irecx, ot, nt, dt, d)
    return shot
end

function spt2shot!(shot::ShotGather, spt::SnapShot, Nz::Int64, Nx::Int64, ext::Int64, zupper::Int64)
    nrec = length(shot.irx)
    it = spt.it
    for irec = 1 : nrec
        iz = shot.irz[irec] + zupper
        ix = shot.irx[irec] + ext
        ind= (ix-1)*Nz + iz
        shot.d[it, irec] = spt.pz[ind] + spt.px[ind]
    end
    return nothing
end

function OneStepForward!(spt2::SnapShot, spt1::SnapShot, fidMtx::FidMtx)
    spt2.it  = spt1.it + 1
    spt2.vz = fidMtx.MvzBvz * spt1.vz + fidMtx.MvzBp  * (spt1.pz+spt1.px)
    spt2.vx = fidMtx.MvxBvx * spt1.vx + fidMtx.MvxBp  * (spt1.pz+spt1.px)
    spt2.pz = fidMtx.MpzBpz * spt1.pz + fidMtx.MpzBvz *  spt2.vz
    spt2.px = fidMtx.MpxBpx * spt1.px + fidMtx.MpxBvx *  spt2.vx
    return nothing
end

function MultiStepForward(pos::Array{Int64,2}, src::Source, fidMtx::FidMtx; tmax=2.0, interval=500, print_flag=false)
    nz = fidMtx.nz ; nx = fidMtx.nx;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    dt = fidMtx.dt
    nt = round(Int64, tmax/dt)+1
    isz= src.iz; isx= src.ix;
    shot = InitShotGather(isz, isx, pos, 0.0, nt, dt)
    tl = 0.0; tu = tl + (src.nt-1)*dt;
    spt1 = InitSnapShot(src.iz, src.ix, nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(src.iz, src.ix, nz, nx, ext, iflag, dt, 2)
    AddSource!(spt1, src)
    spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
    for it = 2: nt
        OneStepForward!(spt2, spt1, fidMtx)
        if tl <= (it-1)*dt <= tu   #Add sources
           AddSource!(spt2, src)
        end
        copy(spt1, spt2)
        spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    return shot
end

"""
    fidMtx = AcousticSetup(v, dz, dx, dt, f0, ext=20, iflag=1)

Setup sparse matrix for computing spatial derivation, return a composite type FidMtx

# Arguments
* `v   :: Array{Float64,2}`: P-wave velocity model
* `dz  :: Float64`: vertical grid size
* `dx  :: Float64`: horizontal grid size
* `dt  :: Float64`: size of time step
* `f0  :: Float64`: dominant frequency of Ricker wavelet

# keywords Arguments
* `ext=20` : number of PML absorbing boundry layers
* `iflag=1`: four sides(iflag=1) or three sides(iflag=2) absorbing boundries

# Output
* `fidMtx :: FidMtx`: composite type of sparse matrix
"""
function AcousticSetup(v::Array{Float64,2}, dz::Float64, dx::Float64, dt::Float64, f0::Float64; ext=30, iflag=1)
    (nz, nx) = size(v)
    vmax = maximum(v); vmin = minimum(v)
    fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v)
    return fidMtx
end

"""
    shot = SeisAcousticWave(fidMtx, pos, isz, isx, f0, dt, tmax=2.0)

finite difference modeling of acoustic wave field, generate a common shot gather

# Arguments
* `fidMtx :: FidMtx`        : composite type of sparse matrix
* `pos    :: Array{Int64,2}`: index of receivers, first column is the vertical index, second column contains horizontal index.
* `isz    :: Int64`         : vertical index of source
* `isz    :: Int64`         : horizontal index of source
* `f0     :: Float64`       : dominant frequency of Ricker wavelet
* `dt     :: Float64`       : size of time step

# keywords Arguments
* `tmax=1.0` : length of simulation

# Output
* `shot :: ShotGather`: composite type for common shot gather
"""
function SeisAcousticWave(fidMtx::FidMtx, pos::Array{Int64,2}, isz::Int64, isx::Int64, f0::Float64, dt::Float64; tmax=1.0, print_flag=false, interval=100)
    src  = InitSource(isz, isx, f0, dt)
    shot = MultiStepForward(pos, src, fidMtx, tmax=tmax, print_flag=print_flag, interval=interval)
    return shot
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
#=
The general outline of each tile is shown in beautiful ascii "art".

                  Tile
|----------------------------------------|
|                      | ^               |
|   Mappable           | |               |
|   |-----|            | |               |
|   |     |<-\         | |               |
|   |     |\  \        | |               |
|   |-----| \  \       | |               |
|            \  \      | |               |
|             \  \     | |               |
|              \  \    | |  Crossbar     |
|               \  \_|-----|             |
|                \-->|     |------------>|
|------------------->|     |<------------|
|<-------------------|-----|             |
|                      | ^               |
|                      | |               |
|                      | |               |
|                      | |               |
|----------------------------------------|

=#

# First, we have to build a simple Mappable component.
function build_mappable()
    # Start with a bare component
    # Need to indicate that this component is capable of holding a task, so we
    # initialize the metadata dictionary with some data we can use later.
    component = Component("simple"; metadata = Dict("mappable" => true))

    # Now all we need to do is add an input and output port.
    add_port(component, "in", Input)
    add_port(component, "out", Output)

    # And we're done, a simple component is build with just a single input and
    # output port.
    return component
end

# We now have to provide Mapper2 with some information regarding which components
# are mappable and which arent.
#
# We check the metadata dictionary for the "mappable" key. If it exists and is
# "true", than the compnent is mappable. Otherwise, it isn't. Simple :)
#
# Note that we MUST make the first argument of this function our locally defined
# subtype of Architecture to correctly extend the "ismappable" method.
Mapper2.ismappable(::Type{TutorialRuleset}, x) = get(x.metadata, "mappable", false)

# Since Crossbars are Mapper primitives, we are now ready to build our tile.
function build_tile()
    # Tiles are Components as well. Components are hierarchical in that
    # component may have subcomponent, each of which may have further more sub-
    # components etc.
    tile = Component("tile")

    # Declare the 4 input and 4 output ports
    for direction in ("north", "east", "south", "west")
        add_port(tile, "$(direction)_in", Input)
        add_port(tile, "$(direction)_out", Output)
    end

    # Instantiate the mappable primitive and add it as a subcomponent of the tile.
    # The instance name of the mappable component is "mappable"
    add_child(tile, build_mappable(), "mappable")

    # Next, we want to add the crossbar, which is unfortunately called "mux" in
    # the Mapper. Referencing the picture above, we want 5 inputs and 5 outputs.
    #
    # Note that the instance name of this component is "mux".
    add_child(tile, build_mux(5, 5), "mux")

    # Finally, start connecting ports together.
    # We first connect the "mappable" to "mux". Note that port names are
    # hierarchical based on instance names and port names.
    add_link(tile, "mappable.out", "mux.in[4]")
    add_link(tile, "mux.out[4]", "mappable.in")

    # Next, we connect the crossbar to the input and output ports of the tile.
    add_link(tile, "mux.out[0]", "north_out")
    add_link(tile, "mux.out[1]", "east_out")
    add_link(tile, "mux.out[2]", "south_out")
    add_link(tile, "mux.out[3]", "west_out")

    add_link(tile, "north_in", "mux.in[0]")
    add_link(tile, "east_in", "mux.in[1]")
    add_link(tile, "south_in", "mux.in[2]")
    add_link(tile, "west_in", "mux.in[3]")

    # And that's it!
    return tile
end

# Now that we've defined how to build a tile, we have to define how to build
# an architecture. This example will generate a 2D rectangular, rectilinear
# architecture with parameterized width and height.
function build_arch(width, height)
    # Initialize an empty TopLevel.
    # The parameter for the TopLevel is the dimensionality of the
    # architecture, which is 2.
    arch = TopLevel{2}("architecture")

    # Build a tile. We'll use this same tile at all locations of the grid.
    tile = build_tile()

    # Iterate over all rows and columns needed to achieve this height and width.
    # Build an address for this row and column and instantiate the tile.
    #
    # Since we don't provide the "add_child" method with a name for the
    # component, it will automatically generate one.
    for (row, col) in Iterators.product(1:width, 1:height)
        address = Address(row, col)
        add_child(arch, tile, address)
    end

    # Finally, we have to connect all these tiles together.
    # Simply define the Address offsets and source -> destination port pairs.
    offsets = [
        Offset((-1, 0), "north_out", "south_in"),
        Offset((1, 0), "south_out", "north_in"),
        Offset((0, 1), "east_out", "west_in"),
        Offset((0, -1), "west_out", "east_in"),
    ]

    connection_rule(arch, offsets)

    # And we're done!.
    return arch
end

# We don't want to have to worry about congestion in the final mapping, so
# lets set the capacity of all our routing resources to 100.
Mapper2.getcapacity(::Type{TutorialRuleset}, x) = 100


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

X = Dict(
    1 => [2.0, 3.0, 4.0],
    2 => [6.0, 3.0, 2.0],
    3 => [2.0, 5.0, 1.0],
    4 => [4.0, 3.0, 6.0])
u = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [0.0, 0.0, 0.0],
    3 => [0.0, 0.0, 0.0],
    4 => [0.25, 0.25, 0.25])
element1 = Element(Tet4, (1, 2, 3, 4))
element2 = Element(Tri3, (1, 2, 3))
update!((element1, element2), "geometry", X)
update!((element1, element2), "displacement", 0.0 => u)
update!(element1, "youngs modulus", 96.0)
update!(element1, "poissons ratio", 1.0/3.0)
update!(element1, "density", 420.0)
update!(element2, "displacement 1", 0.0)
update!(element2, "displacement 2", 0.0)
update!(element2, "displacement 3", 0.0)
problem1 = Problem(Elasticity, "test problem", 3)
problem1.properties.finite_strain = false
problem1.properties.geometric_stiffness = false
problem2 = Problem(Dirichlet, "boundary condition", 3, "displacement")
add_elements!(problem1, element1)
add_elements!(problem2, element2)
analysis = Analysis(Modal)
analysis.properties.which = :LM
add_problems!(analysis, problem1, problem2)

# test eigenvalues for single tet4 element
run!(analysis)
@test isapprox(analysis.properties.eigvals, [4/3, 1/3])

# test eigenvalues for single tet4 element, with geometric stiffness
problem1.properties.geometric_stiffness = true
analysis.properties.geometric_stiffness = true
run!(analysis)
@test isapprox(analysis.properties.eigvals, [5/3, 2/3])

# test poisson problem modal analysis without tie
X = Dict(
    1 => [0.0, 0.0],
    2 => [1.0, 0.0],
    3 => [1.0, 3.0],
    4 => [0.0, 3.0],
    5 => [0.0, 3.0],
    6 => [1.0, 3.0],
    7 => [1.0, 9.0],
    8 => [0.0, 9.0])
element1 = Element(Quad4, (1, 2, 3, 4))
element2 = Element(Quad4, (4, 3, 7, 8))
element3 = Element(Seg2, (1, 2))
element4 = Element(Seg2, (7, 8))
update!((element1, element2, element3, element4), "geometry", X)
update!((element1, element2), "density", 6.0)
update!((element1, element2), "thermal conductivity", 36.0)
update!((element3, element4), "temperature 1", 0.0)
problem1 = Problem(PlaneHeat, "combined body", 1)
problem2 = Problem(Dirichlet, "fixed ends", 1, "temperature")
add_elements!(problem1, element1, element2)
add_elements!(problem2, element3, element4)
analysis = Analysis(Modal)
add_problems!(analysis, problem1, problem2)
run!(analysis)
@test isapprox(first(analysis.properties.eigvals), 1.0)

#=
@testset "test poisson modal problem with mesh tie" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 3.0],
        4 => [0.0, 3.0],
        5 => [0.0, 3.0],
        6 => [1.0, 3.0],
        7 => [1.0, 9.0],
        8 => [0.0, 9.0])
    el1 = Element(Quad4, [1, 2, 3, 4])
    el2 = Element(Quad4, [5, 6, 7, 8])
    el3 = Element(Seg2, [1, 2])
    el4 = Element(Seg2, [7, 8])
    el5 = Element(Seg2, [3, 4])
    el6 = Element(Seg2, [5, 6])
    update!([el1, el2, el3, el4, el5, el6], "geometry", X)
    update!([el1, el2], "density", 6.0)
    update!([el1, el2], "thermal conductivity", 36.0)
    update!([el3, el4], "temperature 1", 0.0)
    update!(el5, "master elements", [el6])
    p1 = Problem(PlaneHeat, "body 1", 1)
    add_elements!(p1, [el1])
    p2 = Problem(PlaneHeat, "body 2", 1)
    add_elements!(p2, [el2])
    p3 = Problem(Dirichlet, "fixed ends", 1, "temperature")
    add_elements!(p3, [el3, el4])
    p4 = Problem(Mortar2D, "interface between bodies", 1, "temperature")
    add_slave_elements!(p4, [el5])
    add_master_elements!(p4, [el6])

    solver = Solver(Modal)
    push!(solver, p1, p2, p3, p4)
    solver()
    @test isapprox(solver.properties.eigvals[1], 1.0)
end
=#


text:
const pywrapcallback = pynew()
const pyjlcallbacktype = pynew()

pyjlcallback_repr(self) = Py("<jl $(repr(self))>")

pyjlcallback_str(self) = Py(string(self))

function pyjlcallback_call(self, args_::Py, kwargs_::Py)
    if pylen(kwargs_) > 0
        args = pyconvert(Vector{Py}, args_)
        kwargs = pyconvert(Dict{Symbol,Py}, kwargs_)
        Py(self(args...; kwargs...))
    elseif pylen(args_) > 0
        args = pyconvert(Vector{Py}, args_)
        Py(self(args...))
    else
        Py(self())
    end
end
pyjl_handle_error_type(::typeof(pyjlcallback_call), self, exc::MethodError) = exc.f === self ? pybuiltins.TypeError : PyNULL

function init_jlwrap_callback()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class CallbackValue(ValueBase):
        __slots__ = ()
        __module__ = "juliacall"
        def __repr__(self):
            if self._jl_isnull():
                return "<jl NULL>"
            else:
                return self._jl_callmethod($(pyjl_methodnum(pyjlcallback_repr)))
        def __str__(self):
            if self._jl_isnull():
                return "NULL"
            else:
                return self._jl_callmethod($(pyjl_methodnum(pyjlcallback_str)))
        def __call__(self, *args, **kwargs):
            return self._jl_callmethod($(pyjl_methodnum(pyjlcallback_call)), args, kwargs)
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlcallbacktype, jl.CallbackValue)
    pycopy!(pywrapcallback, pybuiltins.eval("lambda f: lambda *args, **kwargs: f(*args, **kwargs)", pydict()))
end

pyjlcallback(f) = pyjl(pyjlcallbacktype, f)
export pyjlcallback

function pycallback(f; name=nothing, doc=nothing)
    f2 = pyjlcallback(f)
    f3 = pywrapcallback(f2)
    pydel!(f2)
    if name !== nothing
        f3.__name__ = f3.__qualname__ = name
    end
    if doc !== nothing
        f3.__doc__ = doc
    end
    return f3
end
export pycallback


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.58980225336104 <= q <= 16.332536463331582)

                     add_NL_constraint(m, :((0.1097699546178803 + 0.5720816873107069*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + 0.11808627755560508*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.12730937872071246*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.6166642606025916*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + 0.45996993866072433*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)*(1 + erf((0.1097699546178803 + 0.5720816873107069*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + 0.11808627755560508*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.12730937872071246*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.6166642606025916*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + 0.45996993866072433*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.8288330741680832 + -0.5580861219922526*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + -0.8786358511296881*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.6915444540383735*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.5554450106282389*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + -0.20937059223327026*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.8288330741680832 + -0.5580861219922526*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + -0.8786358511296881*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.6915444540383735*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.5554450106282389*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + -0.20937059223327026*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (0.13615042538554878 + 0.08502572267265318*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + -0.9980047609363689*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.5678523898140919*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.29618367003382495*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + -0.7256587322261439*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)*(1 + erf((0.13615042538554878 + 0.08502572267265318*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + -0.9980047609363689*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.5678523898140919*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.29618367003382495*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + -0.7256587322261439*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.16139551949919406 + -0.5886735107938796*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + -0.02630736114123211*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.03441572289730743*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.3935924617025486*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + 0.3425056157482871*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.16139551949919406 + -0.5886735107938796*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + -0.02630736114123211*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.03441572289730743*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + 0.3935924617025486*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + 0.3425056157482871*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.5329713534130147 + -0.3183971019485696*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + 0.5847814320797045*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.7330877632484647*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + -0.37731199331975374*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + -0.08817423127621415*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.5329713534130147 + -0.3183971019485696*(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))*(1 + erf((-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))/sqrt(2)))/2 + 0.5847814320797045*(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))*(1 + erf((0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))/sqrt(2)))/2 + -0.7330877632484647*(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))*(1 + erf((0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))/sqrt(2)))/2 + -0.37731199331975374*(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))*(1 + erf((0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))/sqrt(2)))/2 + -0.08817423127621415*(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))*(1 + erf((-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using FeedbackParticleFilters, Distributions, Random, PDMats

println("Testing DifferentialRKHSMethodS1.jl:")

@testset "DifferentialRKHSMethodS1.jl" begin
    print("  solver")
    N = 6
    x_pf=[-2.422480820086937 -0.592332203303167 -2.017301296096984 -1.5151245392598531 0.02565906919199346 0.15161614796874012;]
    testens=UnweightedParticleEnsemble(x_pf)
    eq = PoissonEquation(x->[x[1]^2,exp(x[1])], testens)
    update!(eq, testens)
    solve!(eq, DifferentialRKHSMethodS1(1E1, 1E-6));
    print("-")
    @test maximum(abs.(eq.gain[1,:,1] - [ -0.03630508403333434
                                          -0.5438751742317391
                                          -1.0013005320529875
                                          -0.8498038869177289
                                          -1.0778306319564457
                                           1.0096426629620914 ])) < 1e-6
    print("-")
    @test maximum(abs.(eq.gain[1,:,2] - [ -0.04446454607643974
                                           0.15511332541765455
                                           0.0646397390242831
                                           0.10240796498343294
                                           0.35929499989174624
                                          -0.2929181164799246 ])) < 1e-6
    print("-")
    @test maximum(abs.(eq.potential[1,:] - [ 96.14908939499902
                                           -140.78275506175999
                                             56.7520800164404
                                            122.01620383184002
                                            -66.39285991954131
                                            -67.74175826197815 ])) < 1e-6
    print("-")
    @test maximum(abs.(eq.potential[2,:] - [ -17.371718290088566
                                              24.208561463587788
                                              -9.588138488155298
                                             -23.06405207506813
                                              12.786152526976196
                                              13.02919486274801 ])) < 1e-6
    println("DONE")
end#DifferentialRKHSMethodS1.jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
type Hand
  cards::Array{AbstractString}

  function Hand(h::Array)
    @assert length(unique(h)) == length(h) """a hand cannot contain duplicate cards"""
    @assert length(setdiff(h, DECK)) == 0 """hand contains invalid cards"""
    return orderhand!(new(h))
  end
end

function Hand(h::AbstractString)
  return Hand(split(h, ","))
end

function Hand(h::Integer)
  idx = find(n -> in(n, keys(factor(h))), primes(239))
  return Hand(DECK[idx])
end

function Base.convert(::Type{Integer}, h::Hand)
  # first 52 primes
  p = primes(239)
  idx = find(card -> in(card, h.cards), DECK)
  return prod(p[idx])
end

Base.convert(::Type{AbstractString}, h::Hand) = join(h.cards, ",")

Base.hash(hand::Hand) = hash(hand.cards)

Base.isequal{T<:Hand}(x::T, y::T) = hash(x) == hash(y)

Base.(:(==)){T<:Hand}(x::T, y::T) = hash(x) == hash(y)


text:
"""
	Newton
"""
function newton(res::Function, x;
		tol_r = 1.0e-8, tol_d = 1.0e-6)
	y = copy(x)
	Δy = copy(x)

    r = res(y)

    iter = 0

    while norm(r, 2) > tol_r && iter < 25
        ∇r = ForwardDiff.jacobian(res, y)

		try
        	Δy = -1.0 * ∇r \ r
		catch
			@warn "implicit-function failure"
			return y
		end

        α = 1.0

		iter_ls = 0
        while α > 1.0e-8 && iter_ls < 10
            ŷ = y + α * Δy
            r̂ = res(ŷ)

            if norm(r̂) < norm(r)
                y = ŷ
                r = r̂
                break
            else
                α *= 0.5
				iter_ls += 1
            end

			iter_ls == 10 && (@warn "line search failed")
        end

        iter += 1
    end

	iter == 10 && (@warn "Newton failure")

    return y
end


text:
using JuMP
using Gurobi
using Plasmo
using PlasmoAlgorithms

mp = Model(solver = GurobiSolver())
sp1 = Model(solver = GurobiSolver())
sp2 = Model(solver = GurobiSolver())
sp3 = Model(solver = GurobiSolver())

@variable(mp,x[1:2]>=0)
@variable(sp1,y[1:2]>=0)
@variable(sp2,z[1:2]>=0)
@variable(sp3,t[1:2]>=0)
@constraint(mp,x[1]+x[2]>=27)
@objective(mp,Min, 3x[1]+2x[2])

@variable(sp1,x[1:2]>=0)
@constraint(sp1,x[1]+y[2]>=34)
@constraint(sp1,x[2]+y[1]>=52)
@objective(sp1,Min, 4y[1]+6y[2])

@variable(sp2,x[1:2]>=0)
@constraint(sp2,x[1]+z[2]>=61)
@constraint(sp2,x[2]+z[1]>=29)
@objective(sp2,Min, 7z[1]+8z[2])

@variable(sp3,x[1:2]>=0)
@constraint(sp3,x[1]+t[2]>=43)
@constraint(sp3,x[2]+t[1]>=75)
@objective(sp3,Min, 10t[1]+5t[2])

g = ModelGraph()
setsolver(g, GurobiSolver())
n1 = add_node(g)
n2 = add_node(g)
n3 = add_node(g)
n4 = add_node(g)

setmodel(n1, mp)
setmodel(n2, sp1)
setmodel(n3, sp2)
setmodel(n4, sp3)

edge1 = Plasmo.add_edge(g,n1,n2)
edge2 = Plasmo.add_edge(g,n1,n3)
edge3 = Plasmo.add_edge(g,n1,n4)

@linkconstraint(g,[i in 1:2], n1[:x][i] == n2[:x][i])
@linkconstraint(g,[i in 1:2], n1[:x][i] == n3[:x][i])
@linkconstraint(g,[i in 1:2], n1[:x][i] == n4[:x][i])

bendersolve(g,max_iterations = 10)


text:
"""
    $(TYPEDEF)
Simple file data store with a specified data type and a specified key type.
"""
struct FileDataStore <: AbstractDataStore{UInt8}
    folder::String
    FileDataStore(folder::String) = begin
        if !isdir(folder)
            @warn "Folder '$folder' doesn't exist - creating."
            mkpath(folder)
        end
        return new(folder)
    end
end


text:
# SPDX-FileCopyrightText: 2021 Johannes Hampp <johannes.hampp@zeu.jlug.de>
#
# SPDX-License-Identifier: MIT

# Activate julia environment
using Pkg
Pkg.activate("./workflow/envs")
using GlobalEnergyGIS

scenario = snakemake.wildcards["ssp_scenario"]
ssp_year = parse(Int64, snakemake.wildcards["ssp_year"])

# Prepare the scenario data
create_scenario_datasets(scenario, ssp_year)

max_stars_repo_name: euronion/synde

max_issues_repo_name: euronion/synde

max_forks_repo_name: euronion/synde


text:
# ====================================================================
# Component structure
mutable struct OffsetSlope_1D <: AbstractComponent
    offset::Parameter
    x0::Parameter
    slope::Parameter

    function OffsetSlope_1D(offset::Number, x0::Number, slope::Number)
        out = new(Parameter(offset), Parameter(x0), Parameter(slope))
        out.x0.fixed = true
        return out
    end
end

mutable struct OffsetSlope_2D <: AbstractComponent
    offset::Parameter
    x0::Parameter
    y0::Parameter
    slopeX::Parameter
    slopeY::Parameter

    function OffsetSlope_2D(offset::Number, x0::Number, y0::Number, slopeX::Number, slopeY::Number)
        out = new(Parameter(offset), Parameter(x0), Parameter(y0), Parameter(slopeX), Parameter(slopeY))
        out.x0.fixed = true
        out.y0.fixed = true
        return out
    end
end

OffsetSlope(offset, x0, slope) = OffsetSlope_1D(offset, x0, slope)
OffsetSlope(offset, x0, y0, slopeX, slopeY) = OffsetSlope_2D(offset, x0, y0, slopeX, slopeY)

# ====================================================================
# Evaluate component
function evaluate!(buffer, comp::OffsetSlope_1D, x::AbstractDomain{1},
                   offset, x0, slope)
    @. (buffer = slope * (x - x0) + offset)
end

function evaluate!(buffer, comp::OffsetSlope_2D, domain::AbstractDomain{2},
                   offset, x0, y0, slopeX, slopeY)
    x = domain[1]
    y = domain[2]
    @. (buffer =
        slopeX * (x - x0) +
        slopeY * (y - y0) +
        offset)
end


text:
export parse_header_value

const NOOP = (args...) -> nothing
const IDENTITY = (v, args...) -> v

const Param = Pair{AbstractString,AbstractString}
const Params = Vector{Param}
const Entry = Pair{AbstractString,Params}
const Entries = Vector{Entry}

const EOF = nothing

function parse_header_value(header::AbstractString; headername::AbstractString="", parameterized=true, multi=true, paramsep=UInt8(';'), permittedchars=TOKEN, allowempty=false, filterempty=false)
    paramsep = UInt8(paramsep)

    entries = Entries()
    header = string(header)
    value = ""
    paramkey = ""
    paramval = ""
    params = Params()

    buff = Vector{UInt8}()

    lws = true
    tws = false
    inquote = false
    inescape = false
    parsingparamkey = false
    parsingparamval = false
    maxind = length(header)+1
    for i in 1:maxind
        char = i == maxind ? EOF : UInt8(header[i])
        if (lws || tws) && in(char, LWS)
            continue
        elseif inquote && char != EOF
            if inescape
                # TODO: not any character can be quoted...
                push!(buffer, char)
                inescape = false
            elseif b == QUOTE
                inquote = false
                tws = true
            elseif b == BACKSLASH
                inescape = true
            elseif in(b, TEXT)
                push!(buffer, char)
            else
                throw(HTTPError(400; message="Invalid $headername header (charcode $char not allowed in quoted-string)"))
            end
        elseif char == EOF || (multi && char == COMMA)
            if parsingparamval
                paramval = String(buff)
                push!(params, paramkey => paramval)
            elseif parsingparamkey
                paramkey = String(buff)
                push!(params, paramkey => paramval)
            else
                value = String(buff)
            end

            add = true
            if length(value) == 0
                if !allowempty
                    throw(HTTPError(400; message="Invalid $headername header (empty value)"))
                elseif filterempty
                    add = false
                end
            end

            if add
                push!(entries, value => params)
            end
            buff = Vector{UInt8}()
            params = Params()
            value = ""
            paramkey = ""
            paramval = ""
            parsingparamkey = false
            parsingparamval = false
            tws = false
            lws = true
        elseif parameterized && char == paramsep
            if parsingparamval
                paramval = String(buff)
                push!(params, paramkey => paramval)
            elseif parsingparamkey
                paramkey = String(buff)
                push!(params, paramkey => paramval)
            else
                value = String(buff)
            end
            if length(value) == 0 && !allowempty
                throw(HTTPError(400; message="Invalid $headername header (empty value)"))
            end
            buff = Vector{UInt8}()
            parsingparamkey = true
            parsingparamval = false
            paramkey = ""
            paramval = ""
            lws = true
            tws = false
        elseif parsingparamkey
            if char == EQUAL
                paramkey = String(buff)
                buff = Vector{UInt8}()
                parsingparamkey = false
                parsingparamval = true
                tws = false
                lws = true
            else
                if tws
                    throw(HTTPError(400; message="Invalid $headername header (invalid parameter key)"))
                elseif in(char, TOKEN)
                    push!(buff, char)
                elseif in(char, LWS)
                    tws = true
                else
                    throw(HTTPError(400; message="Invalid $headername header (invalid parameter key)"))
                end
                lws=false
            end
        elseif parsingparamval
            if lws && char == QUOTE
                inquote = true
            elseif tws
                throw(HTTPError(400; message="Invalid $headername header (invalid parameter value)"))
            elseif in(char, TOKEN)
                push!(buff, char)
            elseif in(char, LWS)
                tws = true
            else
                throw(HTTPError(400; message="Invalid $headername header (invalid parameter value)"))
            end
            lws=false
        elseif tws
            throw(HTTPError(400; message="Invalid $headername header (invalid value)"))
        elseif in(char, permittedchars)
            push!(buff, char)
            lws = false
            tws = false
        elseif in(char, LWS)
            tws = true
        else
            throw(HTTPError(400; message="Invalid $headername header (unpermitted bytecode)"))
        end
    end

    entries
end


text:
module PlanetaryEphemeris

# __precompile__(false)

export au, yr, sundofs, earthdofs,
    c_au_per_day, μ, NBP_pN_A_J23E_J23M_J2S!,
    NBP_pN_A_J23E_J23M_J2S_threads!, DE430!,
    semimajoraxis, eccentricity, inclination,
    longascnode, argperi, longperi,
    trueanomaly, ecanomaly, meananomaly,
    timeperipass, lrlvec, eccentricanomaly,
    meanan2truean, meanmotion, time2truean,
    su, ea, mo, au, yr, daysec, clightkms,
    c_au_per_day, c_au_per_sec, c_cm_per_sec,
    J2000, R_sun, α_p_sun, δ_p_sun, au,
    UJ_interaction, de430_343ast_ids, Rx, Ry, Rz,
    ITM_und, ITM1, ITM2, R_moon, τ_M, k_2M,
    JSEM, CM, SM, n1SEM, n2M, J2E, J2EDOT, RE,
    k_20E, k_21E, k_22E, τ_0p, τ_1p, τ_2p, τ_0, τ_1, τ_2, ω_E, EMRAT,
    TaylorInterpolant, selecteph2jld, ssb_posvel_pN, nbodyind

using AutoHashEquals
using TaylorIntegration, LinearAlgebra
using Printf
using Dates: DateTime, julian2datetime, datetime2julian, year
using DelimitedFiles
using JLD
using Quadmath

import Base.reverse

# integration parameters
const order = 30
const abstol = 1.0E-20

nbodyind(N::Int, i::Int) = union(3i-2:3i, 3*(N+i)-2:3*(N+i))

function nbodyind(N::Int, ivec::AbstractVector{Int})
    a = Int[]
    for i in ivec
        i > N && continue
        a = union(a, nbodyind(N,i))
    end
    return sort(a)
end

include("constants.jl")
include("jpl-de-430-431-earth-orientation-model.jl")
include("initial_conditions.jl")
include("dynamical_model.jl")
include("interpolation.jl")
include("plephinteg.jl")
include("propagation.jl")
include("osculating.jl")
include("barycenter.jl")

end # module


text:
alg_order(alg::MIRK4) = 4
alg_order(alg::GeneralMIRK4) = 4


text:
using RecipesBase
using ColorTypes
import PlotUtils: cgrad, palette, color_list
using LaTeXStrings

@userplot NozzlePlot

@recipe function f(h::NozzlePlot;fields=(),gas=DefaultPerfectGas)
  nfields = length(fields) + 1

  if nfields == 1
    length(h.args) == 1 || error("`nozzleplot` should be given one argument.  Got: $(typeof(h.args))")
    noz, = h.args
  else
    length(h.args) == 4 || error("`nozzleplot` should be given four arguments.  Got: $(typeof(h.args))")
    noz, pb, p0, T0 = h.args
    nozproc = NozzleProcess(noz,pb,p0,T0,gas=gas)
  end

  linecolor := 1
  if plotattributes[:plot_object].n > 0
    lc = plotattributes[:plot_object].series_list[end][:linecolor]
    idx = findall(x -> RGBA(x) == lc,color_list(palette(:default)))
    linecolor := idx[1] + 1
  end

  layout := (nfields,1)
  size := (500,200*nfields)
  xticks := (0:0,["Throat"])

  xline = [0,0]

  fields_lower = lowercase.(fields)

  subnum = 0
  if in("pressure",fields_lower)
    subnum += 1

    p = ustrip.(pressure(nozproc),u"kPa")
    pmax = maximum(p)+50
    @series begin
      subplot := subnum
      linestyle --> :dash
      linecolor := :black
      #label := ""
      xline, [0,pmax]
    end
    @series begin
      subplot := subnum
      ylims --> (0,pmax)
      xlims := (-Inf,Inf)
      yguide := "Pressure (KPa)"
      #legend := true
      #annotations := (positions(noz)[end],p[end],nozzle_quality(noz,pb,p0))
      positions(noz), p
    end
  end

  if in("temperature",fields_lower)
    subnum += 1

    T = ustrip.(temperature(nozproc),u"K")
    Tmax = maximum(T)+100.0
    @series begin
      subplot := subnum
      linestyle --> :dash
      linecolor := :black
      xline, [0,Tmax]
    end
    @series begin
      subplot := subnum

      ylims --> (0,Tmax)
      xlims := (-Inf,Inf)
      yguide := "Temperature (K)"
      positions(noz), T
    end
  end

  if in("density",fields_lower)
    subnum += 1

    ρ = ustrip.(density(nozproc))
    ρmax = maximum(ρ)+1.0
    @series begin
      subplot := subnum
      linestyle --> :dash
      linecolor := :black
      xline, [0,ρmax]
    end
    @series begin
      subplot := subnum

      ylims --> (0,ρmax)
      xlims := (-Inf,Inf)
      yguide := "Density (kg/m3)"
      positions(noz), ρ
    end
  end

  if in("machnumber",fields_lower) || in("mach",fields_lower)
    subnum += 1

    M = ustrip.(machnumber(nozproc))
    Mmax = maximum(M)+1.0
    @series begin
      subplot := subnum
      linestyle --> :dash
      linecolor := :black
      xline, [0,Mmax]
    end
    @series begin
      subplot := subnum
      linestyle --> :dash
      linecolor := :black
      positions(noz), ones(length(positions(noz)))
    end
    @series begin
      subplot := subnum

      ylims --> (0,Mmax)
      xlims := (-Inf,Inf)
      yguide := "Mach Number"
      positions(noz), M
    end
  end

  subnum += 1
  A = ustrip.(areas(noz),u"cm^2")
  Amax = maximum(A)+10
  @series begin
    subplot := subnum
    linestyle --> :dash
    linecolor := :black
    xline, [0,Amax]
  end
  #@series begin
     subplot := subnum

     ylims --> (0,Amax)
     xlims := (-Inf,Inf)

     yguide := "Area (sq cm)"
     positions(noz),A
  #end

end


text:
function preprocess_img(img, shape)
    resized_img = Images.imresize(img, (shape[2], shape[1]))
    sample = channelview(resized_img) * 256
    sample[1,:,:] -= 123.68
    sample[2,:,:] -= 116.779
    sample[3,:,:] -= 103.939
    sample = permutedims(sample, [2,3,1])
    return reshape(sample, (size(sample)[1], size(sample)[2], 3, 1))
end

function postprocess_img(img)
    img = reshape(img, (size(img)[1], size(img)[2], 3))
    img[:,:,1] += 123.68
    img[:,:,2] += 116.779
    img[:,:,3] += 103.939
	img = permutedims(img, [3,1,2])
    img = clamp.(img, 0, 255) / 255
    colorview(RGB{Float32}, img)
end


text:
module TestReplaceCommits

using ReleaseFlow: replace_commits_since_impl, message_readme_change
using Test

mkurl(tag) =
    "https://img.shields.io/github/commits-since/USER/PACKAGE.jl/$tag.svg"

@testset begin
    readme = replace_commits_since_impl((
        tag = "TAG",
        orig = mkurl("v1.2.3"),
    ))
    @test readme.content == mkurl("TAG")

    readme = replace_commits_since_impl((
        tag = "TAG",
        orig = "ignored",
        other = 123,
    ))
    @test readme.content == "ignored"
    @test readme.other == 123

    @info "↓↓↓ (to be ignored)"
    message_readme_change(((
        origmatch = (
            captures = [
                "https://hostname/paths/"
                "..."
                ".svg"
            ],
        ),
        content = mkurl("v1.2.3"),
    )))
    @info "↑↑↑ (to be ignored)"
end

end  # module


text:
using BinaryBuilder

# Collection of sources required to build SHTOOLS
name = "SHTOOLS"
version = v"4.9.1"
sources = [
    ArchiveSource("https://github.com/SHTOOLS/SHTOOLS/releases/download/v4.9.1/SHTOOLS-4.9.1.tar.gz",
                  "5c22064f9daf6e9aa08cace182146993aa6b25a6ea593d92572c59f4013d53c2"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SHTOOLS-*

# Patch source code
# Don't use libtool
atomic_patch -p0 $WORKSPACE/srcdir/patches/no-libtool.patch
# Correct C interface for MakeGradientDH (see <https://github.com/SHTOOLS/SHTOOLS/issues/328>)
atomic_patch -p1 $WORKSPACE/srcdir/patches/correct-cMakeGradientDH.patch

# Build and install static libraries
make fortran -j${nproc} F95FLAGS="-fPIC -O3 -std=gnu"
make fortran-mp -j${nproc} F95FLAGS="-fPIC -O3 -std=gnu"
make install PREFIX=${prefix}

# Create shared libraries
# TODO: Do we need to add " | cut -d' ' -f1"?
whole_archive=$(flagon --whole-archive)
if [ -n "${whole_archive}" ]; then
    whole_archive="-Wl,${whole_archive}"
fi
no_whole_archive=$(flagon --no-whole-archive)
if [ -n "${no_whole_archive}" ]; then
    no_whole_archive="-Wl,${no_whole_archive}"
fi
gfortran -shared -o ${libdir}/libSHTOOLS.${dlext} ${whole_archive} ${prefix}/lib/libSHTOOLS.a ${no_whole_archive} -lfftw3 -lopenblas -lm
gfortran -fopenmp -shared -o ${libdir}/libSHTOOLS-mp.${dlext} ${whole_archive} ${prefix}/lib/libSHTOOLS-mp.a ${no_whole_archive} -lfftw3 -lopenblas -lm
"""

platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libSHTOOLS", :libSHTOOLS),
    LibraryProduct("libSHTOOLS-mp", :libSHTOOLS_mp),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("FFTW_jll"),
    Dependency("OpenBLAS32_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module DiscreteFunctions

using SimpleTools, Permutations, LinearAlgebra

import Base: length, show, getindex, setindex!, *, ==, hash, ^, inv, +, Matrix
import Permutations: fixed_points, Permutation, cycles

import LinearAlgebra: eigvals

export DiscreteFunction, IdentityFunction, RandomFunction, has_inv
export fixed_points, image, is_permutation, Permutation, eigvals

"""
`DiscreteFunction` is a function from `{1,2,...,n}` to itself.
It can be created by `DiscreteFunction(list)` where `list`
is a one-dimensional array of positive integers. Alternatively,
it can be created using positive integer arguments.

The following are equivalent:
```
DiscreteFunction([1,4,2,3])
DiscreteFunction(1,4,2,3)
```
"""
struct DiscreteFunction
    data::Vector{Int}
    function DiscreteFunction(a::Vector{Int})
        if minimum(a)<1
            error("All function values must be positive")
        end
        n = max(length(a), maximum(a))
        data_ = ones(Int,n)
        for j=1:length(a)
            data_[j] = a[j]
        end
        new(data_)
    end
end

function DiscreteFunction(a::Int, args...)
    n = length(args)+1
    data = zeros(Int,n)
    data[1] = a
    for j=1:n-1
        data[j+1] = args[j]
    end
    return DiscreteFunction(data)
end

DiscreteFunction(p::Permutation) = DiscreteFunction(p.data)

function Permutation(f::DiscreteFunction)
    if is_permutation(f)
        return Permutation(f.data)
    end
    error("This function cannot be converted to a Permutation")
end

matrix2vec(A::Matrix) = matrix2vec(A .!= 0)

function matrix2vec(A::BitArray{2})
    r,c = size(A)
    if r != c
        @error "Matrix must be square"
    end
    u = sum(A,dims=2)  # get row sums
    if !all(u.==1)
        @error "Matrix should have a single nonzero entry in each row"
    end

    vec = zeros(Int,r)
    for i=1:r
        vec[i] = first(findall(A[i,:]))
    end

    return vec

end

"""
`DiscreteFunction(A)` creates a new
discrete function based on the square matrix `A`. It is required
that each row of `A` have exactly one nonzero element. Note that
if `A=Matrix(f)` then `DiscreteFunction(A)==f`.
"""
DiscreteFunction(A::Matrix) = DiscreteFunction(matrix2vec(A))
DiscreteFunction(A::BitArray{2}) = DiscreteFunction(matrix2vec(A))

"""
`IdentityFunction(n)` creates the identity `DiscreteFunction` on
the set `{1,2,...,n}`.
"""
function IdentityFunction(n::Int)::DiscreteFunction
    @assert n>0 "Argument must be positive"
    data = collect(1:n)
    return DiscreteFunction(data)
end

"""
`RandomFunction(n)` creates a random `DiscreteFunction`
on `{1,2,...,n}`.
"""
function RandomFunction(n::Int)
    @assert n>0 "Argument must be positive"
    data = rand(1:n,n)
    return DiscreteFunction(data)
end

length(f::DiscreteFunction) = length(f.data)

(f::DiscreteFunction)(x::Int) = f.data[x]
getindex(f::DiscreteFunction, x::Int) = f.data[x]

function setindex!(f::DiscreteFunction, val::Int, x::Int)
    n = length(f)
    if val < 1 || val > n
        error("Invalid value $val for x=$x; must be between 1 and $n")
    end
    f.data[x] = val
end

(==)(f::DiscreteFunction,g::DiscreteFunction) = f.data == g.data

function show(io::IO, f::DiscreteFunction)
    n = length(f)
    print(io,"DiscreteFunction on [$n]\n")
    w = 4  # window size for printing; shame to hard code this
    for i=1:n
        print(io,flush_print(i,w))
    end
    print(io,"\n")
    for i=1:n
        print(io,flush_print(f(i),w))
    end
    print(io,"\n")
end

function (*)(f::DiscreteFunction, g::DiscreteFunction)::DiscreteFunction
    n1 = length(f)
    n2 = length(g)
    if n1 != n2
        error("Cannot compose `DiscreteFunction`s of different lengths")
    end
    data = zeros(Int,n1)
    for j=1:n1
        data[j] = f(g(j))
    end
    return DiscreteFunction(data)
end

hash(f::DiscreteFunction,h::UInt64) = hash(f.data, h)
hash(f::DiscreteFunction) = hash(f.data)

"""
`has_inv(f::DiscreteFunction)` tests if `f` is invertible.
"""
function has_inv(f::DiscreteFunction)::Bool
    n = length(f)
    return collect(1:n) == sort(f.data)
end

function inv(f::DiscreteFunction)::DiscreteFunction
    @assert has_inv(f) "This function is not invertible"
    n = length(f)
    data = zeros(Int,n)
    for j=1:n
        data[f(j)] = j
    end
    return DiscreteFunction(data)
end

function (^)(f::DiscreteFunction, t::Integer)
    if t<0
        return inv(f)^t
    end
    if t==0
        return IdentityFunction(length(f))
    end
    if t==1
        return f
    end
    half_t = Int(floor(t/2))
    g = f^half_t
    if t%2 == 0
        return g*g
    end
    return f*g*g
end

"""
If `f` and `g` are of type `DiscreteFunction`, then `f+g`
is their *disjoint sum*.
"""
function (+)(f::DiscreteFunction, g::DiscreteFunction)::DiscreteFunction
    n1 = length(f)
    d1 = f.data
    d2 = g.data .+ n1
    d = vcat(d1,d2)
    return DiscreteFunction(d)
end

"""
`fixed_points(f::DiscreteFunction)` returns a list of fixed points of
the function, i.e., those values `x` such that `f(x)==x`.
"""
fixed_points(f::DiscreteFunction) = [ i for i in 1:length(f) if f(i)==i ]

"""
`image("f::DiscreteFunction")` returns a `Set` containing all the output
values of `f`.
"""
image(f::DiscreteFunction)::Set{Int} = Set{Int}(f.data)

"""
`is_permutation(f::DiscreteFunction)` returns `true` if `f` is a
bijection on its domain.
"""
is_permutation(f::DiscreteFunction)::Bool = has_inv(f)

"""
`Matrix(f::DiscreteFunction)` returns an `n`-by-`n` zero-one matrix
in which there is a `1` in position `i,j` exactly when `f(i)==j`.
"""
function Matrix(f::DiscreteFunction)::Array{Int,2}
    n = length(f)
    A = zeros(Int,n,n)
    for i=1:n
        j = f(i)
        A[i,j]=1
    end
    return A
end

"""
`eigvals(f::DiscreteFunction)` returns the eigenvalues of `Matrix(f)`.
"""
eigvals(f::DiscreteFunction) = LinearAlgebra.eigvals(Matrix(f))

include("all_functions.jl")
include("cycles.jl")
include("sources.jl")
include("opt_sqrt.jl")

end  # end of module