text:
function update_rewiring_parameters(parameters::Dict{Symbol,Any}, biomass, t)
  S = size(parameters[:A], 1)

  if parameters[:rewire_method] == :ADBM
    # Is there any extinction and/or rewiring happening?
    A = parameters[:A] #interaction matrix
    Anew = ADBM(S,parameters,biomass) #potentially new interaction  matrix
    if A != Anew #different matrices => rewired
      append!(parameters[:tmpA], [parameters[:A]]) #store the old matrix
      if parameters[:adbm_trigger] == :interval
        append!(parameters[:rewiretime], t) #store the time t whichrewiring happened
      end
    end
    #Is there any extinction?
    extinct = findall(biomass .< 100*eps())
    for i in extinct
      if i ∉ parameters[:extinctions]
        append!(parameters[:extinctions], i) ;
        append!(parameters[:extinctionstime], [(t, i)])
      end
    end
    sort!(parameters[:extinctions])

    #assign new array
    parameters[:A] = Anew

    #update the parameters
    get_herbivores(parameters) #
    getW_preference(parameters) #
    get_efficiency(parameters)

  elseif parameters[:rewire_method] ∈ [:Gilljam, :DS]
    #add extinction
    extinct = findall(biomass .< 100*eps())
    for i in extinct
      if i ∉ parameters[:extinctions]
        append!(parameters[:extinctions], i) ;
        append!(parameters[:extinctionstime], [(t, i)])
        append!(parameters[:tmpA], [parameters[:A]])
      end
    end
    sort!(parameters[:extinctions])

    #assign new array and update costs
    parameters[:A] , parameters = Gilljam(S,parameters,biomass)

    #update rewiring parameters
    if parameters[:preferenceMethod] == :specialist
      parameters = update_specialist_preference(parameters,S)
    end

    #update parameters
    get_herbivores(parameters)
    getW_preference(parameters)
    parameters[:w][(LinearIndices(parameters[:w]))[findall(x -> x .== Inf, parameters[:w])]] .= 1
    get_efficiency(parameters)

  elseif parameters[:rewire_method] ∈ [:stan, :DO]
    #add extinction
    extinct = findall(biomass .< 100*eps())
    for i in extinct
      if i ∉ parameters[:extinctions]
        append!(parameters[:extinctions], i) ;
        append!(parameters[:extinctionstime], [(t, i)])
        append!(parameters[:tmpA], [parameters[:A]])
      end
    end
    sort!(parameters[:extinctions])

    parameters[:A] = Staniczenko_rewire(parameters)
    get_herbivores(parameters)
    getW_preference(parameters)
    get_efficiency(parameters)

  else
    #add extinction
    extinct = findall(biomass .< 100*eps())
    for i in extinct
      if i ∉ parameters[:extinctions]
        append!(parameters[:extinctions], i) ;
        append!(parameters[:extinctionstime], [(t, i)])
        append!(parameters[:tmpA], [parameters[:A]])
      end
    end
    sort!(parameters[:extinctions])
  end

  return parameters
end

function update_specialist_preference(parameters::Dict{Symbol,Any})
  S = size(parameters[:A], 1)
  #find those that have lost their specialistPref and are not extinct
  lost = falses(S)
  for i = 1:S
    if !parameters[:is_producer][i] && !in(i,parameters[:extinctions])
      if parameters[:A][i,parameters[:specialistPref][i]] == 0
        #assign new random species from diet
        parameters[:specialistPref][i] = sample(findall(!iszero, parameters[:A][i,:]))
      end
    else
      parameters[:specialistPref][i] = 0
    end
  end
  #return(p_r)
end


text:
function resource_url(app::DashApp, resource::AppRelativeResource)
    prefix = get_setting(app, :requests_pathname_prefix)
    return string(prefix,
         "_dash-component-suites/",
         resource.namespace,
         "/",
         build_fingerprint(resource.relative_path, resource.version, resource.ts)
     )
end

#TODO move to somewhere
function asset_path(app::DashApp, path::AbstractString)
    return string(
        get_setting(app, :requests_pathname_prefix),
        lstrip(get_setting(app, :assets_url_path), '/'),
        "/",
        path
    )
end

resource_url(app::DashApp, resource::AppExternalResource) = resource.url
function resource_url(app::DashApp, resource::AppAssetResource)
    return string(
        asset_path(app, resource.path),
        "?m=", resource.ts
    )
end

make_css_tag(url::String) = """<link rel="stylesheet" href="$(url)">"""
make_css_tag(app::DashApp, resource::AppCustomResource) = format_tag("link", resource.params, opened = true)

make_css_tag(app::DashApp, resource::AppResource) = make_css_tag(resource_url(app, resource))

make_script_tag(url::String) = """"""
make_script_tag(app::DashApp, resource::AppCustomResource) = format_tag("script", resource.params)
make_script_tag(app::DashApp, resource::AppResource) = make_script_tag(resource_url(app, resource))
make_inline_script_tag(script::String) = """"""

function metas_html(app::DashApp)
    meta_tags = app.config.meta_tags
    has_ie_compat = any(meta_tags) do tag
        get(tag, "http-equiv", "") == "X-UA-Compatible"
    end
    has_charset = any(tag -> haskey(tag, "charset"), meta_tags)

    result = String[]
    !has_ie_compat && push!(result, "<meta http-equiv=\"X-UA-Compatible\" content=\"IE=edge\">")
    !has_charset && push!(result, "<meta charset=\"UTF-8\">")

    append!(result, format_tag.("meta", meta_tags, opened = true))
    return join(result, "\n        ")

end

function css_html(app::DashApp, resources::ApplicationResources)
     join(
         make_css_tag.(Ref(app), resources.css), "\n       "
         )
end

function scripts_html(app::DashApp, resources::ApplicationResources)
    include_string = join(
         vcat(
            make_script_tag.(Ref(app), resources.js),
            make_inline_script_tag.(app.inline_scripts)
        ),
         "\n       "
        )

end

app_entry_html() = """

        Loading...

"""

function config_html(app::DashApp)
    config = Dict{Symbol, Any}(
        :url_base_pathname => get_setting(app, :url_base_pathname),
        :requests_pathname_prefix => get_setting(app, :requests_pathname_prefix),
        :ui => get_devsetting(app, :ui),
        :props_check => get_devsetting(app, :props_check),
        :show_undo_redo => get_setting(app, :show_undo_redo),
        :suppress_callback_exceptions => get_setting(app, :suppress_callback_exceptions),
        :update_title => get_setting(app, :update_title)
    )
    if get_devsetting(app, :hot_reload)
        config[:hot_reload] = (
            interval = trunc(Int, get_devsetting(app, :hot_reload_interval) * 1000),
            max_retry = get_devsetting(app, :hot_reload_max_retry)
        )
    end
    return """"""
end

renderer_html() = """"""

favicon_html(app::DashApp) = ""

function index_page(app::DashApp, resources::ApplicationResources)

    result = interpolate_string(app.index_string,
        metas = metas_html(app),
        title = app.title,
        favicon = favicon_html(app),
        css = css_html(app, resources),
        app_entry = app_entry_html(),
        config = config_html(app),
        scripts = scripts_html(app, resources),
        renderer = renderer_html()
    )

    validate_index("index", result,
        [
            "#react-entry-point" => r"id=\"react-entry-point\"",
            "#_dash_config" => r"id=\"_dash-config\"",
            "dash-renderer" => r"src=\"[^\"]*dash[-_]renderer[^\"]*\"",
            "new DashRenderer" => r"id=\"_dash-renderer",
        ]
    )
    return result
end

ext: jl


text:
if false                       # if 从句需要 Bool 测试
  z = 1
elseif 1 == 2
  z = 2
else
  a = 3
end                            # 此后 a = 3，z 未定义

1==2 ? "A" : "B"               # 标准三元操作符

i = 1
while true
  i += 1
  if i > 10
    break
  end
end

for x in 1:10                  # x 在区间中，这里可以用 = 代替 in
  if 3 < x < 6
    continue                   # 跳过一次迭代
  end
  println(x)
end                            # x 在此循环的内部作用域被定义

[10i+j for i in 1:3, j in 1:3] # 隐式循环构造矩阵

x = (a = 1; 2 * a)             # 用 ; 合并多行语句此后  x = 2; a = 1
y = begin                      # 也可以用 begin...end 来封装多个语句
  b = 3                        # 多个语句的总体结果是最后一个语句的值
  3 * b
end                            # 此后 y = 9; b = 3

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
using FastHenryHelper: Context, xyz1, nodes_xyz1, width_height,
                        wxyz, corners_xyz1_thick

function contexttest1()
  (g1,(t,sp,pn1,pn2,u1,u2,u3,u4,u5,u6,u7,n1,n2,n3,n4,n5,n6,n7,d1,p1,s1)) =
    groupfortests()

  context = Context(g1)

  #@test context.dict[t].default == Default()
  #@test context.dict[t].units == Units()
  #@test context.dict[t].autoname == 0
  @test context.firstunits == Units("km")
  @test context.title == "Test Title"

  @test context.dict[n1].default == Default()
  @test context.dict[n1].units == Units("km")
  @test context.dict[n1].autoname == 1

  @test context.dict[n2].default == Default()
  @test context.dict[n2].units == Units("m")
  @test context.dict[n2].autoname == 2

  @test context.dict[n3].default == Default()
  @test context.dict[n3].units == Units("cm")
  @test context.dict[n3].autoname == 2

  @test context.dict[n4].default == Default()
  @test context.dict[n4].units == Units("mm")
  @test context.dict[n4].autoname == 2

  @test context.dict[n5].default == Default(w=3)
  @test context.dict[n5].units == Units("um")
  @test context.dict[n5].autoname == 3

  @test context.dict[n6].default == Default(w=3)
  @test context.dict[n6].units == Units("in")
  @test context.dict[n6].autoname == 4

  @test context.dict[n7].default == Default(w=3)
  @test context.dict[n7].units == Units("mils")
  @test context.dict[n7].autoname == 5

  @test context.dict[p1].default == Default(w=3)
  @test context.dict[p1].units == Units("mils")
  @test context.dict[p1].autoname == 6

  @test context.dict[s1].default == Default(w=3)
  @test context.dict[s1].units == Units("mils")
  @test context.dict[s1].autoname == 7

  @test xyz1(n1, context) ≈ [0.0, 0.0, 1.0, 1.0]
  @test xyz1(n2,context) ≈ [0.0,0.0, 0.001  ,1.0]
  @test xyz1(n3,context) ≈ [0.0,0.0, 1.0e-5 ,1.0]
  @test xyz1(n4,context) ≈ [0.0,0.0, 1.0e-6 ,1.0]
  @test xyz1(n5,context) ≈ [0.0,0.0, 1.0e-9 ,1.0]
  @test xyz1(n6,context) ≈ [0.0,0.0, 2.54e-5,1.0]
  @test xyz1(n7,context) ≈ [0.0,0.0, 2.54e-8,1.0]

  (n1xyz1,n2xyz1) = nodes_xyz1(s1,context)
  @test n1xyz1 ≈ [0.0,0.0,1.0e-5,1.0]
  @test n2xyz1 ≈ [0.0,0.0,1.0e-6,1.0]

  @test [width_height(s1,context)...] ≈ [3.0e-9,5.08e-8]

  @test wxyz(s1,context) ≈ [-1.0,0.0,0.0]

  (c1,c2,c3,thick) = corners_xyz1_thick(p1,context)
  @test c1 ≈ [0.0,0.0,0.0,1.0]
  @test c2 ≈ [2.54e-7,0.0,0.0,1.0]
  @test c3 ≈ [2.54e-7,2.54e-7,0.0,1.0]
  @test thick ≈ 2.54e-8

  (pn1xyz1,pn2xyz1)=nodes_xyz1(p1,context)
  @test pn1xyz1 ≈ [2.54e-8,5.08e-8,0.0,1.0]
  @test pn2xyz1 ≈ [1.016e-7,1.27e-7,0.0,1.0]
end
contexttest1()

function contexttest2()
  n1 = Node(1,2,3)
  n2 = Node(4,5,6)
  sp = SegmentParameters(w=1,h=2)
  d = Default(w=7,h=8)
  s1 = Segment(n1,n2,sp)
  s2 = Segment(n1,n2)
  mm = Units("mm")
  cm  = Units("cm")

  g2 = Group([])
  context = Context(g2)

  g3 = Group([mm,s1])
  context = Context(g3)
  (n1xyz1,n2xyz1) = nodes_xyz1(s1,context)
  @test n1xyz1 ≈ [1.0,2.0,3.0,1.0]
  @test n2xyz1 ≈ [4.0,5.0,6.0,1.0]
  @test [width_height(s1,context)...] ≈ [1.0,2.0]

  g4 = Group([cm,n1,n2,d,mm,s2])
  context = Context(g4)
  (n1xyz1,n2xyz1) = nodes_xyz1(s2,context)
  @test n1xyz1 ≈ [1.0,2.0,3.0,1.0]
  @test n2xyz1 ≈ [4.0,5.0,6.0,1.0]
  @test [width_height(s2,context)...] ≈ [7.0,8.0]

  g5 = Group([cm,n1,n2,mm,d,s2])
  context = Context(g5)
  (n1xyz1,n2xyz1) = nodes_xyz1(s2,context)
  @test n1xyz1 ≈ [1.0,2.0,3.0,1.0]
  @test n2xyz1 ≈ [4.0,5.0,6.0,1.0]
  @test [width_height(s2,context)...] ≈ [0.7,0.8]

  g6 = Group([cm,n1,mm,n2,d,s2])
  context = Context(g6)
  (n1xyz1,n2xyz1) = nodes_xyz1(s2,context)
  @test n1xyz1 ≈ [1.0,2.0,3.0,1.0]
  @test n2xyz1 ≈ [0.4,0.5,0.6,1.0]
  @test [width_height(s2,context)...] ≈ [0.7,0.8]

  g7 = Group([cm,mm,n1,n2,d,s2])
  context = Context(g7)
  (n1xyz1,n2xyz1) = nodes_xyz1(s2,context)
  @test n1xyz1 ≈ [0.1,0.2,0.3,1.0]
  @test n2xyz1 ≈ [0.4,0.5,0.6,1.0]
  @test [width_height(s2,context)...] ≈ [0.7,0.8]
end
contexttest2()

function testrecursioncheck()
  g = Group([Node(:a,2,3,4)])
  push!(g,g)
  @test_throws(ErrorException,print(g))
end
testrecursioncheck()


text:
using GrayCoding
using Test

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


text:
function square_distance(coords1, coords2)
    """Compute the square distances between coords1 and coords2.

    Note: Depending on the input value and the dtype, negative values may
    be mixed in the return value.

    Args:
        coords1: coordinates (N, C)
        coords2: coordinates (M, C)

    Returns:
        square distances:
            square distances between coords1 and coords2 (N, M)
    """
    dot_product = -2 * coords1 * coords2'
    column = repeat(sum(coords1 .* coords1, dims=2), 1, size(coords2, 1))
    row = repeat(sum(coords2 .* coords2, dims=2), 1, size(coords1, 1))'
    square_dist = column + dot_product + row
    return square_dist
end

if abspath(PROGRAM_FILE) == @__FILE__
    arr1 = reshape(range(0, 23, step=1), 8, 3)
    arr2 = reshape(range(24, 47, step=1), 8, 3)
    @time square_distance(arr1, arr2)
    @time square_distance(arr1, arr2)
end


text:
using PackageCompilerX
using Test
using Libdl

@testset "PackageCompilerX.jl" begin
    # Write your own tests here.
    sysimage_path = "sys." * Libdl.dlext
    PackageCompilerX.create_sysimage(:Example; sysimage_path=sysimage_path,
                                     precompile_execution_file="precompile_execution.jl",
                                     precompile_statements_file="precompile_statements.jl")
    run(`$(Base.julia_cmd()) -J $(sysimage_path) -e 'println(1337)'`)

    # Test creating an app without bundling
    app_dir = joinpath(@__DIR__, "..", "examples/MyApp/")
    PackageCompilerX.create_app(app_dir)
    app_path = abspath(app_dir, "MyApp" * (Sys.iswindows() ? ".exe" : ""))
    app_output = read(`$app_path`, String)
    @test occursin("Example.domath", app_output)
end


text:
workspace()
using HoTRG_lemon.SimulatorModule
using HoTRG_lemon.SpinModule
using HoTRG_lemon.LatticeModule
using TensorOperations
using TensorMatrices_lemon

dimM = 5
trotterparameter = 0.01
trotteriteration = 30

isinginfo = SpinInfo("quantum_ising_2",1.0e-13,0.0)
trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
ising = QuantumIsingModel(isinginfo, trotterinfo)
lattice = buildLattice("quantum_2d_fractal", ising)
simulator = buildSimulator(lattice,dimM)

# fieldrange = linspace(0.1, 3.5, 10)
# fieldrange = 0.1
print("Input Float for fieldrange: ")
fieldrange = parse(Float64, readline())
simulatorQuantum(fieldrange, simulator; filename = "log_mag.txt", printlog="mag")


text:
abstract type SafeInteger  <: Integer     end
abstract type SafeUnsigned <: SafeInteger end
abstract type SafeSigned   <: SafeInteger end

primitive type SafeInt8    <: SafeSigned     8 end
primitive type SafeInt16   <: SafeSigned    16 end
primitive type SafeInt32   <: SafeSigned    32 end
primitive type SafeInt64   <: SafeSigned    64 end
primitive type SafeInt128  <: SafeSigned   128 end

primitive type SafeUInt8   <: SafeUnsigned   8 end
primitive type SafeUInt16  <: SafeUnsigned  16 end
primitive type SafeUInt32  <: SafeUnsigned  32 end
primitive type SafeUInt64  <: SafeUnsigned  64 end
primitive type SafeUInt128 <: SafeUnsigned 128 end

if Sys.WORD_SIZE == 32
    const SafeInt  = SafeInt32
    const SafeUInt = SafeUInt32
else
    const SafeInt  = SafeInt64
    const SafeUInt = SafeUInt64
end

const UnsafeSigned = Union{Int128, Int64, Int32, Int16, Int8}
const UnsafeUnsigned = Union{UInt128, UInt64, UInt32, UInt16, UInt8}
const UnsafeInteger = Union{UnsafeSigned, UnsafeUnsigned}


text:
function __init__()
    ScientificTypes.set_convention(MLJ())
    ScientificTypes.TRAIT_FUNCTION_GIVEN_NAME[:table] = Tables.istable
end


text:
function goToBot!(r)
    moves = []
    pushMoves!(r, moves)
    while moves[end] != 0 && moves[end-1] != 0
        pushMoves!(r, moves)
    end
    return moves
end

function pushMoves!(r, moves)
    push!(moves, goToWall!(r, Sud))
    push!(moves, goToWall!(r, West))
end

function goToWall!(r, side)
    k = 0
    while !isborder(r, side)
        k+=1
        move!(r, side)
    end
    return k
end

function setMarkers(r)
    sides = [Nord, Ost, Sud, West]
    for side = sides
        putmarker!(r)
        goToAllLine!(r, side)
    end
end

function goToAllLine!(r, side)
    while !isborder(r, side)
        move!(r, side)
    end
end

function goToStart!(r, moves)
    reverse!(moves)
    for n = 1:length(moves)
        if n%2==1
            moves!(r, Ost, moves[n])
        else
            moves!(r, Nord, moves[n])
        end
    end
end

function moves!(r, side, move)
    for _ = 1:move
        move!(r, side)
    end
end

function drawCorners(r)
    moves = goToBot!(r)

    setMarkers(r)

    goToStart!(r, moves)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
################
# License Type #
################

@ghdef mutable struct License
    key::Union{String, Nothing}
    name::Union{String, Nothing}
    spdx_id::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    description::Union{String, Nothing}
    permissions::Union{Vector{String}, Nothing}
    conditions::Union{Vector{String}, Nothing}
    limitations::Union{Vector{String}, Nothing}
    body::Union{String, Nothing}
    featured::Union{Bool, Nothing}
end

License(spdx_id::AbstractString) = License(Dict("spdx_id" => spdx_id))

namefield(license::License) = license.spdx_id

###############
# API Methods #
###############

@api_default function licenses(api::GitHubAPI; options...)
    results, page_data = gh_get_paged_json(api, "/licenses"; options...)
    results = [License(result) for result in results]
    return results, page_data
end

@api_default function license(api::GitHubAPI, license_id; options...)
    result = gh_get_json(api, "/licenses/$license_id"; options...)
    return License(result)
end

@api_default function repo_license(api::GitHubAPI, repo_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo_obj))/license"; options...)
    return Content(result)
end


text:
using DSGE, ModelConstructors, Distributed
using Nullables, DataFrames, OrderedCollections, Dates

#############################
# To use:
# Just run in the Julia REPL
# include("run_default.jl")
# Note that the estimation
# step will take 2-3 hours.
#############################

##############
# Model Setup
##############
# Instantiate the FRBNY DSGE model object
m = Model1002()

# estimate as of 2015-Q3 using the default data vintage from 2015 Nov 27
m <= Setting(:data_vintage, "181115")
m <= Setting(:date_forecast_start, quartertodate("2018-Q4"))

# The following settings ensure that this script runs in
# a short amount of time. To properly estimate and
# forecast, we recommend either using the default settings
# (i.e. comment out the settings below) or
# changing the settings yourself.
m <= Setting(:n_mh_simulations, 100) # Do 100 MH steps during estimation
m <= Setting(:use_population_forecast, false) # Population forecast not available as data to turn off
m <= Setting(:forecast_block_size, 5) # adjust block size to run on small number of estimations

#############
# Estimation
#############
# Reoptimize parameter vector, compute Hessian at mode, and full posterior
# parameter sampling.
#
# Note some columns will have missing data because not all our data
# is publicly available. By default, `load_data` will error
# if any of the columns in the loaded DataFrame is empty,
# but we turn this feature off by setting the keyword `check_empty_columns = false`.
# Warnings will be still thrown after calling load_data indicating which columns
# are empty. However, estimate will still run when data is missing.
@time begin
    df = load_data(m, try_disk = false, check_empty_columns = false, summary_statistics = :none)
    data = df_to_matrix(m, df)
    estimate(m, data)
end

# produce LaTeX tables of parameter moments
moment_tables(m)

############
# Forecasts
############
# (1) forecast_one produces forecasts of the variables specified in output_vars
# and stores the forecasts in model units (log deviations from steady state).
# (2) compute_meansbands stores the forecast in a MeansBands type, defined in DSGE.jl,
# which transforms the forecast from model units to "human readable" units, i.e.
# the GDP forecast is transformed from annualized quarterly log per-capita growth rates
# to annualized quarterly aggregate percent change.

# The variables we want to forecast. In this case, all of the model observables
output_vars = [:histobs, :forecastobs]

# Modal forecast (point forecast)
forecast_one(m, :mode, :none, output_vars; check_empty_columns = false)
compute_meansbands(m, :mode, :none, output_vars; check_empty_columns = false)

# Full-distribution forecast (point forecast (mean) and uncertainty bands)
# Optionally add 10 processes to run the forecast in parallel (uncomment the 3 lines below).
# Alternatively, you can load the ClusterManagers package and add processes
# using one of the schedulers such as SGE or Slurm.
# addprocs(10)
# @everywhere using DSGE
# m <= Setting(:use_parallel_workers, true)

forecast_one(m, :full, :none, output_vars; check_empty_columns = false)
compute_meansbands(m, :full, :none, output_vars; check_empty_columns = false)

# Comment out the line below if you did not run the forecast in parallel.
# rmprocs(procs())


text:
# This file includes examples of symbol transmissions

using DigitalCommunications

# Settings
k = 2
M = 2^k
nsymbols = 10
nbits = nsymbols * k

# Blocks
gen = BitGenerator(nbits)
for modulator in [
        Modulator(PAM(M)),
        Modulator(ASK(M)),
        Modulator(PSK(M)),
        Modulator(QAM(M))
    ]
    tx = gen.bits |> modulator
    @show(tx)
end


text:
# Simulation

type Simulation
    YAML_filename::String
    grid::Grids.Grid
    flow_data :: Array

    configure::Function
    setup::Function

    function Simulation(YAML_filename="IBFS_config_def.yml")
        this = new()
        this.YAML_filename = YAML_filename

        # ------------------- Configuration ------------------- #
        this.configure = function()
            println("------ Configure ...")
            # Configuration file
            println("Configuration file: ", this.YAML_filename)
            config_options = YAML.load(open(this.YAML_filename))
            config_options["solver"]["cg_eps"] = float(config_options["solver"]["cg_eps"])
            println(config_options)

            # Grid instance and configuration
            if config_options["grid"]["implementation"] == "StaggeredGrid2D"
                this.grid = Grids.StaggeredGrid2D()
            else
                println("Not implemented option.")
            end
            this.grid.configure(config_options["grid"])

            # Flow data instances of length nt
            this.flow_data = Array(DataTypes.FlowSnapshot,config_options["solver"]["nt"])
        end #fun config

        # ------------------- Initialization ------------------- #
        this.setup = function()
            println("------ Setup grid and initialize variables ...")
        end #fun init

        return this
    end

end #type


text:
import Makie

#Construct isaPlot3d from type AMFMmodel called with a time index using Makie backend
function isaPlot3d_Makie(z::AMFMmodel,t::Vector{Float64})
    a_max = 1 #need to finish
    myScene = Makie.Scene(resolution = (1080,1080),
    camera = Makie.cam3d!)
    for i in 1:length(z.comps)
        myScene = Makie.lines!(-t,
         -z.comps[i].ω.(t),
          real.(z.comps[i](t)),
           color = ISA.cmap[max.(min.(round.(Int, abs.(z.comps[i].a.(t)) .* 256 / a_max), 256), 1)])
    end
              axis = myScene[Makie.Axis] # get the axis object from the scene
              axis.names.axisnames = ("t", "ω(t)","x(t)")
              myScene
end

#lims = FRect3D(Vec3f0(-1, -1, -1), Vec3f0(1, 1, 2))

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
#=

This file contains all types that are used in the Scattering module:

- `AbstractAerosolTypes` specify aerosol properties
- `AbstractFourierDecompositionTypes` specify the decomposition method to use (NAI2 vs PCW)
- `AbstractPolarizationTypes` specify the polarization type (I/IQU/IQUV)
- `AbstractTruncationTypes` specify the type of truncation for legendre terms
- `GreekCoefs` holds all greek coefficients
- `ScatteringMatrix` holds all computed phase function elements
- `AerosolOptics` holds all computed aerosol optical properties

=#

"""
    type AbstractAerosolType
Abstract aerosol type
"""
abstract type AbstractAerosolType end

"Aerosol type with its properties (size distribution and refractive index)"
mutable struct Aerosol{}
    "Univariate size distribution"
    size_distribution::ContinuousUnivariateDistribution
    "Real part of refractive index"
    nᵣ
    "Imag part of refractive index"
    nᵢ
end

# TODO: struct MultivariateAerosol{FT,FT2} <: AbstractAerosolType

#=

Types of Fourier Decomposition (NAI2 or PCW)

=#

"""
    type AbstractFourierDecompositionType
Abstract aerosol Fourier Decomposition computation type (NAI2 and PCW)
"""
abstract type AbstractFourierDecompositionType end

"""
    type NAI2

Perform Siewart's numerical integration method, NAI-2, to compute aerosol phase function
decomposition. See: http://adsabs.harvard.edu/full/1982A%26A...109..195S
"""
struct NAI2  <: AbstractFourierDecompositionType end

"""
    type PCW

Perform Domke's Precomputed Wigner Symbols method, PCW, to compute aerosol phase function
decomposition. See: http://adsabs.harvard.edu/full/1984A%26A...131..237D
"""
struct PCW <: AbstractFourierDecompositionType end

#=

Types of Polarization (which Stokes vector to use)

=#

"""
    type AbstractPolarizationType

Abstract Polarization type
"""
abstract type AbstractPolarizationType  end

"""
    struct Stokes_IQUV{FT<:AbstractFloat}

A struct which defines full Stokes Vector ([I,Q,U,V]) RT code

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Stokes_IQUV{FT<:AbstractFloat} <: AbstractPolarizationType
    "Number of Stokes components (int)"
    n::Int = 4
    "Vector of length `n` for ... (see eq in Sanghavi )"
    D::Array{FT}  = [1.0, 1.0, -1.0, -1.0]
    "Incoming Stokes vector for scalar only"
    I₀::Array{FT} = [1.0, 0.0, 0.0, 0.0] #assuming completely unpolarized incident stellar radiation
end

"""
    struct Stokes_IQU{FT<:AbstractFloat}

A struct which defines Stokes Vector ([I,Q,U]) RT code

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Stokes_IQU{FT<:AbstractFloat} <: AbstractPolarizationType
    "Number of Stokes components (int)"
    n::Int = 3
    "Vector of length `n` for ... (see eq in Sanghavi )"
    D::Array{FT}  = [1.0, 1.0, -1.0]
    "Incoming Stokes vector for scalar only"
    I₀::Array{FT} = [1.0, 0.0, 0.0] #assuming linearly unpolarized incident stellar radiation
end

"""
    struct Stokes_I{FT<:AbstractFloat}

A struct which define scalar I only RT code

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Stokes_I{FT<:AbstractFloat} <: AbstractPolarizationType
    "Number of Stokes components (int)"
    n::Int = 1
    "Vector of length `n` for ... (see eq in Sanghavi )"
    D::Array{FT} = [1.0]
    "Incoming Stokes vector for scalar only"
    I₀::Array{FT} = [1.0]
end

#=

Types of Truncation (for Legendre terms)

=#

"""
    type AbstractTruncationType

Abstract greek coefficient truncation type
"""
abstract type AbstractTruncationType end

"""
    type δBGE{FT} <: AbstractTruncationType

# Fields
$(DocStringExtensions.FIELDS)
"""
struct δBGE{FT} <: AbstractTruncationType
    "Trunction length for legendre terms"
    l_max::Int
    "Exclusion angle for forward peak (in fitting procedure) `[degrees]`"
    Δ_angle::FT
end

#=

Model that specifies the Mie computation details

=#

"""
    type MieModel

Model to hold all Mie computation details for NAI2 and PCW

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct MieModel{FDT<:AbstractFourierDecompositionType, FT}

    computation_type::FDT
    aerosol::Aerosol
    λ::Real
    polarization_type::AbstractPolarizationType
    truncation_type::AbstractTruncationType

    "Maximum radius `[μm]`"
    r_max::FT
    "Number of quadrature points for integration over size distribution"
    nquad_radius::Int

    wigner_A = zeros(1, 1, 1)
    wigner_B = zeros(1, 1, 1)

end

#=

Types that are needed for the output of the Fourier decomposition

=#

"""
    struct GreekCoefs{FT}

A struct which holds all Greek coefficient lists (over l) in one object.
See eq 16 in Sanghavi 2014 for details.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct GreekCoefs{FT<:Union{AbstractFloat, ForwardDiff.Dual}}
    "Greek matrix coefficient α, is in B[2,2]"
    α::Array{FT,1}
    "Greek matrix coefficient β, is in B[1,1] (only important one for scalar!)"
    β::Array{FT,1}
    "Greek matrix coefficient γ, is in B[2,1],B[1,2]"
    γ::Array{FT,1}
    "Greek matrix coefficient δ, is in B[4,4]"
    δ::Array{FT,1}
    "Greek matrix coefficient ϵ, is in B[3,4] and - in B[4,3]"
    ϵ::Array{FT,1}
    "Greek matrix coefficient ζ, is in B[3,3]"
    ζ::Array{FT,1}
end

""" Extend Base.isapprox (≈) to compare two GreekCoefs """
function Base.:isapprox(greek_coefs_a::GreekCoefs, greek_coefs_b::GreekCoefs)
    field_names = fieldnames(GreekCoefs)
    return all([getproperty(greek_coefs_a, field) ≈ getproperty(greek_coefs_b, field) for field in field_names])
end

"""
    struct ScatteringMatrix

A struct which holds all computed phase function elements.
f₁₁ represents the phase function p for the Intensity (first Stokes Vector element) and is normalized as follows:
1/4π ∫₀²⁽ᵖⁱ⁾ dϕ ∫₋₁¹ p(μ) dμ  = 1

# Fields
$(DocStringExtensions.FIELDS)
"""
struct ScatteringMatrix{FT}
    f₁₁::FT
    f₁₂::FT
    f₂₂::FT
    f₃₃::FT
    f₃₄::FT
    f₄₄::FT
end

"""
    struct AerosolOptics

A struct which holds all computed aerosol optics

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct AerosolOptics{FT<:Union{AbstractFloat, ForwardDiff.Dual}}
    "Greek matrix"
    greek_coefs::GreekCoefs
    "Single Scattering Albedo"
    ω̃::FT
    "Extinction cross-section"
    k::FT
    "Truncation factor"
    fᵗ::FT
    "Derivatives"
    derivs = zeros(1)
end

""" Extend Base.isapprox (≈) to compare two AerosolOptics """
function Base.:isapprox(aerosol_optics_a::AerosolOptics, aerosol_optics_b::AerosolOptics)
    field_names = fieldnames(AerosolOptics)
    return all([getproperty(aerosol_optics_a, field) ≈ getproperty(aerosol_optics_b, field) for field in field_names])
end


text:
@testset "Discrete Explicit Tiger" begin
    S = [:left, :right]
    A = [:left, :right, :listen]
    O = [:left, :right]
    γ = 0.95

    function T(s, a, sp)
        if a == :listen
            return s == sp
        else # a door is opened
            return 0.5 #reset
        end
    end

    function Z(a, sp, o)
        if a == :listen
            if o == sp
                return 0.85
            else
                return 0.15
            end
        else
            return 0.5
        end
    end

    function R(s, a)
        if a == :listen
            return -1.0
        elseif s == a # the tiger was found
            return -100.0
        else # the tiger was escaped
            return 10.0
        end
    end

    m = DiscreteExplicitPOMDP(S,A,O,T,Z,R,γ)

    solver = FunctionSolver(x->:listen)
    policy = solve(solver, m)
    updater = DiscreteUpdater(m)

    rsum = 0.0
    for (s,b,a,o,r) in stepthrough(m, policy, updater, "s,b,a,o,r", max_steps=10)
        println("s: $s, b: $([pdf(b,s) for s in S]), a: $a, o: $o")
        rsum += r
    end
    println("Undiscounted reward was $rsum.")
    @test rsum == -10.0

    dm = DiscreteExplicitPOMDP(S,A,O,T,Z,R,γ,Deterministic(:left))
    @test rand(Random.GLOBAL_RNG, initialstate(dm)) == :left
    tm = DiscreteExplicitPOMDP(S,A,O,T,Z,R,γ,terminals=Set(S))
    @test isterminal(tm, rand(Random.GLOBAL_RNG, initialstate(tm)))

    for s in states(m)
        @test convert_s(statetype(m), convert_s(Vector{Float64}, s, m), m) == s
    end
    for a in actions(m)
        @test convert_a(actiontype(m), convert_a(Vector{Float64}, a, m), m) == a
    end
    for o in observations(m)
        @test convert_o(obstype(m), convert_o(Vector{Float64}, o, m), m) == o
    end
end

@testset "Discrete Explicit MDP" begin
    S = 1:5
    A = [-1, 1]
    γ = 0.95
    p₀ = Deterministic(1)

    function T(s, a, sp)
        if sp == clamp(s+a,1,5)
            return 0.8
        elseif sp == clamp(s-a,1,5)
            return 0.2
        else
            return 0.0
        end
    end

    function R(s, a)
        if s == 5
            return 1.0
        else
            return -1.0
        end
    end

    m = DiscreteExplicitMDP(S,A,T,R,γ)
    m = DiscreteExplicitMDP(S,A,T,R,γ,p₀)
    m = DiscreteExplicitMDP(S,A,T,R,γ,p₀,terminals=Set(5))
    @test isterminal(m, 5)

    solver = FunctionSolver(x->1)
    policy = solve(solver, m)

    rsum = 0.0
    for (s,a,r) in stepthrough(m, policy, "s,a,r", max_steps=10)
        println("s: $s, a: $a")
        rsum += r
    end
    println("Undiscounted reward was $rsum.")

    for s in states(m)
        @test convert_s(statetype(m), convert_s(Vector{Float64}, s, m), m) == s
    end
    for a in actions(m)
        @test convert_a(actiontype(m), convert_a(Vector{Float64}, a, m), m) == a
    end
end

@testset "convert Number POMDP" begin
    s = [1,47]
    T(s, a, sp) = s == sp
    Z(a, sp, o) = o == sp
    R(s, a) = 1.0
    m = DiscreteExplicitPOMDP(s,s,s,T,Z,R,0.95)
    for s in states(m)
        @test convert_s(statetype(m), convert_s(Vector{Float64}, s, m), m) == s
    end
    for a in actions(m)
        @test convert_a(actiontype(m), convert_a(Vector{Float64}, a, m), m) == a
    end
    for o in observations(m)
        @test convert_o(obstype(m), convert_o(Vector{Float64}, o, m), m) == o
    end
end

@testset "convert Vector POMDP" begin
    s = [[1,2],[3,4]]
    T(s, a, sp) = s == sp
    Z(a, sp, o) = o == sp
    R(s, a) = 1.0
    m = DiscreteExplicitPOMDP(s,s,s,T,Z,R,0.95)
    for s in states(m)
        @test convert_s(statetype(m), convert_s(Vector{Float64}, s, m), m) == s
    end
    for a in actions(m)
        @test convert_a(actiontype(m), convert_a(Vector{Float64}, a, m), m) == a
    end
    for o in observations(m)
        @test convert_o(obstype(m), convert_o(Vector{Float64}, o, m), m) == o
    end
end


text:
using SinGAN

using FileIO
using Flux

function main()
    # hp = HyperParams()
    hp = load_hyperparams("hyperparams.json")
    # save_hyperparams("hyperparams.json", hp)
    dscrp, genp = setup_models(hp)
    @info dscrp
    @info genp

    img_name = "floral_shoppe.jpg"
    orig_rgb_img = load(img_name)
    orig_img = SinGAN.rgb_to_array(orig_rgb_img) |> gpu
    img_shapes = image_shapes(hp)
    real_img_p = SinGAN.build_image_pyramid(orig_img, img_shapes, img_shapes)
    # load_model_params!(dscrp, genp, Base.length(img_shapes))

    amplifiers, z_rec = train!(dscrp, genp, real_img_p, hp)
    generate_animation(genp, 0.1f0, 0.9f0, first(amplifiers), z_rec)
end

main()

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
module FileTools
filesep = Base.Filesystem.path_separator;
export filesep
# Write your package code here.
include("datalist.jl");

include("openit.jl");
export openit

include("filelist.jl")
export filelist

include("filelistall.jl")
export filelistall

include("gettag.jl")
export gettag

include("folderlist.jl")
export folderlist, folderlistall

using CSV
include("skipwrite.jl")
export skipwrite

using Printf
include("pathnorepeat.jl")
export pathnorepeat

include("trymkpath.jl")
export trymkpath

include("movefiles.jl")
export movefiles
end


text:
# fg! -- objective function and gradient (avoiding recomputations)
# Hv! -- application of the Hessian
# smooth_fg! -- objective associated with the smooth part of the objective

# ----------------------- #
#  -- Ridge Regression -- #
# ----------------------- #
# ->  f(θ)  = |Xθ - y|₂²/2 + λ|θ|₂²/2
# -> ∇f(θ)  = X'(Xθ - y) + λθ
# -> ∇²f(θ) = X'X + λI
# NOTE:
# * Hv! used in iterative solution
# ---------------------------------------------------------

function Hv!(glr::GLR{L2Loss,<:L2R}, X, y, scratch)
    n, p = size(X)
    λ    = get_penalty_scale(glr, n)
    if glr.fit_intercept
        # H = [X 1]'[X 1] + λ I
        # rows a 1:p = [X'X + λI | X'1]
        # row  e end = [1'X      | n+λι] where ι is 1 if glr.penalize_intercept
        ι = float(glr.penalize_intercept)
        (Hv, v) -> begin
            # view on the first p rows
            a     = 1:p
            Hvₐ   = view(Hv, a)
            vₐ    = view(v,  a)
            Xt1   = view(scratch.p, a)
            Xt1 .*= 0
            @inbounds for i in a, j in 1:n
                Xt1[i] += X[j, i]           # -- X'1
            end
            vₑ  = v[end]
            # update for the first p rows   -- (X'X + λI)v[1:p] + (X'1)v[end]
            Xvₐ = scratch.n
            mul!(Xvₐ, X, vₐ)
            mul!(Hvₐ, X', Xvₐ)
            Hvₐ .+= λ .* vₐ .+ Xt1 .* vₑ
            # update for the last row       -- (X'1)'v + n v[end]
            Hv[end] = dot(Xt1, vₐ) + (n + λ_if_penalize_intercept(glr, λ)) * vₑ
        end
    else
        (Hv, v) -> begin
            Xv = scratch.n
            mul!(Xv, X, v)       # -- Xv
            mul!(Hv, X', Xv)     # -- X'Xv
            Hv .+= λ .* v        # -- X'Xv + λv
        end
    end
end

# ----------------------------- #
#  -- Lasso/Elnet Regression -- #
# ----------------------------- #
# ->  J(θ)  = f(θ) + r(θ)
# ->  f(θ)  = |Xθ - y|₂²/2 + λ|θ|₂²  // smooth
# ->  r(θ)  = γ|θ|₁                  // non-smooth with prox
# -> ∇f(θ)  = X'(Xθ - y) + λθ
# -> ∇²f(θ) = X'X + λI
# -> prox_r = soft-thresh
# ---------------------------------------------------------

function smooth_fg!(glr::GLR{L2Loss,<:ENR}, X, y, scratch)
    λ = get_penalty_scale_l2(glr, length(y))
    (g, θ) -> begin
        # cache contains the residuals (Xθ-y)
        r = scratch.n
        get_residuals!(r, X, θ, y) # -- r = Xθ-y
        apply_Xt!(g, X, r)
        g .+= λ .* θ
        glr.fit_intercept && (glr.penalize_intercept || (g[end] -= λ * θ[end]))
        return glr.loss(r) + get_l2(glr.penalty)(view_θ(glr, θ))
    end
end


text:
@defcomp co2forcing begin
    c_CO2concentration=Parameter(index=[time],unit="ppbv")
    f0_CO2baseforcing=Parameter(unit="W/m2", default=1.68)
    fslope_CO2forcingslope=Parameter(unit="W/m2", default=5.5)
    c0_baseCO2conc=Parameter(unit="ppbv", default=400859.5833)
    f_CO2forcing=Variable(index=[time],unit="W/m2")

    function run_timestep(p, v, d, t)

        #eq.13 in Hope 2006
        v.f_CO2forcing[t]=p.f0_CO2baseforcing+p.fslope_CO2forcingslope*log(p.c_CO2concentration[t]/p.c0_baseCO2conc)
    end
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using Jecco
using HDF5
#using Plots
#gr()
include("./Fourier_modes.jl")

dirname = "/Users/apple/Documents/Jecco.jl/data/old_potential_e_10/"
#Specify the initial iteration, the final and the desired step.
initial_it = 0
step_it = 500
final_it = 100000

#Phi0 and inverse of phiM squared
phi0 = 1.0
phiM2 = -(2.3^2)

To_hdf5 = "no"
To_mathematica = "yes"
Fourier_modes = "no"
order_x, order_y = (4,0)
integration_method = "trapezoidal"
#Read_tensor = "yes"

files = OpenPMDTimeSeries(dirname,prefix="boundary")
files_xi = OpenPMDTimeSeries(dirname,prefix="gauge")

u, x, y = get_field(files,it=initial_it,field="a4")[2][:]
nt = Int(floor((final_it-initial_it)/step_it))+1
# Add the last point
append!(x, x[end]+(x[2]-x[1]))
append!(y, y[end]+(y[2]-y[1]))
nx, ny = length(x), length(y)

e, Jx, Jy = zeros(nt,nx,ny), zeros(nt,nx,ny), zeros(nt,nx,ny)
px, py, pxy, pz = zeros(nt,nx,ny), zeros(nt,nx,ny), zeros(nt,nx,ny), zeros(nt,nx,ny)
O_phi = zeros(nt,nx,ny)
t = zeros(nt)
#Matrix of coefficients of fourier modes of the energy density.
if Fourier_modes == "yes"
    a, b = zeros(nt,order_x+1,order_y+1), zeros(nt,order_x+1,order_y+1)
end

for it in initial_it:step_it:final_it
    #We first read all the files
    a4 = get_field(files,it=it,field="a4")[1][1,:,:]
    b14 = get_field(files,it=it,field="b14")[1][1,:,:]
    b24 = get_field(files,it=it,field="b24")[1][1,:,:]
    fx2 = get_field(files,it=it,field="fx2")[1][1,:,:]
    fy2 = get_field(files,it=it,field="fy2")[1][1,:,:]
    g4 = get_field(files,it=it,field="g4")[1][1,:,:]
    phi2 = get_field(files,it=it,field="phi2")[1][1,:,:]
    #xi = get_field(files_xi,it=it,field="xi")[1][1,:,:]

    i = Int(floor((it-initial_it)/step_it))+1
    #phi2 = phi0^3 .*phi3 .- phi0.*xi.^2

    e[i,1:end-1,1:end-1] = -3/4 .*a4./phi0^4 - phi2./phi0.^3 .- (9-7*phiM2)/(36*phiM2)
    Jx[i,1:end-1,1:end-1] = fx2./phi0^4
    Jy[i,1:end-1,1:end-1] = fy2./phi0^4
    px[i,1:end-1,1:end-1] = -a4./(4*phi0^4) - b14./phi0^4 - b24./phi0^4 + phi2./(3*phi0^3) .+ (-5/108 + 1/(4*phiM2))
    py[i,1:end-1,1:end-1] = -a4./(4*phi0^4) + b14./phi0^4 - b24./phi0^4 + phi2./(3*phi0^3) .+ (-5/108 + 1/(4*phiM2))
    pxy[i,1:end-1,1:end-1] = -g4./phi0^4
    pz[i,1:end-1,1:end-1] = -a4./(4*phi0^4) + 2 .*b24./phi0^4 + phi2./(3*phi0^3) .+ (-5/108 + 1/(4*phiM2))
    O_phi[i,1:end-1,1:end-1] = -2 .*phi2./phi0^3 .+ (1/3 - 1/phiM2)

    t[i] = files.current_t

    #Add the last points
    e[i,end,:], e[i,:,end], e[i,end,end] = e[i,1,:], e[i,:,1], e[i,1,1]
    Jx[i,end,:], Jx[i,:,end], Jx[i,end,end] = Jx[i,1,:], Jx[i,:,1], Jx[i,1,1]
    Jy[i,end,:], Jy[i,:,end], Jy[i,end,end] = Jy[i,1,:], Jy[i,:,1], Jy[i,1,1]
    px[i,end,:], px[i,:,end], px[i,end,end] = px[i,1,:], px[i,:,1], px[i,1,1]
    py[i,end,:], py[i,:,end], px[i,end,end] = py[i,1,:], py[i,:,1], py[i,1,1]
    pxy[i,end,:], pxy[i,:,end], pxy[i,end,end] = pxy[i,1,:], pxy[i,:,1], pxy[i,1,1]
    pz[i,end,:], pz[i,:,end], pz[i,end,end] = pz[i,1,:], pz[i,:,1], pz[i,1,1]
    O_phi[i,end,:], O_phi[i,:,end], O_phi[i,end,end] = O_phi[i,1,:], O_phi[i,:,1], O_phi[i,1,1]

    if Fourier_modes == "yes"
        a[i,:,:], b[i,:,:] = Fourier_modes_2D(e[i,:,:], x, y, order_x, order_y, integration_method)
    end

end

#Writing the data into a h5df file.
if To_hdf5 == "yes"

    output = string(dirname,"Stress-Tensor.h5df")
    fid = h5open(output,"w")

    group_coord = g_create(fid, "Coordinates")
    group_st = g_create(fid, "Stress-Tensor")
    group_svev = g_create(fid, "Scalar_VEV")
    group_fourier = g_create(fid, "Fourier modes")

    Jecco.write_dataset(group_coord, "t",t)
    Jecco.write_dataset(group_coord, "x",x)
    Jecco.write_dataset(group_coord, "y",y)

    Jecco.write_dataset(group_st, "e",e)
    Jecco.write_dataset(group_st, "Jx",Jx)
    Jecco.write_dataset(group_st, "Jy",Jy)
    Jecco.write_dataset(group_st, "px",px)
    Jecco.write_dataset(group_st, "py",py)
    Jecco.write_dataset(group_st, "pxy",pxy)
    Jecco.write_dataset(group_st, "pz",pz)

    Jecco.write_dataset(group_svev, "O_phi",O_phi)

    Jecco.write_dataset(group_fourier, "a",a)
    Jecco.write_dataset(group_fourier, "b",b)

    close(fid)
end

#This writes into a h5df file that has the mathematica array format directly
if To_mathematica == "yes"
    T_m = zeros(11,nx*ny*nt)
    if Fourier_modes == "yes"
        f_m = zeros(2*(order_x+1)*(order_y+1)+1,nt)
    end
    n = 0
    for i in 1:nt
        for j in 1:nx
            for k in 1:ny
                global n += 1
                T_m[:,n] = [t[i] x[j] y[k] e[i,j,k] Jx[i,j,k] Jy[i,j,k] px[i,j,k] pxy[i,j,k] py[i,j,k] pz[i,j,k] O_phi[i,j,k]]
            end
        end
        if Fourier_modes == "yes"
            f_m[:,i] = hcat(t[i], reshape(transpose(a[i,:,:]),1,:), reshape(transpose(b[i,:,:]),1,:))
        end
    end

    output = string(dirname,"data_mathematica.h5df")
    fid = h5open(output,"w")

    group_st = g_create(fid, "data")

    Jecco.write_dataset(group_st, "VEVs",T_m)
    if Fourier_modes == "yes"
        Jecco.write_dataset(group_st, "Fourier modes",f_m)
    end

    close(fid)

end


text:
using SparseArrays
using LinearAlgebra
using jInv.Mesh
using Multigrid
using Multigrid.DomainDecomposition
using Helmholtz

domain = [0.0, 1.0, 0.0, 1.0];
n = [64,64];
Mr = getRegularMesh(domain,n)

pad = 8;
m = ones(tuple(n.+1...));

maxOmega = getMaximalFrequency(m,Mr);
ABLamp = maxOmega;

w = 0.9*maxOmega;
println("omega*h:");
println(w*Mr.h*sqrt(maximum(m)));
pad = pad*ones(Int64,Mr.dim);

println("************************************************* Hybrid Keczmarz for Acoustic Helmholtz 2D ******************************************************");

Sommerfeld = true;
NeumannAtFirstDim = true;
H,gamma = GetHelmholtzOperator(Mr,m,w,w*ones(size(m))*0.01,NeumannAtFirstDim,pad,ABLamp,Sommerfeld);
shift = 0.5;
SHT = sparse((H .+ GetHelmholtzShiftOP(m, real(w),shift))');

numCores = 2;
innerKatzIter = 3;
inner = 3;

q = getAcousticPointSource(Mr,ComplexF64)[1];
x = zeros(ComplexF64,size(q));
N = size(SHT,2);

HKparam = getHybridKaczmarz(ComplexF64,Int64, SHT,Mr,[4,4], getNodalIndicesOfCell,0.95,numCores,innerKatzIter);
prec = getHybridKaczmarzPrecond(HKparam,SHT);
x[:] .= FGMRES_relaxation(SHT,q[:],x[:],inner,prec,1e-3,true,true,numCores)[1]
println("Relative norm is: ",norm(SHT'*x[:] - q[:])./norm(q[:]))

println("In float precision:")

x[:].=0.0;
x = convert(Array{ComplexF32},x);
q = convert(Array{ComplexF32},q);
SHT = convert(SparseMatrixCSC{ComplexF32,Int64},SHT);

HKparam = getHybridKaczmarz(ComplexF32,Int64, SHT,Mr,[4,4], getNodalIndicesOfCell,0.95,numCores,innerKatzIter);
prec = getHybridKaczmarzPrecond(HKparam,SHT);
x[:] .= FGMRES_relaxation(SHT,q[:],x[:],inner,prec,1e-3,true,true,numCores)[1]
println("Relative norm is: ",norm(SHT'*x[:] - q[:])./norm(q[:]))


text:
#####
# UM
#####
export UMOracle

#returns zstar, ustar
function suppFcnUM(xs, lquants, uquants, cut_sense=:Max)
    const d = length(lquants)

    zstar = 0.0
    for i in 1:d
        mn, mx = minmax(lquants[i] * xs[i], uquants[i] * xs[i])
        zstar += ifelse(cut_sense == :Min, mn, mx)
    end

    ustar = Array(Float64, d)
    posxopt = (cut_sense == :Min) ? lquants : uquants
    negxopt = (cut_sense == :Min) ? uquants : lquants
    for i in 1:d
       ustar[i] = xs[i] > 0 ? posxopt[i] : negxopt[i]
    end

    return zstar, ustar
end

function calc_s(data, eps_, alpha)
	N, d = size(data)
	if (1-eps_/d)^N  > alpha / 2d
		return N + 1
	else
		dBin = Binomial(N, 1-eps_/d)
		return quantile(dBin, 1-alpha/2d)
	end
end

######################
type UMOracle <: AbstractOracle
    lquants::Vector{Float64}
    uquants::Vector{Float64}
    cut_tol::Float64  ##defaults to 1e-6

    # Other options
    debug_printcut::Bool
end

suppFcn(xs::Vector, w::UMOracle, cut_sense) =
        suppFcnUM(xs, w.lquants, w.uquants, cut_sense)

#Preferred Interface
function UMOracle(data, lbounds, ubounds, eps_, alpha; cut_tol=1e-6, debug_printcut=false)
	s = calc_s(data, eps_, alpha)
	N = size(data, 1)
	@assert N- s + 1 < s "UM: N not sufficiently big N: $N \t s: $s"
	if s == N + 1
		return UMOracle(lbounds, ubounds, cut_tol, debug_printcut)
	else
		data_sort = sort_data_cols(data)
		return UMOracle(vec(data_sort[N-s+1, :]), vec(data_sort[s, :]), cut_tol, debug_printcut)
	end
end

function setup(w::UMOracle, rm::Model, prefs)
    # Extract preferences we care about
    w.debug_printcut = get(prefs, :debug_printcut, false)
    w.cut_tol        = get(prefs, :cut_tol, w.cut_tol)
    rd = JuMPeR.getRobust(rm)
    @assert (rd.numUncs == length(w.lquants)) "Num Uncertainties $(rd.numUncs) doesn't match columns in lquants $(w.lquants)"
    @assert (rd.numUncs == length(w.uquants)) "Num Uncertainties $(rd.numUncs) doesn't match columns in uquants $(w.uquants)"
    @assert (length(rd.uncertaintyset) == 0) "Auxiliary constraints on uncertainties not yet supported"
end


text:
using Pkg
Pkg.activate(pwd())

# # Tuples

t = (1, 2.0, "3")
t = 1, 2.0, "3"

#+

typeof(t)

#+

t[1] # the first element
t[end] # the last element
t[1:2] # the first two elements

#+

a, b, c = t
println("The values stored in the tuple are: $a, $b and $c")

# ### Exercise:
# Create a tuple that contains the first four letters of the alphabet (these letters should
# be of type `String`). Then unpack this tuple into four variables `a`, `b`, `c` and `d`.
#
# ---
# ### Solution:

# ---
#
# # Named Tuples

t = (a = 1, b = 2.0, c = "3")

#+

a = 1;
b = 2.0;
c = "3";

t = (; a, b, c)

#+

t[1] # the first element
t[end] # the last element
t[1:2] # error

#+

t.a
t.c
a, b, c = t

println("The values stored in the tuple are: a = $a, b = $b")


text:
using SimulatedNeuralMoments, Statistics, DelimitedFiles
results = readdlm("mcresults.txt")
sofar = 500
println("__________ replication: ", sofar, "_______________")
clabels = ["99%", "95%", "90%"]
prettyprint(reshape(mean(results[1:sofar,8:end],dims=1),7,3),clabels)
dstats(results[1:sofar,1:7])

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
# LLVM IR optimization

function optimize!(ctx::CompilerContext, mod::LLVM.Module, entry::LLVM.Function)
    tm = AMDGPUnative.machine(ctx.agent, triple(mod))

    if ctx.kernel
        entry = promote_kernel!(ctx, mod, entry)
    end

    let pm = ModulePassManager()
        global global_ctx
        global_ctx = ctx

        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        internalize!(pm, [LLVM.name(entry)])

        # lower intrinsics
        add!(pm, FunctionPass("LowerGCFrame", lower_gc_frame!))
        aggressive_dce!(pm) # remove dead uses of ptls
        add!(pm, ModulePass("LowerPTLS", lower_ptls!))

        ccall(:jl_add_optimization_passes, Cvoid,
              (LLVM.API.LLVMPassManagerRef, Cint),
              LLVM.ref(pm), Base.JLOptions().opt_level)

        # NVPTX's target machine info enables runtime unrolling,
        # but Julia's pass sequence only invokes the simple unroller.
        loop_unroll!(pm)
        instruction_combining!(pm)  # clean-up redundancy
        licm!(pm)                   # the inner runtime check might be outer loop invariant

        # the above loop unroll pass might have unrolled regular, non-runtime nested loops.
        # that code still needs to be optimized (arguably, multiple unroll passes should be
        # scheduled by the Julia optimizer). do so here, instead of re-optimizing entirely.
        early_csemem_ssa!(pm) # TODO: gvn instead? see NVPTXTargetMachine.cpp::addEarlyCSEOrGVNPass
        dead_store_elimination!(pm)

        constant_merge!(pm)

        # NOTE: if an optimization is missing, try scheduling an entirely new optimization
        # to see which passes need to be added to the list above
        #     LLVM.clopts("-print-after-all", "-filter-print-funcs=$(LLVM.name(entry))")
        #     ModulePassManager() do pm
        #         add_library_info!(pm, triple(mod))
        #         add_transform_info!(pm, tm)
        #         PassManagerBuilder() do pmb
        #             populate!(pm, pmb)
        #         end
        #         run!(pm, mod)
        #     end

        cfgsimplification!(pm)

        run!(pm, mod)
        dispose(pm)
    end

    # we compile a module containing the entire call graph,
    # so perform some interprocedural optimizations.
    #
    # for some reason, these passes need to be distinct from the regular optimization chain,
    # or certain values (such as the constant arrays used to populare llvm.compiler.user ad
    # part of the LateLowerGCFrame pass) aren't collected properly.
    #
    # these might not always be safe, as Julia's IR metadata isn't designed for IPO.
    let pm = ModulePassManager()
        dead_arg_elimination!(pm)   # parent doesn't use return value --> ret void

        run!(pm, mod)
        dispose(pm)
    end

    return entry
end

## kernel-specific optimizations

# promote a function to a kernel
# FIXME: sig vs tt (code_llvm vs rocfunction)
function promote_kernel!(ctx::CompilerContext, mod::LLVM.Module, entry_f::LLVM.Function)
    kernel = wrap_entry!(ctx, mod, entry_f)

    # property annotations
    annotations = LLVM.Value[kernel]

    ## kernel metadata
    append!(annotations, [MDString("kernel"), ConstantInt(Int32(1), JuliaContext())])

    return kernel
end

function wrapper_type(julia_t::Type, codegen_t::LLVMType)::LLVMType
    if !isbitstype(julia_t)
        # don't pass jl_value_t by value; it's an opaque structure
        return codegen_t
    elseif isa(codegen_t, LLVM.PointerType) && !(julia_t <: Ptr)
        # we didn't specify a pointer, but codegen passes one anyway.
        # make the wrapper accept the underlying value instead.
        return eltype(codegen_t)
    else
        return codegen_t
    end
end

# generate a kernel wrapper to fix & improve argument passing
function wrap_entry!(ctx::CompilerContext, mod::LLVM.Module, entry_f::LLVM.Function)
    entry_ft = eltype(llvmtype(entry_f)::LLVM.PointerType)::LLVM.FunctionType
    @compiler_assert return_type(entry_ft) == LLVM.VoidType(JuliaContext()) ctx

    # filter out ghost types, which don't occur in the LLVM function signatures
    sig = Base.signature_type(ctx.f, ctx.tt)::Type
    julia_types = Type[]
    for dt::Type in sig.parameters
        if !isghosttype(dt)
            push!(julia_types, dt)
        end
    end

    # generate the wrapper function type & definition
    wrapper_types = LLVM.LLVMType[wrapper_type(julia_t, codegen_t)
                                  for (julia_t, codegen_t)
                                  in zip(julia_types, parameters(entry_ft))]
    wrapper_fn = replace(LLVM.name(entry_f), r"^.+?_"=>"roccall_") # change the CC tag
    wrapper_ft = LLVM.FunctionType(LLVM.VoidType(JuliaContext()), wrapper_types)
    wrapper_f = LLVM.Function(mod, wrapper_fn, wrapper_ft)

    # emit IR performing the "conversions"
    let builder = Builder(JuliaContext())
        entry = BasicBlock(wrapper_f, "entry", JuliaContext())
        position!(builder, entry)

        wrapper_args = Vector{LLVM.Value}()

        # perform argument conversions
        codegen_types = parameters(entry_ft)
        wrapper_params = parameters(wrapper_f)
        param_index = 0
        for (julia_t, codegen_t, wrapper_t, wrapper_param) in
            zip(julia_types, codegen_types, wrapper_types, wrapper_params)
            param_index += 1
            if codegen_t != wrapper_t
                # the wrapper argument doesn't match the kernel parameter type.
                # this only happens when codegen wants to pass a pointer.
                @compiler_assert isa(codegen_t, LLVM.PointerType) ctx
                @compiler_assert eltype(codegen_t) == wrapper_t ctx

                # copy the argument value to a stack slot, and reference it.
                ptr = alloca!(builder, wrapper_t)
                if LLVM.addrspace(codegen_t) != 0
                    ptr = addrspacecast!(builder, ptr, codegen_t)
                end
                store!(builder, wrapper_param, ptr)
                push!(wrapper_args, ptr)
            else
                push!(wrapper_args, wrapper_param)
                for attr in collect(parameter_attributes(entry_f, param_index))
                    push!(parameter_attributes(wrapper_f, param_index), attr)
                end
            end
        end

        call!(builder, entry_f, wrapper_args)

        ret!(builder)

        dispose(builder)
    end

    # early-inline the original entry function into the wrapper
    push!(function_attributes(entry_f), EnumAttribute("alwaysinline", 0, JuliaContext()))
    linkage!(entry_f, LLVM.API.LLVMInternalLinkage)

    # set calling convention to AMDGPU_KERNEL
    # NOTE: only valid to set this on the outermost function
    callconv!(wrapper_f, LLVM.API.LLVMCallConv(91))

    fixup_metadata!(entry_f)
    ModulePassManager() do pm
        always_inliner!(pm)
        verifier!(pm)
        run!(pm, mod)
    end

    return wrapper_f
end

# HACK: get rid of invariant.load and const TBAA metadata on loads from pointer args,
#       since storing to a stack slot violates the semantics of those attributes.
# TODO: can we emit a wrapper that doesn't violate Julia's metadata?
function fixup_metadata!(f::LLVM.Function)
    for param in parameters(f)
        if isa(llvmtype(param), LLVM.PointerType)
            # collect all uses of the pointer
            worklist = Vector{LLVM.Instruction}(user.(collect(uses(param))))
            while !isempty(worklist)
                value = popfirst!(worklist)

                # remove the invariant.load attribute
                md = metadata(value)
                if haskey(md, LLVM.MD_invariant_load)
                    delete!(md, LLVM.MD_invariant_load)
                end
                if haskey(md, LLVM.MD_tbaa)
                    delete!(md, LLVM.MD_tbaa)
                end

                # recurse on the output of some instructions
                if isa(value, LLVM.BitCastInst) ||
                   isa(value, LLVM.GetElementPtrInst) ||
                   isa(value, LLVM.AddrSpaceCastInst)
                    append!(worklist, user.(collect(uses(value))))
                end

                # IMPORTANT NOTE: if we ever want to inline functions at the LLVM level,
                # we need to recurse into call instructions here, and strip metadata from
                # called functions (see CUDAnative.jl#238).
            end
        end
    end
end

# lower object allocations to to PTX malloc
#
# this is a PoC implementation that is very simple: allocate, and never free. it also runs
# _before_ Julia's GC lowering passes, so we don't get to use the results of its analyses.
# when we ever implement a more potent GC, we will need those results, but the relevant pass
# is currently very architecture/CPU specific: hard-coded pool sizes, TLS references, etc.
# such IR is hard to clean-up, so we probably will need to have the GC lowering pass emit
# lower-level intrinsics which then can be lowered to architecture-specific code.
function lower_gc_frame!(fun::LLVM.Function)
    ctx = global_ctx::CompilerContext
    mod = LLVM.parent(fun)
    changed = false

    # plain alloc
    if haskey(functions(mod), "julia.gc_alloc_obj")
        alloc_obj = functions(mod)["julia.gc_alloc_obj"]
        alloc_obj_ft = eltype(llvmtype(alloc_obj))
        T_prjlvalue = return_type(alloc_obj_ft)
        T_pjlvalue = convert(LLVMType, Any, true)

        for use in uses(alloc_obj)
            call = user(use)::LLVM.CallInst

            # decode the call
            ops = collect(operands(call))
            sz = ops[2]

            # replace with PTX alloc_obj
            let builder = Builder(JuliaContext())
                position!(builder, call)
                ptr = call!(builder, Runtime.get(:gc_pool_alloc), [sz])
                replace_uses!(call, ptr)
                dispose(builder)
            end

            unsafe_delete!(LLVM.parent(call), call)

            changed = true
        end

        @compiler_assert isempty(uses(alloc_obj)) ctx
    end

    # we don't care about write barriers
    if haskey(functions(mod), "julia.write_barrier")
        barrier = functions(mod)["julia.write_barrier"]

        for use in uses(barrier)
            call = user(use)::LLVM.CallInst
            unsafe_delete!(LLVM.parent(call), call)
            changed = true
        end

        @compiler_assert isempty(uses(barrier)) ctx
    end

    return changed
end

# lower the `julia.ptls_states` intrinsic by removing it, since it is GPU incompatible.
#
# this assumes and checks that the TLS is unused, which should be the case for most GPU code
# after lowering the GC intrinsics to TLS-less code and having run DCE.
#
# TODO: maybe don't have Julia emit actual uses of the TLS, but use intrinsics instead,
#       making it easier to remove or reimplement that functionality here.
function lower_ptls!(mod::LLVM.Module)
    ctx = global_ctx::CompilerContext
    changed = false

    if haskey(functions(mod), "julia.ptls_states")
        ptls_getter = functions(mod)["julia.ptls_states"]

        for use in uses(ptls_getter)
            val = user(use)
            if !isempty(uses(val))
                error("Thread local storage is not implemented")
            end
            unsafe_delete!(LLVM.parent(val), val)
            changed = true
        end

        @compiler_assert isempty(uses(ptls_getter)) ctx
     end

    return changed
end


text:
using Documenter, LifeContingencies

makedocs(;
    modules=[LifeContingencies],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaActuary/LifeContingencies.jl/blob/{commit}{path}#L{line}",
    sitename="LifeContingencies.jl",
    authors="Alec Loudenback",
)

deploydocs(;
    repo="github.com/JuliaActuary/LifeContingencies.jl",
)


text:
#
# Least Squares Estimation
#

"""
    LeastSquaresEstimator(sys::LinearSystem, obs::LinearObserver,
                          estimate::AbstractUncertainState)

Linear batch least squares estimator. The internal estimator model takes on a
specific form depending on the type of the initial estimate. For a
`DiscreteState`,

\$x_k = A x_{k-1} + w_k\$
\$y_k = H x_k + v_k\$

where \$w_k \\sim N(0, Q)\$ and \$v_k \\sim N(0, R)\$. For a `ContinuousState`,

\$\\dot{x}(t_k) = A x(t_k) + w(t_k)\$
\$y(t_k) = H x(t_k) + v(t_k)\$

Construction of a LeastSquaresEstimator requires an initial estimate.
Internally, the initial estimate is an uncertain state type. The constructor
automatically converts absolute initial estimates to uncertain states with zero
covariance.

`LeastSquaresEstimator` does not currently use the value of the initial estimate
or its covariance in solving for a new estimate.
"""
immutable LeastSquaresEstimator{T,S<:AbstractUncertainState{T},
                            M<:AbstractAbsoluteState{T}} <: BatchEstimator{T,S}
    sys::LinearSystem{T}
    obs::LinearObserver{T}
    estimate::S
    measurements::Vector{M}

    function LeastSquaresEstimator(sys::LinearSystem{T}, obs::LinearObserver{T},
        estimate::S) where {T,S<:AbstractUncertainState{T}}
        assert_compatibility(sys, estimate)
        assert_compatibility(obs, estimate)
        M = absolute_type(estimate)
        new{T,S,M}(sys, obs, estimate, Vector{M}([]))
    end
end
function LeastSquaresEstimator(sys::LinearSystem, obs::LinearObserver,
                               estimate::AbstractAbsoluteState)
    LeastSquaresEstimator(sys, obs, make_uncertain(estimate))
end

"""
    NonlinearLeastSquaresEstimator(sys::NonlinearSystem, obs::NonlinearObserver,
        estimate::AbstractState[, tolerance::AbstractFloat, max_iterations::Integer])

Nonlinear least squares estimator type. The internal estimator model takes on a
specific form depending on the type of the initial estimate. For a `DiscreteState`,

\$x_k = F(k, x_{k-1}) + w_k\$
\$y_k = H(k, x_k) + v_k\$

where \$w_k \\sim N(0, Q)\$ and \$v_k \\sim N(0, R)\$. For a `ContinuousState`,

\$\\dot{x}(t_k) = F(t_k, x(t_k)) + w(t_k)\$
\$y(t_k) = H(t_k, x(t_k)) + v(t_k)\$

Construction of a NonlinearLeastSquaresEstimator requires an initial estimate.
Internally, the initial estimate is an uncertain state type. The constructor
automatically converts absolute initial estimates to uncertain states with zero
covariance.

Convergence tolerance and maximum iteration settings are optional. The default
tolerance is 0.01 and the default maximum number of iterations is 15. Convergence
occurs when

\$\\left\\| \\frac{\\Delta x_i}{x_i} \\right\\| \\leq \\epsilon\$

for \$i = 1, \\ldots, n\$ where \$n\$ is the dimension of the state and \$\\epsilon\$
is the convergence tolerance.
"""
immutable NonlinearLeastSquaresEstimator{T,S<:AbstractUncertainState{T},
                            M<:AbstractAbsoluteState{T}} <: BatchEstimator{T,S}
    sys::NonlinearSystem{T}
    obs::NonlinearObserver{T}
    estimate::S
    measurements::Vector{M}

    tolerance::T
    max_iterations::UInt16

    function NonlinearLeastSquaresEstimator(sys, obs, estimate::S, tol=1e-2,
        max_iterations=15) where {T,S<:AbstractUncertainState{T}}
        M = absolute_type(estimate)
        new{T,S,M}(NonlinearSystem(sys), NonlinearObserver(obs), estimate,
                   Vector{M}([]), tol, max_iterations)
    end
end
function NonlinearLeastSquaresEstimator(sys, obs, estimate::AbstractAbsoluteState,
                                        tol=1e-2, max_iterations=15)
    NonlinearLeastSquaresEstimator(sys, obs, make_uncertain(estimate), tol)
end

"""
    add!(estimator::BatchEstimator, measurement::AbstractAbsoluteState)
    add!(estimator::BatchEstimator, measurements::Vector{AbstractAbsoluteState})

Add one or more measurements to a `BatchEstimator` for future processing. The
`BatchEstimator` will store these measurements internally.
"""
function add!(::BatchEstimator) end
function add!(estimator::BatchEstimator, measurement::AbstractAbsoluteState)
    assert_compatibility(measurement, estimator.obs)
    push!(estimator.measurements, measurement)
    return nothing
end
function add!(estimator::BatchEstimator, measurements::Vector)
    for idx = 1:length(measurements)
        add!(estimator, measurements[idx])
    end
    return nothing
end

"""
    compute_residuals!(residuals::Vector, lse::LeastSquaresEstimator[,
                       initial_state::AbstractUncertainState])

Compute residuals for all internal measurements, returning a Vector of
residual states.
"""
function compute_residuals!{T,S}(residuals::Vector{S},
                                 estimator::BatchEstimator{T,S},
                                 initial_state::S)
    for idx = 1:length(estimator.measurements)
        pred_measurement = predict(
            predict(initial_state, estimator.sys, estimator.measurements[idx].t),
            estimator.obs)
        pred_measurement .-= estimator.measurements[idx].x
        push!(residuals, pred_measurement)
    end
    return nothing
end
compute_residuals!(residuals::Vector, estimator::BatchEstimator) =
    compute_residuals!(residuals, estimator, estimator.estimate)

"""
    compute_states(state_history::Vector, lse::LeastSquaresEstimator)

Compute states for all internal measurements, returning a Vector.
"""
function compute_states!{T,S}(state_history::Vector{S},
                              estimator::BatchEstimator{T,S},
                              initial_state::S)
    for idx = 1:length(estimator.measurements)
        pred_state = predict(
            initial_state, estimator.sys, estimator.measurements[idx].t)
        push!(state_history, pred_state)
    end
    return nothing
end
compute_states!(state_history::Vector, estimator::BatchEstimator) =
    compute_states(state_history, estimator, estimator.estimate)

"""
    solve(lse:LeastSquaresEstimator[, archive::EstimatorHistory])

Solve for an updated state estimate using the internal model and measurements
in a `LeastSquaresEstimator`. Returns an uncertain state containing the updated
state estimate and its covariance. Use `solve!` to update the internal estimate
of the `LeastSquaresEstimator` in-place.

Optionally provide an `EstimatorHistory` archive variable to store intermediate
solution data.
"""
function solve(::BatchEstimator) end
function solve{T}(lse::LeastSquaresEstimator{T})
    m,n = size(lse.obs.H)
    A = zeros(T, n, m*length(lse.measurements))
    b = zeros(T, m*length(lse.measurements))
    C = spzeros(T, m*length(lse.measurements), m*length(lse.measurements))

    for idx = 1:length(lse.measurements)
        start_idx = (idx-1)*m + 1
        end_idx = start_idx + m - 1

        A[:, start_idx:end_idx] .= state_transition_matrix(lse.sys,
            lse.estimate, lse.measurements[idx].t)' * lse.obs.H'
        b[start_idx:end_idx] .= lse.measurements[idx].x
        C[start_idx:end_idx, start_idx:end_idx] .= lse.obs.R
    end

    estimate = deepcopy(lse.estimate)
    temp = inv(A*A')*A
    estimate.x = temp * b
    estimate.P = temp * C * temp'

    return estimate
end
function solve{T}(nlse::NonlinearLeastSquaresEstimator{T})
    m,n = size(nlse.measurements[1].x, 1), length(nlse.estimate.x)
    A = zeros(T, n, m*length(nlse.measurements))
    b = zeros(T, m*length(nlse.measurements))
    C = spzeros(T, m*length(nlse.measurements), m*length(nlse.measurements))
    temp = zeros(T, n, m)

    estimate = deepcopy(nlse.estimate)
    convergence = nlse.tolerance+1
    convergence = fill(nlse.tolerance+1, size(nlse.estimate.x))
    iteration_count = nlse.max_iterations

    while any(convergence .> nlse.tolerance) && iteration_count > 0
        for idx = 1:length(nlse.measurements)
            start_idx = (idx-1)*m + 1
            end_idx = start_idx + m - 1

            t = nlse.measurements[idx].t
            predicted_measurement =
                predict(predict(estimate, nlse.sys, t), nlse.obs)

            A[:, start_idx:end_idx] .= state_transition_matrix(nlse.sys,
                estimate, t)' * nlse.obs.dH_dx(t, estimate.x)'
            b[start_idx:end_idx] .=
                nlse.measurements[idx].x - predicted_measurement.x
            C[start_idx:end_idx, start_idx:end_idx] .= nlse.obs.R
        end

        temp = inv(A*A')*A
        delta = temp * b
        estimate.x += delta
        convergence = abs.(delta ./ estimate.x)
        iteration_count -= 1
    end

    estimate.P = temp * C * temp'
    return estimate
end
function solve(estimator::BatchEstimator, archive::EstimatorHistory)
    estimate = solve(estimator)
    compute_residuals!(archive.residuals, estimator, estimate)
    compute_states!(archive.states, estimator, estimate)
    return estimate
end

"""
    solve!(lse:LeastSquaresEstimator[, archive::EstimatorHistory])
"""
function solve!(estimator::BatchEstimator)
    estimator.estimate .= solve(estimator)
    return nothing
end
function solve!(estimator::BatchEstimator, archive::EstimatorHistory)
    estimator.estimate .= solve(estimator, archive)
    return nothing
end


text:
score = 0
n = 5
for i in 1:n
    a = rand(1:10)
    b = rand(1:10)
    print(a, "*", b, " = ")
    c = parse(Int, readline(stdin))
    if a * b == c
        printstyled("Yay! Correct\n", color=:green)
        global score += 1
    else
        printstyled("Dumbass!\n", color=:red)
    end
end
printstyled("Score is $score out of $n\n", color=:yellow)
# this is my change

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
module RK45
import Dates
import Logging
import Printf: @sprintf

#Get Butcher tableau etc from separate file (for convenience of changing if wanted)
include("dopri.jl")

function solve(f!, y0, t, dt, tmax;
               rtol=1e-6, atol=1e-10, safety=0.9, max_dt=Inf, min_dt=0, locextrap=true,
               norm=weaknorm,
               kwargs...)
    stepper = Stepper(f!, y0, t, dt,
                      rtol=rtol, atol=atol, safety=safety, max_dt=max_dt, min_dt=min_dt, locextrap=locextrap, norm=norm)
    return solve(stepper, tmax; kwargs...)
end

function solve_precon(f!, linop, y0, t, dt, tmax;
                    rtol=1e-6, atol=1e-10, safety=0.9, max_dt=Inf, min_dt=0, locextrap=true, norm=weaknorm,
                    kwargs...)
    stepper = PreconStepper(f!, linop, y0, t, dt,
                      rtol=rtol, atol=atol, safety=safety, max_dt=max_dt, min_dt=min_dt, locextrap=locextrap, norm=norm)
    return solve(stepper, tmax; kwargs...)
end

function solve(s, tmax; stepfun=donothing!, output=false, outputN=201,
                        status_period=1, repeat_limit=10)
    if output
        yout = Array{eltype(s.y)}(undef, (size(s.y)..., outputN))
        tout = range(s.t, stop=tmax, length=outputN)
        saved = 1
        yout[fill(:, ndims(s.y))..., 1] = s.y
    end

    steps = 0
    repeated = 0
    repeated_tot = 0

    Logging.@info "Starting propagation"
    start = Dates.now()
    tic = Dates.now()
    while s.tn <= tmax
        ok = step!(s)
        steps += 1
        if Dates.value(Dates.now()-tic) > 1000*status_period
            speed = s.tn/(Dates.value(Dates.now()-start)/1000)
            eta_in_s = (tmax-s.tn)/(speed)
            if eta_in_s > 356400
                Logging.@info @sprintf("Progress: %.2f %%, ETA: XX:XX:XX, stepsize %.2e, err %.2f, repeated %d",
                s.tn/tmax*100, s.dt, s.err, repeated_tot)
            else
                eta_in_ms = Dates.Millisecond(ceil(eta_in_s*1000))
                etad = Dates.DateTime(Dates.UTInstant(eta_in_ms))
                Logging.@info @sprintf("Progress: %.2f %%, ETA: %s, stepsize %.2e, err %.2f, repeated %d",
                    s.tn/tmax*100, Dates.format(etad, "HH:MM:SS"), s.dt, s.err, repeated_tot)
            end
            tic = Dates.now()
        end
        if ok
            if output
                while (saved<outputN) && tout[saved+1] < s.tn
                    ti = tout[saved+1]
                    yout[fill(:, ndims(s.y))..., saved+1] .= interpolate(s, ti)
                    saved += 1
                end
            end
            stepfun(s.yn, s.tn, s.dtn, t -> interpolate(s, t))
            repeated = 0
        else
            repeated += 1
            repeated_tot += 1
            if repeated > repeat_limit
                error("Reached limit for step repetition ($repeat_limit)")
            end
        end
    end
    Logging.@info @sprintf("Propagation finished in %.3f seconds, %d steps",
                           Dates.value(Dates.now()-start)/1000, steps)

    if output
        return collect(tout), yout, steps
    else
        return nothing
    end
end

mutable struct Stepper{T<:AbstractArray, F, nT}
    f!::F  # RHS function
    y::T  # Solution at current t
    yn::T  # Solution at t+dt
    yi::T  # Interpolant array (see interpolate())
    yerr::T  # solution error estimate (from embedded RK)
    ks::NTuple{7, T}  # k values (intermediate solutions for Runge-Kutta method)
    t::Float64  # current time (propagation variable)
    tn::Float64  # next time
    dt::Float64  # time step
    dtn::Float64  # time step for next step
    rtol::Float64  # relative tolerance on error
    atol::Float64  # absolute tolerance on error
    safety::Float64  # safety factor for stepsize control
    max_dt::Float64  # maximum value for dt (default Inf)
    min_dt::Float64  # minimum value for dt (default 0)
    locextrap::Bool  # true if using local extrapolation
    ok::Bool  # true if current step was successful
    err::Float64  # error metric to be compared to tol
    errlast::Float64  # error of the most recent successful step
    norm::nT # function to calculate error metric, defaults to RK45.weaknorm
end

function Stepper(f!, y0, t, dt;
                 rtol=1e-6, atol=1e-10, safety=0.9, max_dt=Inf, min_dt=0,
                 locextrap=true, norm=weaknorm)
    k1 = similar(y0)
    f!(k1, y0, t)
    ks = (k1, similar(k1), similar(k1), similar(k1), similar(k1), similar(k1), similar(k1))
    yerr = similar(y0)
    return Stepper(f!, copy(y0), copy(y0), similar(y0), yerr, ks,
        float(t), float(t), float(dt), float(dt),
        float(rtol), float(atol), float(safety), float(max_dt), float(min_dt),
        locextrap, false, 0.0, 0.0, norm)
end

mutable struct PreconStepper{T<:AbstractArray, F, P, nT}
    fbar!::F  # RHS callable
    prop!::P # linear propagator callable
    y::T  # Solution at current t
    yn::T  # Solution at t+dt
    yi::T  # Interpolant array (see interpolate())
    yerr::T  # solution error estimate (from embedded RK)
    ks::NTuple{7, T}  # k values (intermediate solutions for Runge-Kutta method)
    t::Float64  # current time (propagation variable)
    tn::Float64  # next time
    dt::Float64  # time step
    dtn::Float64  # time step for next step
    rtol::Float64  # relative tolerance on error
    atol::Float64  # absolute tolerance on error
    safety::Float64  # safety factor for stepsize control
    max_dt::Float64  # maximum value for dt (default Inf)
    min_dt::Float64  # minimum value for dt (default 0)
    locextrap::Bool  # true if using local extrapolation
    ok::Bool  # true if current step was successful
    err::Float64  # error metric to be compared to tol
    errlast::Float64  # error of the most recent successful step
    norm::nT  # function to calculate error metric, defaults to RK45.weaknorm
end

function PreconStepper(f!, linop, y0, t, dt;
                       rtol=1e-6, atol=1e-10, safety=0.9, max_dt=Inf, min_dt=0,
                       locextrap=true, norm=weaknorm)
    prop! = make_prop!(linop, y0)
    fbar! = make_fbar!(f!, prop!, y0)
    k1 = similar(y0)
    fbar!(k1, y0, t, t)
    ks = (k1, similar(k1), similar(k1), similar(k1), similar(k1), similar(k1), similar(k1))
    yerr = similar(y0)

    return PreconStepper(fbar!, prop!, copy(y0), copy(y0), similar(y0), yerr, ks,
        float(t), float(t), float(dt), float(dt), float(rtol), float(atol), float(safety),
        float(max_dt), float(min_dt), locextrap, false, 0.0, 0.0, norm)
end

function step!(s)
    evaluate!(s)

    if s.locextrap
        s.yn .= s.y
        for jj = 1:7
            b5[jj] == 0 || (s.yn .+= s.dt*b5[jj].*s.ks[jj])
        end
    end

    fill!(s.yerr, 0)
    for ii = 1:7
        errest[ii] == 0 || (@. s.yerr += s.dt*s.ks[ii]*errest[ii])
    end
    s.err = s.norm(s.yerr, s.y, s.yn, s.rtol, s.atol)
    s.ok = s.err <= 1
    stepcontrolPI!(s)
    if s.ok
        s.tn = s.t + s.dt
        s.ks[1] .= s.ks[end]
    else
        s.yn .= s.y
    end
    prop!_maybe(s) # propagate to new time to pass correct solution to stepfun
    return s.ok
end

function evaluate!(s::Stepper)
    # Set new time and stepsize values -- this happens at the beginning because
    # the interpolant still requires the old values after the step has finished
    s.dt = s.dtn
    s.t = s.tn
    s.y .= s.yn
    for ii = 1:6
        s.yn .= s.y
        for jj = 1:ii
            B[ii][jj] == 0 || (s.yn .+= s.dt*B[ii][jj].*s.ks[jj])
        end
        s.f!(s.ks[ii+1], s.yn, s.t+nodes[ii]*s.dt)
    end
end

function evaluate!(s::PreconStepper)
    # Set new time and stepsize values -- this happens at the beginning because
    # the interpolant still requires the old values after the step has finished
    s.y .= s.yn
    s.prop!(s.ks[1], s.t, s.tn)
    s.dt = s.dtn
    s.t = s.tn
    for ii = 1:6
        s.yn .= s.y
        for jj = 1:ii
            B[ii][jj] == 0 || (s.yn .+= s.dt*B[ii][jj].*s.ks[jj])
        end
        s.fbar!(s.ks[ii+1], s.yn, s.t, s.t+nodes[ii]*s.dt)
    end
end

prop!_maybe(s::PreconStepper) = s.prop!(s.yn, s.t, s.tn)
prop!_maybe(s) = nothing

"Interpolate solution, aka dense output."
function interpolate(s::Stepper, ti::Float64)
    if ti > s.tn
        error("Attempting to extrapolate!")
    end
    if ti == s.t
        return s.y
    elseif ti == s.tn
        return s.yn
    end
    σ = (ti - s.t)/s.dt
    σp = map(p -> σ^p, range(1, stop=4))
    b = sum(σp.*interpC, dims=1)
    fill!(s.yi, 0)
    for ii = 1:7
         s.yi .+= s.ks[ii].*b[ii]
    end
    return @. s.y + s.dt.*s.yi
end

"Interpolate solution, aka dense output."
function interpolate(s::PreconStepper, ti::Float64)
    if ti > s.tn
        error("Attempting to extrapolate!")
    end
    if ti == s.t
        return s.y
    elseif ti == s.tn
        return s.yn
    end
    σ = (ti - s.t)/s.dt
    σp = map(p -> σ^p, range(1, stop=4))
    b = sum(σp.*interpC, dims=1)
    fill!(s.yi, 0)
    for ii = 1:7
         s.yi .+= s.ks[ii].*b[ii]
    end
    out =  @. s.y + s.dt.*s.yi
    s.prop!(out, s.t, ti)
    return out
end

"Make propagator for the case of constant linear operator"
function make_prop!(linop::AbstractArray, y0)
    prop! = let linop=linop
        function prop!(y, t1, t2, bwd=false)
            if bwd
                @. y *= exp(linop*(t1-t2))
            else
                @. y *= exp(linop*(t2-t1))
            end
        end
    end
end

"Make propagator for the case of non-constant linear operator"
function make_prop!(linop!, y0)
    linop_int = similar(y0)
    lastt2 = [typemin(Float64)]
    function prop!(y, t1, t2, bwd=false)
        #= linop is always evaluated at later time, even for backward propagation
            therefore, linop is often evaluated at the same t2 twice in a row=#
        (lastt2[1] != t2) && linop!(linop_int, t2)
        lastt2[1] = t2
        dt = bwd ? (t1-t2) : (t2-t1)
        @. y *= exp(linop_int*dt)
    end
    return prop!
end

"Make closure for the pre-conditioned RHS function."
function make_fbar!(f!, prop!, y0)
    y = similar(y0)
    fbar! = let f! = f!, prop! = prop!, y=y
        function fbar!(out, ybar, t1, t2)
            y .= ybar
            prop!(y, t1, t2) # propagate to t2
            f!(out, y, t2) # evaluate RHS function
            prop!(out, t1, t2, true) # propagate back to t1
        end
    end
end

"Max-ish norm (from Dane Austin's code, no idea where he got it from)."
function maxnorm(yerr, y, yn, rtol, atol)
    maxerr = 0
    maxy = 0
    for ii in eachindex(yerr)
        maxerr = max(maxerr, abs(yerr[ii]))
        maxy = max(maxy, max(abs(y[ii]), abs(yn[ii])))
    end
    return maxerr/(atol + rtol*maxy)
end

"Alternative form of max-ish norm."
function maxnorm_ratio(yerr, y, yn, rtol, atol)
    m = 0
    for ii in eachindex(yerr)
        den = atol + rtol*max(abs(y[ii]), abs(yn[ii]))
        m = max(abs(yerr[ii])/den, m)
    end
    return m
end

"Semi-norm as used in DifferentialEquations.jl, see Hairer, Solving Ordinary Differential
Equations: Nonstiff Problems, eq. (4.11) (p.168 of the second revised edition)."
function normnorm(yerr, y, yn, rtol, atol)
    s = 0
    for ii in eachindex(yerr)
        s += abs2(yerr[ii]/(atol + rtol*max(abs(y[ii]), abs(yn[ii]))))
    end
    sqrt(s/length(yerr))
end

"'Weak' norm as used in fnfep."
function weaknorm(yerr, y, yn, rtol, atol)
    sy = 0
    syn = 0
    syerr = 0
    for ii in eachindex(yerr)
        sy += abs2(y[ii])
        syn += abs2(yn[ii])
        syerr += abs2(yerr[ii])
    end
    errwt = max(max(sqrt(sy), sqrt(syn)), atol)
    return sqrt(syerr)/rtol/errwt
end

"Simple proportional error controller, see e.g. Hairer eq. (4.13)."
function stepcontrolP!(s)
    if s.ok
        # if error is zero, there is no nonlinearity: increase step size by a lot
        s.dtn = s.err == 0 ? 1.5*s.dt : s.dt * min(5, s.safety*(s.err)^(-1/5))
    else
        if !isfinite(s.err) # check for NaN or Inf
            s.dtn = s.dt/2  # if we have one then we're in big trouble so halve the step size
        else
            s.dtn = s.dt * max(0.1, s.safety*(s.err)^(-1/5))
        end
    end
    steplims!(s)
end

"Proportional-integral error controller, aka Lund stabilisation.
See G. Söderlind and L. Wang, J. Comput. Appl. Math. 185, 225 (2006).
"
function stepcontrolPI!(s)
    β1 = 3/5 / 5
    β2 = -1/5 / 5
    ε = 0.8
    if s.ok
        s.errlast == 0 && (s.errlast = s.err) # if last error is zero, use current error instead
        if s.err == 0
            fac = 1.5 # zero error means no nonlinearity: increase step size by a lot
        else
            fac = s.safety * (ε/s.err)^β1 * (ε/s.errlast)^β2
        end
        # (0.99 <= fac <= 1.01) && (fac = 1.0)
        s.dtn = fac * s.dt
        s.errlast = s.err
    else
        if !isfinite(s.err) # check for NaN or Inf
            s.dtn = s.dt/2  # if we have one then we're in big trouble so halve the step size
        else
            s.dtn = s.dt * max(0.1, s.safety*(s.err)^(-1/5))
        end
    end
    steplims!(s)
end

"Apply user-defined limits on step size."
function steplims!(s)
    if s.dtn > s.max_dt
        s.dtn = s.max_dt
    elseif s.dtn < s.min_dt
        s.dtn = s.min_dt
        s.ok = true
    end
end

function donothing!(y, z, dz, interpolant)
end

end


text:
using POMDPPolicies
using DataFrames

@testset "evaluation.jl" begin
    mg, me = simple_mdps()
    pp = PMDPs.pp(mg)
    # trace = PMDPs.simulate_trace(mg, MersenneTwister(123))

    # Three non-zero requests with budget 10
    trace = simple_short_trace(pp)

    hrpl = PMDPs.HistoryReplayer(mg, trace)
    hrec = HistoryRecorder(max_steps = PMDPs.selling_period_end(mg), rng = MersenneTwister(4321))

    # test replay()
    reject =  FunctionPolicy(x->1000.)
    hᵣ = PMDPs.replay(hrpl, reject, MersenneTwister(123))
    @test isa(hᵣ, SimHistory)
    @test [(s = r.s, info=r.info) for r in hᵣ] == trace.hist

    # test get_metrics()
    (r,u,nₛ, nᵣ) = PMDPs.get_metrics(hrpl, hᵣ)
    @test (r,u,nₛ, nᵣ) == (0.,0,0,3)

    accept = FunctionPolicy(x->5.)
    hₐ = PMDPs.replay(hrpl, accept, MersenneTwister(123))
    @test PMDPs.get_metrics(hrpl, hₐ) == (r=15., u=3, nₛ=3, nᵣ=3)

    # test eval()
    policies = (vi = PMDPs.get_VI_policy(me),
                mcts = PMDPs.get_MCTS_planner(mg),
                hind = PMDPs.LP.get_MILP_hindsight_policy(mg, trace),
                flat =  PMDPs.get_flatrate_policy(mg, [trace, trace])
                )

    # Every policy should give optimal allocation for provided trace
    for policy in policies
        @test PMDPs.get_metrics(hrpl, PMDPs.replay(hrpl, policy, MersenneTwister(1))) == (r=30., u=3, nₛ=3, nᵣ=3)
    end

    metrics = PMDPs.eval_policy(mg, trace, policies, MersenneTwister(1))
    @test isa(metrics, DataFrame)
    @test size(metrics)==(length(policies),8)

    metrics_all = PMDPs.eval_policy(mg, [trace, trace], policies, MersenneTwister(1))
    @test isa(metrics_all, DataFrame)
    @test size(metrics_all)==(length(policies)*2,8)
end


text:
@testset "log" begin
    all_types = (RotMatrix{3}, AngleAxis, RotationVec,
                 UnitQuaternion, RodriguesParam, MRP,
                 RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
                 RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ,
                 RotX, RotY, RotZ,
                 RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY,
                 RotMatrix{2}, Angle2d)

    @testset "$(T)" for T in all_types, F in (one, rand)
        R = F(T)
        @test R ≈ exp(log(R))
        @test log(R) isa SMatrix
    end
end


text:
# overloadable for other types that might want to offer similar interactions
function interactions end

interactions(ax::Axis) = ax.interactions
interactions(ax3::Axis3) = ax3.interactions

"""
    register_interaction!(parent, name::Symbol, interaction)

Register `interaction` with `parent` under the name `name`.
The parent will call `process_interaction(interaction, event, parent)`
whenever suitable events happen.

The interaction can be removed with `deregister_interaction!` or temporarily
toggled with `activate_interaction!` / `deactivate_interaction!`.
"""
function register_interaction!(parent, name::Symbol, interaction)
    haskey(interactions(parent), name) && error("Interaction $name already exists.")
    registration_setup!(parent, interaction)
    push!(interactions(parent), name => (true, interaction))
    return interaction
end

"""
    register_interaction!(interaction::Function, parent, name::Symbol)

Register `interaction` with `parent` under the name `name`.
The parent will call `process_interaction(interaction, event, parent)`
whenever suitable events happen.
This form with the first `Function` argument is especially intended for `do` syntax.

The interaction can be removed with `deregister_interaction!` or temporarily
toggled with `activate_interaction!` / `deactivate_interaction!`.
"""
function register_interaction!(interaction::Function, parent, name::Symbol)
    haskey(interactions(parent), name) && error("Interaction $name already exists.")
    registration_setup!(parent, interaction)
    push!(interactions(parent), name => (true, interaction))
    return interaction
end

"""
    deregister_interaction!(parent, name::Symbol)

Deregister the interaction named `name` registered in `parent`.
"""
function deregister_interaction!(parent, name::Symbol)
    !haskey(interactions(parent), name) && error("Interaction $name does not exist.")
    _, interaction = interactions(parent)[name]

    deregistration_cleanup!(parent, interaction)
    pop!(interactions(parent), name)
    return interaction
end

function registration_setup!(parent, interaction)
    # do nothing in the default case
end

function deregistration_cleanup!(parent, interaction)
    # do nothing in the default case
end

"""
    activate_interaction!(parent, name::Symbol)

Activate the interaction named `name` registered in `parent`.
"""
function activate_interaction!(parent, name::Symbol)
    !haskey(interactions(parent), name) && error("Interaction $name does not exist.")
    interactions(parent)[name] = (true, interactions(parent)[name][2])
    return nothing
end

"""
    deactivate_interaction!(parent, name::Symbol)

Deactivate the interaction named `name` registered in `parent`.
It can be reactivated with `activate_interaction!`.
"""
function deactivate_interaction!(parent, name::Symbol)
    !haskey(interactions(parent), name) && error("Interaction $name does not exist.")
    interactions(parent)[name] = (false, interactions(parent)[name][2])
    return nothing
end

function process_interaction(@nospecialize args...)
    # do nothing in the default case
    return Consume(false)
end

# a generic fallback for functions to have one really simple path to getting interactivity
# without needing to define a special type first
function process_interaction(f::Function, event, parent)
    # in case f is only defined for a specific type of event
    if applicable(f, event, parent)
        # TODO this is deprecation code, make this just `return f(event, parent)` eventually
        return f(event, parent)
    end
    return Consume(false)
end

############################################################################
#                            Axis interactions                            #
############################################################################

function _chosen_limits(rz, ax)

    r = positivize(FRect2D(rz.from, rz.to .- rz.from))
    lims = ax.finallimits[]
    # restrict to y change
    if rz.restrict_x || !ax.xrectzoom[]
        r = FRect2D(lims.origin[1], r.origin[2], widths(lims)[1], widths(r)[2])
    end
    # restrict to x change
    if rz.restrict_y || !ax.yrectzoom[]
        r = FRect2D(r.origin[1], lims.origin[2], widths(r)[1], widths(lims)[2])
    end
    return r
end

function _selection_vertices(outer, inner)
    _clamp(p, plow, phigh) = Point2f0(clamp(p[1], plow[1], phigh[1]), clamp(p[2], plow[2], phigh[2]))

    outer = positivize(outer)
    inner = positivize(inner)

    obl = bottomleft(outer)
    obr = bottomright(outer)
    otl = topleft(outer)
    otr = topright(outer)

    ibl = _clamp(bottomleft(inner), obl, otr)
    ibr = _clamp(bottomright(inner), obl, otr)
    itl = _clamp(topleft(inner), obl, otr)
    itr = _clamp(topright(inner), obl, otr)

    return [obl, obr, otr, otl, ibl, ibr, itr, itl]
end

function process_interaction(r::RectangleZoom, event::MouseEvent, ax::Axis)

    # TODO: actually, the data from the mouse event should be transformed already
    # but the problem is that these mouse events are generated all the time
    # and outside of log axes, you would quickly run into domain errors
    transf = Makie.transform_func(ax)
    inv_transf = Makie.inverse_transform(transf)

    if isnothing(inv_transf)
        @warn "Can't rectangle zoom without inverse transform"
        # TODO, what can we do without inverse?
        return Consume(false)
    end

    if event.type === MouseEventTypes.leftdragstart
        data = Makie.apply_transform(inv_transf, event.data)
        prev_data = Makie.apply_transform(inv_transf, event.prev_data)

        r.from = prev_data
        r.to = data
        r.rectnode[] = _chosen_limits(r, ax)
        r.active[] = true
        return Consume(true)

    elseif event.type === MouseEventTypes.leftdrag
        # clamp mouse data to shown limits
        rect = Makie.apply_transform(transf, ax.finallimits[])
        data = Makie.apply_transform(inv_transf, rectclamp(event.data, rect))

        r.to = data
        r.rectnode[] = _chosen_limits(r, ax)
        return Consume(true)

    elseif event.type === MouseEventTypes.leftdragstop
        try
            r.callback(r.rectnode[])
        catch e
            @warn "error in rectangle zoom" exception=e
        end
        r.active[] = false
        return Consume(true)
    end

    return Consume(false)
end

function rectclamp(p::Point, r::Rect)
    map(p, minimum(r), maximum(r)) do pp, mi, ma
        clamp(pp, mi, ma)
    end |> Point
end

function process_interaction(r::RectangleZoom, event::KeysEvent, ax::Axis)
    r.restrict_y = Keyboard.x in event.keys
    r.restrict_x = Keyboard.y in event.keys
    r.active[] || return Consume(false)

    r.rectnode[] = _chosen_limits(r, ax)
    return Consume(true)
end

function positivize(r::FRect2D)
    negwidths = r.widths .< 0
    newori = ifelse.(negwidths, r.origin .+ r.widths, r.origin)
    newwidths = ifelse.(negwidths, -r.widths, r.widths)
    FRect2D(newori, newwidths)
end

function process_interaction(l::LimitReset, event::MouseEvent, ax::Axis)

    if event.type === MouseEventTypes.leftclick
        if ispressed(ax.scene, Keyboard.left_control)
            if ispressed(ax.scene, Keyboard.left_shift)
                autolimits!(ax)
            else
                reset_limits!(ax)
            end
            return Consume(true)
        end
    end

    return Consume(false)
end

function process_interaction(s::ScrollZoom, event::ScrollEvent, ax::Axis)
    # use vertical zoom
    zoom = event.y

    tlimits = ax.targetlimits
    xzoomlock = ax.xzoomlock
    yzoomlock = ax.yzoomlock
    xzoomkey = ax.xzoomkey
    yzoomkey = ax.yzoomkey

    scene = ax.scene
    e = events(scene)
    cam = camera(scene)

    if zoom != 0
        pa = pixelarea(scene)[]

        z = (1f0 - s.speed)^zoom

        mp_axscene = Vec4f0((e.mouseposition[] .- pa.origin)..., 0, 1)

        # first to normal -1..1 space
        mp_axfraction =  (cam.pixel_space[] * mp_axscene)[1:2] .*
            # now to 1..-1 if an axis is reversed to correct zoom point
            (-2 .* ((ax.xreversed[], ax.yreversed[])) .+ 1) .*
            # now to 0..1
            0.5 .+ 0.5

        xscale = ax.xscale[]
        yscale = ax.yscale[]

        transf = (xscale, yscale)
        tlimits_trans = Makie.apply_transform(transf, tlimits[])

        xorigin = tlimits_trans.origin[1]
        yorigin = tlimits_trans.origin[2]

        xwidth = tlimits_trans.widths[1]
        ywidth = tlimits_trans.widths[2]

        newxwidth = xzoomlock[] ? xwidth : xwidth * z
        newywidth = yzoomlock[] ? ywidth : ywidth * z

        newxorigin = xzoomlock[] ? xorigin : xorigin + mp_axfraction[1] * (xwidth - newxwidth)
        newyorigin = yzoomlock[] ? yorigin : yorigin + mp_axfraction[2] * (ywidth - newywidth)

        timed_ticklabelspace_reset(ax, s.reset_timer, s.prev_xticklabelspace, s.prev_yticklabelspace, s.reset_delay)

        newrect_trans = if ispressed(scene, xzoomkey[])
            FRect(newxorigin, yorigin, newxwidth, ywidth)
        elseif ispressed(scene, yzoomkey[])
            FRect(xorigin, newyorigin, xwidth, newywidth)
        else
            FRect(newxorigin, newyorigin, newxwidth, newywidth)
        end

        inv_transf = Makie.inverse_transform(transf)
        tlimits[] = Makie.apply_transform(inv_transf, newrect_trans)
    end

    # NOTE this might be problematic if if we add scrolling to something like Menu
    return Consume(true)
end

function process_interaction(dp::DragPan, event::MouseEvent, ax)

    if event.type !== MouseEventTypes.rightdrag
        return Consume(false)
    end

    tlimits = ax.targetlimits
    xpanlock = ax.xpanlock
    ypanlock = ax.ypanlock
    xpankey = ax.xpankey
    ypankey = ax.ypankey
    panbutton = ax.panbutton

    scene = ax.scene
    cam = camera(scene)
    pa = pixelarea(scene)[]

    mp_axscene = Vec4f0((event.px .- pa.origin)..., 0, 1)
    mp_axscene_prev = Vec4f0((event.prev_px .- pa.origin)..., 0, 1)

    mp_axfraction, mp_axfraction_prev = map((mp_axscene, mp_axscene_prev)) do mp
        # first to normal -1..1 space
        (cam.pixel_space[] * mp)[1:2] .*
        # now to 1..-1 if an axis is reversed to correct zoom point
        (-2 .* ((ax.xreversed[], ax.yreversed[])) .+ 1) .*
        # now to 0..1
        0.5 .+ 0.5
    end

    xscale = ax.xscale[]
    yscale = ax.yscale[]

    transf = (xscale, yscale)
    tlimits_trans = Makie.apply_transform(transf, tlimits[])

    movement_frac = mp_axfraction .- mp_axfraction_prev

    xscale = ax.xscale[]
    yscale = ax.yscale[]

    transf = (xscale, yscale)
    tlimits_trans = Makie.apply_transform(transf, tlimits[])

    xori, yori = tlimits_trans.origin .- movement_frac .* widths(tlimits_trans)

    if xpanlock[] || ispressed(scene, ypankey[])
        xori = tlimits_trans.origin[1]
    end

    if ypanlock[] || ispressed(scene, xpankey[])
        yori = tlimits_trans.origin[2]
    end

    timed_ticklabelspace_reset(ax, dp.reset_timer, dp.prev_xticklabelspace, dp.prev_yticklabelspace, dp.reset_delay)

    inv_transf = Makie.inverse_transform(transf)
    newrect_trans = FRect(Vec2f0(xori, yori), widths(tlimits_trans))
    tlimits[] = Makie.apply_transform(inv_transf, newrect_trans)

    return Consume(true)
end

function process_interaction(dr::DragRotate, event::MouseEvent, ax3d)
    if event.type !== MouseEventTypes.leftdrag
        return Consume(false)
    end

    dpx = event.px - event.prev_px

    ax3d.azimuth[] += -dpx[1] * 0.01
    ax3d.elevation[] = clamp(ax3d.elevation[] - dpx[2] * 0.01, -pi/2 + 0.001, pi/2 - 0.001)

    return Consume(true)
end


text: include("unknowns.jl")


text:
using LinearAlgebra
points = [[1, 2, 3], [3, 4, 5], [5, 6, 7]]
qpoint = [3, 6, 5]
argmin([norm(p - qpoint) for p in points])

max_stars_repo_name: fvcalderan/Langs-Comparisons

max_issues_repo_name: fvcalderan/Langs-Comparisons

max_forks_repo_name: fvcalderan/Langs-Comparisons