text:
using Test
import AoC2019

testdata1 =
"""#########
#b.A.@.a#
#########"""

testdata2 =
"""########################
#f.D.E.e.C.b.A.@.a.B.c.#
######################.#
#d.....................#
########################"""
testdata3 =
"""########################
#...............b.C.D.f#
#.######################
#.....@.a.B.c.d.A.e.F.g#
########################"""
testdata4 =
"""#################
#i.G..c...e..H.p#
########.########
#j.A..b...f..D.o#
########@########
#k.E..a...g..B.n#
########.########
#l.F..d...h..C.m#
#################"""
testdata5 =
"""########################
#@..............ac.GI.b#
###d#e#f################
###A#B#C################
###g#h#i################
########################"""

testdata6 =
"""#######
#a.#Cd#
##@#@##
#######
##@#@##
#cB#Ab#
#######"""

function can_pass(x, keys = [])
  if x == '.'
    return true
  elseif x == '@'
    return true
  elseif islowercase(x)
    return true
  elseif lowercase(x) in keys
    return true
  else
    return false
  end
end

function search_labyrinth(labyrinth_arr, start_pos, keys = [])
  paths = AoC2019.bfs_depth(labyrinth_arr, start_pos, x -> can_pass(x, keys))
  key_filter = map(x -> islowercase(x) && x ∉ keys , labyrinth_arr)
  paths .*= key_filter
  return paths
end

function generate_combinations(items, num_samples)
  combinations = collect(Iterators.product((items for i=1:num_samples)...))
  combinations = filter((x->allunique(x)), combinations)
  combinations = map(x->Tuple(sort!(collect(x))), combinations)
  return unique(combinations)
end

function do_labyrinth(arr)
  queue = []
  #start_pos = findall(x -> x == '@', arr)
  start_pos = findfirst(x -> x == '@', arr)
  num_robots = length(start_pos)
  #println(start_pos)
  number_of_keys = length(filter(x->islowercase(x), arr))
  println("Number of keys: $number_of_keys")
  stack = [(start_pos=start_pos, found_keys = [], steps_done = 0)]
  counter = 0
  #combinations = generate_combinations('a':'z', num_robots)
  #known_paths = Dict(x => Set() for x in combinations)
  known_paths = Dict(x => Set() for x in 'a':'z')
  max_length = 0
  while length(stack) > 0
    @label start
    (start_pos, found_keys, steps_done) = pop!(stack)
    if length(found_keys) > max_length
      max_length = length(found_keys)
    end
    if length(found_keys) > 0
      dict_entry = sort!(found_keys[1:end-1])
      for entry in known_paths[found_keys[end]]
        if intersect(entry, dict_entry) == dict_entry
          @goto start
        end
      end
      push!(known_paths[found_keys[end]], dict_entry)
    end
    #println(length(known_paths))
    if length(found_keys) == number_of_keys
      return steps_done
    end
    paths = search_labyrinth(arr, start_pos, found_keys)
    indices = findall(x -> x != 0, paths)
    for index in indices
      new_start_pos = index
      new_found_keys = push!(copy(found_keys), arr[index])
      new_steps_done = steps_done + paths[index]
      push!(stack, (start_pos=new_start_pos, found_keys=new_found_keys, steps_done=new_steps_done))
    end
    sort!(stack, by=x->x.steps_done, rev=true)
    counter += 1
    println(max_length)
  end
end

function multiply_robots!(arr)
  start_pos = Tuple(findfirst(x -> x == '@', arr))
  arr[start_pos...] = '#'
  for movement in [(-1,0), (1,0), (0,-1), (0,1)]
    new_pos = start_pos .+ movement
    arr[new_pos...] = '#'
  end
  for movement in [(-1,-1), (1,1), (1,-1), (-1,1)]
    new_pos = start_pos .+ movement
    arr[new_pos...] = '@'
  end
  return arr
end

function do_labyrinth_multiple(arr)
  queue = []
  start_pos = map(x -> Tuple(x), findall(x -> x == '@', arr))
  #start_pos = findfirst(x -> x == '@', arr)
  num_robots = length(start_pos)
  #println(start_pos)
  number_of_keys = length(filter(x->islowercase(x), arr))
  println("Number of keys: $number_of_keys")
  stack = [(start_pos=start_pos, found_keys = [], steps_done = 0)]
  counter = 0
  #combinations = generate_combinations('a':'z', num_robots)
  #known_paths = Dict(x => Set() for x in combinations)
  known_paths = Dict(x => Set() for x in 'a':'z')
  max_length = 0
  while length(stack) > 0
    @label start
    (start_pos, found_keys, steps_done) = pop!(stack)
    if length(found_keys) > max_length
      max_length = length(found_keys)
    end
    if length(found_keys) > 0
      dict_entry = sort!(found_keys[1:end-1])
      for entry in known_paths[found_keys[end]]
        if intersect(entry, dict_entry) == dict_entry
          @goto start
        end
      end
      push!(known_paths[found_keys[end]], dict_entry)
    end
    #println(length(known_paths))
    if length(found_keys) == number_of_keys
      return steps_done
    end
    for (i, robot_pos) in enumerate(start_pos)
      paths = search_labyrinth(arr, robot_pos, found_keys)
      indices = findall(x -> x != 0, paths)
      for index in indices
        #new_start_pos = index
        new_start_pos = vcat(start_pos[1:i-1], start_pos[i+1:end], [index])
        new_found_keys = push!(copy(found_keys), arr[index])
        new_steps_done = steps_done + paths[index]
        push!(stack, (start_pos=new_start_pos, found_keys=new_found_keys, steps_done=new_steps_done))
      end
    end
    sort!(stack, by=x->x.steps_done, rev=true)
    println(max_length)
  end
end
@test do_labyrinth(AoC2019.text2arr(testdata1)) == 8
@test do_labyrinth(AoC2019.text2arr(testdata2)) == 86
@test do_labyrinth(AoC2019.text2arr(testdata3)) == 132
@test do_labyrinth(AoC2019.text2arr(testdata4)) == 136
@test do_labyrinth(AoC2019.text2arr(testdata5)) == 81
#steps_done = do_labyrinth(labyrinth_arr)
#steps_done = do_labyrinth_multiple(AoC2019.text2arr(testdata6))

input_data = read(open("input18.txt"), String)
labyrinth_arr = AoC2019.text2arr(input_data)

labyrinth_arr = multiply_robots!(labyrinth_arr)
steps_done = do_labyrinth_multiple(labyrinth_arr)
println(steps_done)

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
#  This returns data that follow a logit model with the
#  supplied parameters. n is the number of observations.
#  Note how the binary 0/1 variable is generated so that
#  it really follows the logit model
function LogitDGP(n, theta)
	k = size(theta,1)
	x = ones(n,1)
	if k>1
		x = [x  randn(n,k-1)]
	end
	y = Float64.((1.0 ./ (1.0 .+ exp.(-x*theta)) .> rand(n,1)))
    return y, x
end

ext: jl

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
@testset "Main" begin
  @testset "Simulate with trace" begin
    n = 3
    ε = .1
    reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:n]
    instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

    Random.seed!(1)
    n_rounds = 2
    s, t = simulate(instance, ThompsonSampling(), n_rounds, with_trace=true)
    @test s.round == n_rounds

    @test length(t.states) == n_rounds
    @test length(t.arms) == n_rounds
    @test length(t.reward) == n_rounds

    @test t.states[n_rounds].round == s.round
    @test t.states[n_rounds].regret == s.regret
    @test t.states[n_rounds].reward == s.reward
    @test t.states[n_rounds].arm_counts == s.arm_counts
    @test t.states[n_rounds].arm_reward == s.arm_reward
    @test t.states[n_rounds].arm_average_reward == s.arm_average_reward
  end
end


text:
using VariableFrequencyOPF
using Documenter

makedocs(;
    modules=[VariableFrequencyOPF],
    authors="David Sehloff",
    repo="https://github.com/WISPO-POP/VariableFrequencyOPF.jl/blob/{commit}{path}#L{line}",
    sitename="VariableFrequencyOPF.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://WISPO-POP.github.io/VariableFrequencyOPF.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => [
            "Variable Frequency OPF" => "opf_examples.md",
            "Comparing Upgrades" => "comparison_examples.md",
            "Analyzing Frequency Dependence" => "frequency_examples.md",
            "Analyzing LFAC Upgrades Across Multiple Time Points" => "timeseries.md"
            ],
        "Functions" => "functions.md"
    ],
)

deploydocs(;
    repo="github.com/WISPO-POP/VariableFrequencyOPF.jl.git",
)


text:
"""Will load a collection of 90 problems from a passed filename.
Files must be in the
http://people.brunel.ac.uk/~mastjjb/jeb/orlib/mdmkpinfo.html format."""
function parse_file(filename::String, dataset_num::Int, uid=0)
    f = open(filename)

    problems::Vector{MDMKP_Prob} = []

    #the very first item in the array is the amount of problems found in the
    #file.
    amount_of_problems = next_line(f)[1]

	instance_num = 0
    #so now for every problem:
    for problem in 1:amount_of_problems
		instance_num += 1
        n, m = next_line(f)
        lower_than_values::Vector{Vector{Int}} = []
        for i in 1:m
            push!(lower_than_values, next_line(f))
        end
        lower_than_constraints::Vector{Int} =  next_line(f)
        greater_than_values::Vector{Vector{Int}} = []
        for i in 1:m
            push!(greater_than_values, next_line(f))
        end
        greater_than_constraints = next_line(f)
        cost_coefficient_values::Vector{Vector{Int}} = []
        for i in 1:6
            push!(cost_coefficient_values, next_line(f))
        end

        upper_bounds::Vector{Tuple{Vector{Int},Int}} = []
        lower_bounds::Vector{Tuple{Vector{Int},Int}} = []

        for i in 1:m
            push!(lower_bounds, (greater_than_values[i], greater_than_constraints[i]))
            push!(upper_bounds, (lower_than_values[i], lower_than_constraints[i]))
        end

        q = [1, div(m, 2), m, 1, div(m, 2), m]
        for i in 1:6
			uid += 1
			id = Problem_ID(
				uid,
				dataset_num,
				instance_num,
				i,
				get_tightness(instance_num),
				length(cost_coefficient_values[i]),
				length(lower_bounds[1:q[i]]),
				length(upper_bounds),
				is_mixed(i),
				false,
                '?'
			)
            push!(problems, MDMKP_Prob(
                cost_coefficient_values[i],
                upper_bounds,
                lower_bounds[1:q[i]],
				id
            ))
        end
    end

    (problems, uid)
end

function next_line(file::IOStream)
    return parse_line(readline(file))
end

function parse_line(line)
    return map(parse_int, split(line))
end

function parse_int(x)
    return parse(Int, x)
end

function get_tightness(instance)
	[.25, .5, .75][Int(ceil(instance/5))]
end

function is_mixed(case)
	case > 3
end

function load_folder(
		folder_path="MDMKP/benchmark_problems",
		filename="mdmkp_ct{ds}.txt", datasets=1:9)::Vector{MDMKP_Prob}
	collection = []
	id = 0
	for ds in datasets
		fn = replace(filename, "{ds}"=>"$ds")
		problems, id = parse_file(joinpath(folder_path, fn), ds, id)
		push!(collection, problems)
	end
	vcat(collection...)
end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
function add_watch_entry!(state::DebuggerState, cmd::AbstractString)
    expr = Base.parse_input_line(cmd)
    if isexpr(expr, :error) || isexpr(expr, :incomplete)
        printstyled(stderr, "ERROR: syntax: ", expr.args[1], "\n", color=Base.error_color())
        return false
    end
    isexpr(expr, :toplevel) && (expr = expr.args[end])
    push!(state.watch_list, expr)
    return true
end

function show_watch_list(io::IO, state::DebuggerState)
    frame = active_frame(state)
    for (i, expr) in enumerate(state.watch_list)
        vars = filter(v -> v.name != Symbol(""), JuliaInterpreter.locals(frame))
        eval_expr = Expr(:let,
            Expr(:block, map(x->Expr(:(=), x...), [(v.name, maybe_quote(v.value)) for v in vars])...),
            expr)
        errored = false
        res = try
            Core.eval(moduleof(frame), eval_expr)
        catch err
            errored = true
            err
        end
        expr_str = highlight_code(string(expr); context=io) # Should maybe use repr here in some cases (strings)
        if errored
            res_str = sprint((io, args...) -> printstyled(io, args...; color=Base.error_color()), res; context=io)
        else
            res_str = highlight_code(repr(res), context=io)
        end
        println(io, "$i] $(expr_str): $(res_str)")
    end
end

clear_watch_list!(state::DebuggerState) = empty!(state.watch_list)
function clear_watch_list!(state::DebuggerState, i::Int)
    if !checkbounds(Bool, state.watch_list, i)
        printstyled(stderr, "ERROR: watch entry $i does not exist\n"; color=Base.error_color())
        return
    end
    deleteat!(state.watch_list, i)
end


text:
# Move to LLVM.jl

function position_before_terminator!(builder::LLVM.Builder, bb::LLVM.BasicBlock)
    v = LLVM.API.LLVMGetBasicBlockTerminator(LLVM.blockref(bb))
    if v == C_NULL
        LLVM.position!(builder, bb)
    else
        i = LLVM.Instruction(LLVM.API.LLVMGetBasicBlockTerminator(LLVM.blockref(bb)))
        LLVM.position!(builder, i)
    end
end


text:
module GraphqlController
	using Genie.Router, Genie.Requests
	using Diana
	using JSON

	const GQLString = Dict("tipo"=>"String")

	function ping_mutation(root, args, ctx, info)
		print(root, args, ctx, info)
		return string("Hello ", args["id"])
	end

	schema = Dict(
		"query"=> "Query",
		"Query"   => Dict(
			#=
			Currently not working because of current implementation requires some efforts
			"__schema" => Dict(
				"queryType" => Dict(
					"name" => GQLString
				),
				"mutationType" => Dict(
					"name" => GQLString
				),
				"subscriptionType" => Dict(
					"name" => GQLString
				)
			),
			=#
			"hello"=>Dict("tipo"=>"String"),
		),
		"mutation" => "Mutation",
		"Mutation" => Dict(
			"ping" => Dict(
				"args" => Dict(
					"id" => "String"
				),
				"tipo"=>"String"
			)
		)
	)

	resolvers=Dict(
		"Query"=>Dict(
			"hello" => (root,args,ctx,info)->(return "Hello World!"),
		),
		"Mutation" => Dict(
			"ping"=>ping_mutation
		)
	)

	my_schema = Schema(schema, resolvers)

	function playground()
		serve_static_file("/playground.html")
	end

	function graphql()
		body = jsonpayload()
		# print(body, body["query"])
		out = "{}"
		try
			out = my_schema.execute(body["query"], operationName=body["operationName"], Variables=body["variables"])
		catch err
			if isa(err, Diana.GraphQLError)
				out = err.msg
			else
				stack = stacktrace(catch_backtrace())
				stack = map((st)->([st.file, st.line]), stack)
				out = json(Dict(
					"name"=>typeof(err),
					"stacktrace"=>stack,
					"info"=> sprint(showerror, err),
					))
			end
		end
		return out
	end

	function register()
		route("/graphql", playground)
		route("/graphql", method="POST", graphql)
	end

	export register
end


text:
#
# Functions to read the DCD file
#

using FortranFiles

#
# Functions to close and rewind a dcd file
#

closedcd(simulation :: Simulation) = close(simulation.FortranDCD)

rewinddcd(simulation :: Simulation) = rewind(simulation.FortranDCD)

#
# Function to read the header given the simulation data
#

function dcdheader(simulation :: Simulation)

  read_nframes, read_natoms, dcdaxis = dcdheader(simulation.FortranDCD)

  if read_natoms != simulation.natoms
    error("ERROR: Number of atoms in PSF file differs from that in DCD file.")
  end
  nframes = read_nframes
  println(" Number of frames in DCD file: ", nframes)

  if dcdaxis
    println(" DCD contains periodic cell data. ")
  else
    println(" WARNING: DCD does NOT contain periodic cell data! ")
  end

  return read_nframes, read_natoms, dcdaxis

end

#
# Function to read the header given only the file stream
#

function dcdheader(FortranDCD :: FortranFile)

    rewind(FortranDCD)

    # Read header
    IntVec = Vector{Int32}(undef,17)
    hdr, read_nframes, IntVec[1:8], ts, IntVec[9:17] = read(FortranDCD, FString{4}, Int32, (Int32,8), Float64, (Int32,9))
    dummyi, title = read(FortranDCD, Int32, FString{80})
    read_natoms = read(FortranDCD,Int32)

    # Check if dcd file contains axis information
    dcdaxis = false
    x = 0.
    try
      x = read(FortranDCD, [ Float32 for i in 1:read_natoms ])
    catch err
      dcdaxis = true
    end

    firstframe(FortranDCD)

    return read_nframes, read_natoms, dcdaxis

end

# Function that goes to the file point to read the first frame

firstframe(simulation :: Simulation) = firstframe(simulation.FortranDCD)
function firstframe(FortranDCD :: FortranFile)

    # rewind
    rewind(FortranDCD)

    # skip header
    read(FortranDCD)
    read(FortranDCD)
    read(FortranDCD)

end

#
# Function to read next frame
#

function nextframe(simulation :: Simulation;lastatom=0)

  if simulation.natoms == 0
    error(" ERROR: Set natoms before reading frames with M3GTools.init ")
  else
    lastatom = simulation.natoms
  end

  sides, x, y, z = nextframe(simulation.FortranDCD, simulation.dcdaxis, lastatom)
  return sides, x, y, z

end

function nextframe(FortranDCD :: FortranFile, dcdaxis, lastatom)

  if dcdaxis
    sides_read = read(FortranDCD,(Float64,6))
    sides = [ sides_read[1], sides_read[3], sides_read[6] ]
  else
    sides = zeros(3)
  end
  x = read(FortranDCD,(Float32,lastatom))
  y = read(FortranDCD,(Float32,lastatom))
  z = read(FortranDCD,(Float32,lastatom))

  return sides, x, y, z

end

#
# Function to read the number of frames if the header does not inform it
#

function getnframes(FortranDCD :: FortranFile, dcdaxis :: Bool )
  firstframe(FortranDCD)
  nframes = 0
  while true
    try
      if dcdaxis
        read(FortranDCD,Float64)
      end
      read(FortranDCD,Float32)
      read(FortranDCD,Float32)
      read(FortranDCD,Float32)
      nframes = nframes + 1
    catch
      firstframe(FortranDCD)
      return nframes
    end
  end
end

#
# Function that writes a dcd file
#

function writedcd(natoms,nframes,dcdaxis,sides,x,y,z;filename="Namdjl_DCDTEMP.dcd")

  dcdtemp = FortranFile(filename,"w";marker=RECMRK4B)

  IntVec = Vector{Int32}(undef,17)

  #hdr, read_nframes, IntVec[1:8], ts, IntVec[9:17] =
  #  read(dcdfile, FString{4}, Int32, (Int32,8), Float64, (Int32,9))
  #dummyi, title = read(dcdfile, Int32, FString{80})
  #read_natoms = read(dcdfile,Int32)

  hdr = "CORD"
  IntVec = [0 for i in 1:17]
  ts = 0.

#voltar: não funcionou
  write(dcdtemp,FString(4,hdr),Int32(nframes),IntVec[1:8],Float64(ts),IntVec[9:17])
  write(dcdtemp,Int32(1),FString(80,"DCDTEMP created by M3GTools"))
  write(dcdtemp,Int32(natoms))

  if nframes == 1
    if dcdaxis
      write(dcdtemp,Float64(sides[1]),0.,Float64(sides[2]),0.,0.,Float64(sides[3]))
    end
    write(dcdtemp,x)
    write(dcdtemp,y)
    write(dcdtemp,z)
  #else
  #  for i in 1:nframes
  #    if dcdaxis
  #      write(dcdtemp,Float64(sides[i,1]),0.,Float64(sides[i,2]),0.,0.,Float64(sides[i,3]))
  #    end
  #    write(dcdtemp,(x[i,j],j=1,natoms))
  #    write(dcdtemp,(y[i,j],j=1,natoms))
  #    write(dcdtemp,(z[i,j],j=1,natoms))
  #  end
  end
  close(dcdtemp)

end

#
# Function that writes a restart file for NAMD (.coor files)
#

function writecoor(natoms,x,y,z;filename="M3G.coor")
  coortemp = open(filename,"w")
  xyz = Vector{Float64}(undef,3*natoms)
  j = 1
  for i in 1:natoms
    xyz[j] = x[i]
    xyz[j+1] = y[i]
    xyz[j+2] = z[i]
    j = j + 3
  end
  write(coortemp,Int32(natoms),xyz)
  close(coortemp)
end


text:
# Graph to matrices

using Graphs
using Base.Test

# data

eds = [Edge(1,1,2), Edge(2,1,3), Edge(3,2,4), Edge(4,3,4)]

a0 = fill(false, 4, 4)
a0[1, 2] = a0[1, 3] = a0[2, 4] = a0[3, 4] = true

a0u = copy(a0)
a0u[2, 1] = a0u[3, 1] = a0u[4, 2] = a0u[4, 3] = true

a0_s = sparse(a0)
a0u_s = sparse(a0u)

# graphs

gd_elst = simple_edgelist(4, eds)
gu_elst = simple_edgelist(4, eds; is_directed=false)

gd_ilst = simple_inclist(4)
gu_ilst = simple_inclist(4, is_directed=false)

for e in eds
    add_edge!(gd_ilst, e.source, e.target)
    add_edge!(gu_ilst, e.source, e.target)
end

# adjacency matrix

@test adjacency_matrix(gd_elst) == a0
@test adjacency_matrix(gu_elst) == a0u

@test adjacency_matrix_sparse(gd_elst) == a0_s
@test adjacency_matrix_sparse(gu_elst) == a0u_s

@test adjacency_matrix(gd_ilst) == a0
@test adjacency_matrix(gu_ilst) == a0u

@test adjacency_matrix_sparse(gd_ilst) == a0_s
@test adjacency_matrix_sparse(gu_ilst) == a0u_s

# weight matrix

eweights = [1., 2., 3., 4.]
wm0  = [0. 1. 2. 0.; 0. 0. 0. 3.; 0. 0. 0. 4.; 0. 0. 0. 0.]
wm0u = [0. 1. 2. 0.; 1. 0. 0. 3.; 2. 0. 0. 4.; 0. 3. 4. 0.]
wm0_s = sparse(wm0)
wm0u_s = sparse(wm0u)

@test weight_matrix(gd_elst, eweights) == wm0
@test weight_matrix(gu_elst, eweights) == wm0u

@test weight_matrix_sparse(gd_elst, eweights) == wm0_s
@test weight_matrix_sparse(gu_elst, eweights) == wm0u_s

@test weight_matrix(gd_ilst, eweights) == wm0
@test weight_matrix(gu_ilst, eweights) == wm0u

@test weight_matrix_sparse(gd_ilst, eweights) == wm0_s
@test weight_matrix_sparse(gu_ilst, eweights) == wm0u_s

# distance matrix

dm0  = [0. 1. 2. Inf; Inf 0. Inf 3.; Inf Inf 0. 4.; Inf Inf Inf 0.]
dm0u = [0. 1. 2. Inf; 1. 0. Inf 3.; 2. Inf 0. 4.; Inf 3. 4. 0.]

@test distance_matrix(gd_elst, eweights) == dm0
@test distance_matrix(gu_elst, eweights) == dm0u

@test distance_matrix(gd_ilst, eweights) == dm0
@test distance_matrix(gu_ilst, eweights) == dm0u

# Laplacian matrix

L0 = [2. -1. -1. 0.; -1. 2. 0. -1.; -1. 0. 2. -1.; 0. -1. -1. 2.]
L0_s = sparse(L0)

@test laplacian_matrix(gu_elst) == L0
@test laplacian_matrix(gu_ilst) == L0

@test laplacian_matrix_sparse(gu_elst) == L0_s
@test laplacian_matrix_sparse(gu_ilst) == L0_s

Lw = [3. -1. -2. 0.; -1. 4. 0. -3.; -2. 0. 6. -4.; 0. -3. -4. 7.]
Lw_s = sparse(Lw)

@test laplacian_matrix(gu_elst, eweights) == Lw
@test laplacian_matrix(gu_ilst, eweights) == Lw

@test laplacian_matrix_sparse(gu_elst, eweights) == Lw_s
@test laplacian_matrix_sparse(gu_ilst, eweights) == Lw_s


text:
using JuLIP, QMMM2, SHIPs, PrettyTables, LinearAlgebra, Plots, JuLIPMaterials
using JuLIP.MLIPs
include(@__DIR__() * "/swqmmm.jl")

sw_eq = SWqmmm.sw_eq

##
r0 = rnn(:Si)
# domain radius of reference domain ("exact" solution)
Rmax = 75.0 * r0
# domain radius on approximation domains
RDOM = [5.0, 8.0, 12.0, 18.0, 26.0] * r0
RBUF = 2 * cutoff(sw_eq) *  ones(length(RDOM))

atmax, _ = SWqmmm.edge_config(Rmax, sw_eq)
AT, Iinmax, Iinnext = QMMM2.Solve.cluster_sequence(atmax, RDOM, RBUF)
UU, EE, Umax, Emax  =  QMMM2.Solve.solve_sequence(atmax, AT, Iinmax, Iinnext)
err2, errinf = QMMM2.Solve.errors(atmax, AT, Iinmax, UU, Umax)
Eerr = abs.(EE .- Emax)

##
plot(; xaxis = (:log, "domain radius [Å]"),
       yaxis = (:log, "errors"),
       title = "Domain Test, Si, SW",
       legend = :bottomleft )
plot!( RDOM, err2,   c=1, lw=3, m=:o, ms=8, label =  "energy-norm" )
plot!( RDOM, errinf, c=2, lw=3, m=:o, ms=8, label =  "max-norm (strains)" )
plot!( RDOM, Eerr,   c=3, lw=3, m=:o, ms=8, label =  "energy-diff. [eV]" )
t = RDOM[3:5]
plot!(t, 3*t.^(-1),  c=:black, lw=2, ls=:dot, label = "R^-1, R^-2")
plot!(t, 12*t.^(-2), c=:black, lw=2, ls=:dot, label = "")


text:
"""
    $(TYPEDEF) -> TiffFile

Wrap `io` with helper parameters to keep track of file attributes.

$(FIELDS)
"""
mutable struct TiffFile{O <: Unsigned, S <: Stream}
    """A unique identifier for this file"""
    uuid::Union{UUID, Nothing}

    """The relative path to this file"""
    filepath::String

    """The file stream"""
    io::S

    """Location of the first IFD in the file stream"""
    first_offset::Int

    """Whether this file has a different endianness than the host computer"""
    need_bswap::Bool
end

TiffFile{O}(s::Stream) where O <: Unsigned = TiffFile{O, typeof(s)}(nothing, "", s, -1, false)
TiffFile{O}(io::IO) where O <: Unsigned    = TiffFile{O}(getstream(format"TIFF", io))
TiffFile{O}() where O <: Unsigned          = TiffFile{O}(IOBuffer())

function Base.read(io::Stream, ::Type{TiffFile})
    seekstart(io)
    filepath = String(extract_filename(io))::String
    bs = need_bswap(io)
    offset_size = offsetsize(io)
    first_offset_raw = read(io, offset_size)
    first_offset = Int(bs ? bswap(first_offset_raw) : first_offset_raw)
    TiffFile{offset_size, typeof(io)}(nothing, filepath, io, first_offset, bs)
end

Base.read(io::IOStream, t::Type{TiffFile}) = read(getstream(format"TIFF", io, extract_filename(io)), t)

function Base.write(file::TiffFile{O}) where O
    seekstart(file.io)

    if ENDIAN_BOM == 0x04030201 #little endian
        write(file.io, "II")
    else
        write(file.io, "MM")
    end

    ifd_pos = 4 # position where the offset info is
    if O == UInt32
        write(file.io, UInt16(42)) # regular tiff
        write(file.io, UInt32(8)) # first offset is right after header
    elseif O == UInt64
        write(file.io, UInt16(43)) # big tiff
        write(file.io, UInt16(8)) # byte size of offsets
        write(file.io, UInt16(0)) # constant
        ifd_pos = position(file.io)
        write(file.io, UInt64(16)) # first offset is right after header
    else
        error("Unknown offset size")
    end
    ifd_pos
end

"""
    offset(file)

Returns the type of the TIFF file offset. For most TIFFs, it should be UInt32,
while for BigTIFFs it will be UInt64.
"""
offset(file::TiffFile{O}) where {O} = O

function Base.read(file::TiffFile{O}, ::Type{T}) where {O, T}
    value = read(file.io, T)
    file.need_bswap ? bswap(value) : value
end

function Base.read(file::TiffFile{O}, ::Type{String}) where O
    value = read(file.io, O)
    file.need_bswap ? bswap(value) : value
end

function Base.read!(file::TiffFile, arr::AbstractArray)
    _read!(stream(file.io), arr)
    if file.need_bswap
        arr .= bswap.(arr)
    end
end

function Base.read!(file::TiffFile, arr::SubArray{T,N,<:BitArray}) where {T, N}
    error("Strided bilevel TIFFs are not yet supported. Please open an issue against TiffImages.jl.")
end

function Base.read!(file::TiffFile, arr::BitArray)
    Bc = arr.chunks
    n = length(arr)
    nc = length(read!(file.io, Bc))
    if length(Bc) > 0 && Bc[end] & Base._msk_end(n) ≠ Bc[end]
        Bc[end] &= Base._msk_end(n) # ensure that the BitArray is not broken
    end
    for i in 1:nc
       Bc[i] = reversebits(Bc[i])
    end
    arr
end

Base.write(file::TiffFile, t) = write(file.io.io, t)::Int
Base.write(file::TiffFile, arr::AbstractVector{Any}) = write(file.io.io, Array{UInt8}(arr))::Int

Base.seek(file::TiffFile, n::Integer) = seek(file.io, n)

Base.bswap(x::Rational{T}) where {T} = Rational(bswap(x.num), bswap(x.den))

Base.IteratorSize(::TiffFile) = Base.SizeUnknown()
Base.eltype(::TiffFile{O}) where {O} = IFD{O}


text:
module DeconvolveDistribution

using FastGaussQuadrature
import Random
import Distributions
import Distributions: params, partype,
                      logpdf, pdf, cf,
                      logcdf, cdf,
                      rand,
                      UnivariateDistribution,
                      insupport,
                      Normal, Uniform, MixtureModel
using Statistics: var
using StatsBase: Weights
using LinearAlgebra: diagm, norm, dot
using RCall
import Optim

include("fourier.jl")
include("expospline.jl")
include("interface.jl")

export Fhat, fix_CDF!, decon,
       FourierDeconv, EfronDeconv

end # module


text:
@testset "Expected Loss" begin

    @testset "calculate expected loss" begin
        Random.seed!(1234)
        n = 5
        sampleA = [1, 1, 1, 0, 0]
        sampleB = [1, 0, 0, 1, 0]
        @test expectedloss(sampleA, sampleB, upliftloss, 5) == 0.2
        @test expectedloss(sampleB, sampleA, upliftloss, 5) == 0.4
    end

    @testset "Approximate expected loss from posteriors" begin
        Random.seed!(1234)
        n = 10000
        modelA = ConjugateBernoulli(50, 50)
        modelB = ConjugateBernoulli(60, 40)

        exploss = expectedloss(modelA, modelB, [:θ], lossfunc=upliftloss, numsamples=n)
        @test isapprox(exploss, 0.1, atol=0.01)

        exploss = expectedloss(modelB, modelA, [:θ], lossfunc=upliftloss, numsamples=n)
        @test isapprox(exploss, 0.0, atol=0.01)
    end

    @testset "Approximate expected loss from experiment" begin
        Random.seed!(1234)
        n = 10000
        modelA = ConjugateBernoulli(50, 50)
        modelB = ConjugateBernoulli(60, 40)

        exploss = expectedloss(modelA, modelB, lossfunc=upliftloss, numsamples=n)
        @test isapprox(exploss, 0.1, atol=0.01)

        stoppingrule = ExpectedLossThresh(0.001)
        experiment = ExperimentAB([modelA, modelB], stoppingrule)

        (modelnames, explosses) = expectedlosses(experiment, lossfunc=upliftloss, numsamples=n)

        @test modelnames == ["control", "variant 1"]
        @test isapprox(explosses[1], 0.1, atol=0.01)
        @test isapprox(explosses[2], 0.0, atol=0.01)
    end

end

@testset "ExperimentAB" begin
    @testset "ChainedModel with Bernoulli and LogNormal" begin
        Random.seed!(1234)
        eachrunnum = 100000
        truedistA1 = Bernoulli(0.05)
        truedistA2 = LogNormal(1.0, 1.0)
        truedistB1 = Bernoulli(0.05 * 1.05)
        truedistB2 = LogNormal(1.0, 1.0)

        numsims = 50000

        dataA1 = rand(truedistA1, eachrunnum)
        dataA2 = rand(truedistA2, eachrunnum)
        dataB1 = rand(truedistB1, eachrunnum)
        dataB2 = rand(truedistB2, eachrunnum)

        statsA1 = BetaStatistics(dataA1)
        statsA2 = LogNormalStatistics(dataA2)
        statsB1 = BetaStatistics(dataB1)
        statsB2 = LogNormalStatistics(dataB2)

        modelA = ChainedModel(
            [ConjugateBernoulli(1, 1), ConjugateLogNormal(0.0, 1.0, 0.001, 0.001)],
        )
        modelB = ChainedModel(
            [ConjugateBernoulli(1, 1), ConjugateLogNormal(0.0, 1.0, 0.001, 0.001)],
        )

        stoppingrule = ExpectedLossThresh(0.001)
        experiment = ExperimentAB([modelA, modelB], stoppingrule)
        update!(experiment, [[statsA1, statsA2], [statsB1, statsB2]])

        modelnames, explosses = expectedlosses(experiment)

        @testset "Control Group" begin
            model1_list = experiment.models["control"].models
            @test model1_list[1].dist.α == 4972
            @test model1_list[1].dist.β == 95030
            @test model1_list[2].dist.μ ≈ 1.0
            @test model1_list[2].dist.v ≈ 0
            @test model1_list[2].dist.α ≈ 50000
            @test isapprox(model1_list[2].dist.θ, 49373, atol=1.0)
        end

        @testset "Expected loss" begin
            @test modelnames == ["control", "variant 1"]
            @test explosses[1] ≈ 0.016
            @test explosses[2] ≈ 0
        end

        winner, _ = decide!(experiment)
        @test winner == "variant 1"
    end

end

@testset "Number of models" begin
    # Generate sample data
    n = 1000
    dataA = rand(Bernoulli(0.150), n)
    dataB = rand(Bernoulli(0.145), n)
    dataC = rand(Bernoulli(0.180), n)

    # Define the models
    modelA = ConjugateBernoulli(1, 1)
    modelB = ConjugateBernoulli(1, 1)
    modelC = ConjugateBernoulli(1, 1)

    stoppingrule = ProbabilityBeatAllThresh(0.99)
    @test_throws ArgumentError("Number of models needs to be equal to 2.") ExperimentAB([modelA, modelB, modelC], stoppingrule)
end

@testset "ExperiemntBF{NormalEffectSize}" begin
    @testset "Bayes factor calculation" begin
        # the example is taken from
        # https://statswithr.github.io/book/hypothesis-testing-with-normal-populations.html
        n0 = 32.7^2
        m0 = 0.5
        x̄ = 0.500177
        σ = 0.5
        n = 1.0449e8;
        σ0 = 1/sqrt(n0)

        model = NormalEffectSize(m0, σ0)
        stats = NormalStatistics(n=1.0449e8, meanx=x̄, sdx=σ)
        bf10 = bayesfactor(model, stats)
        @test isapprox(bf10, 2.2303, rtol=0.001)
    end

    @testset "Two samples with equal size and sd" begin
        n0 = 32.7^2
        m0 = 0.5 # null hypothesis: difference is 0.5
        σ = 0.5
        n = 1.0449e8;
        σ0 = 1/sqrt(n0)

        # 2*n so that the effective sample size will be `n`
        stats1 = NormalStatistics(n=2*n, meanx=1.000177, sdx=σ)
        stats2 = NormalStatistics(n=2*n, meanx=0.5, sdx=σ)
        stats = TwoNormalStatistics(stats1, stats2)

        model = NormalEffectSize(m0, σ0)
        bf10 = bayesfactor(model, stats)
        @test isapprox(bf10, 2.2303, rtol=0.001)
    end
end

@testset "ExperiemntBF{StudentTEffectSize}" begin
    @testset "Two-Sided null wins, thresh=5" begin
        normalstats = TwoNormalStatistics(
            NormalStatistics(meanx=28.8, sdx=13.5, n=133),
            NormalStatistics(meanx=30.6, sdx=14.3, n=867)
        )
        model = StudentTEffectSize(r=1.0)
        stats = StudentTStatistics(normalstats)
        bf10_model = bayesfactor(model, stats)

        thresh=5
        stoppingrule = TwoSidedBFThresh(thresh)
        experiment = ExperimentBF(model=model, rule=stoppingrule)
        update!(experiment, normalstats)
        bf10_exp = bayesfactor(experiment)

        @test bf10_model == bf10_exp
        @test isapprox(bf10_exp, 1/5.45, rtol=0.01)

        winner, _ = decide!(experiment)
        @test winner == "null"
    end

    @testset "Two-Sided alternative wins, thresh=5" begin
        normalstats = TwoNormalStatistics(
            NormalStatistics(meanx=28.8, sdx=13.5, n=133),
            NormalStatistics(meanx=33.6, sdx=14.3, n=867)
        )
        model = StudentTEffectSize(r=1.0)
        stats = StudentTStatistics(normalstats)
        bf10_model = bayesfactor(model, stats)

        thresh=5
        stoppingrule = TwoSidedBFThresh(thresh)
        experiment = ExperimentBF(model=model, rule=stoppingrule)
        update!(experiment, normalstats)
        bf10_exp = bayesfactor(experiment)

        @test bf10_model == bf10_exp
        @test isapprox(bf10_exp, 46.6077, rtol=0.01)

        winner, _ = decide!(experiment)
        @test winner == "alternative"
    end

    @testset "Two-Sided no winners, thresh=0.01" begin
        normalstats = TwoNormalStatistics(
            NormalStatistics(meanx=28.8, sdx=13.5, n=133),
            NormalStatistics(meanx=30.6, sdx=14.3, n=867)
        )
        model = StudentTEffectSize(r=1.0)
        stats = StudentTStatistics(normalstats)
        bf10_model = bayesfactor(model, stats)

        thresh=20
        stoppingrule = TwoSidedBFThresh(thresh)
        experiment = ExperimentBF(model=model, rule=stoppingrule)
        update!(experiment, normalstats)
        bf10_exp = bayesfactor(experiment)

        @test bf10_model == bf10_exp
        @test isapprox(bf10_exp, 1/5.45, rtol=0.01)

        winner, _ = decide!(experiment)
        @test winner === nothing
    end

    @testset "One-Sided no winner, thresh=5" begin
        normalstats = TwoNormalStatistics(
            NormalStatistics(meanx=28.8, sdx=13.5, n=133),
            NormalStatistics(meanx=30.6, sdx=14.3, n=867)
        )
        model = StudentTEffectSize(r=1.0)
        stats = StudentTStatistics(normalstats)
        bf10_model = bayesfactor(model, stats)

        thresh=5
        stoppingrule = OneSidedBFThresh(thresh)
        experiment = ExperimentBF(model=model, rule=stoppingrule)
        update!(experiment, normalstats)
        bf10_exp = bayesfactor(experiment)

        @test bf10_model == bf10_exp
        @test isapprox(bf10_exp, 1/5.45, rtol=0.01)

        winner, _ = decide!(experiment)
        @test winner === nothing
    end

    @testset "One-Sided alternative wins, thresh=5" begin
        normalstats = TwoNormalStatistics(
            NormalStatistics(meanx=28.8, sdx=13.5, n=133),
            NormalStatistics(meanx=33.6, sdx=14.3, n=867)
        )
        model = StudentTEffectSize(r=1.0)
        stats = StudentTStatistics(normalstats)
        bf10_model = bayesfactor(model, stats)

        thresh=5
        stoppingrule = OneSidedBFThresh(thresh)
        experiment = ExperimentBF(model=model, rule=stoppingrule)
        update!(experiment, normalstats)
        bf10_exp = bayesfactor(experiment)

        @test bf10_model == bf10_exp
        @test isapprox(bf10_exp, 46.6077, rtol=0.01)

        winner, _ = decide!(experiment)
        @test winner == "alternative"
    end
end


text:
total_cost(sim::Simulation) = fixed_cost(sim) + variable_cost(sim)

variable_cost(sim::Simulation) = sum(
    sim.acc.executed_dispatch[a][:cpeti] *
    sum(values(sim.acc.executed_dispatch[a][:flow]))
    for a in arcs(sim.acc.executed_dispatch)
)

fixed_cost(sim::Simulation) = sum(
    sim.acc.executed_dispatch[a][:cjarak] *
    sim.acc.executed_dispatch[a][:trip]
    for a in arcs(sim.acc.executed_dispatch)
)

function lost_sales(sim::Simulation)
    realization = filter(
        p -> p.periode <= sim.params.T, sim.libs.demand_realization
    )
    fulfilled = begin
        sum(abs.(sim.acc.demand_fulfillment.value)) / sum(abs.(realization.value))
    end
    return 1 - fulfilled
end

function network_structure(sim::Simulation)

end


text:
using MPI
using CLIMA
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using Printf
using LinearAlgebra
using Logging
using CLIMA.MPIStateArrays

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  const ArrayTypes = (CuArray, )
else
  const ArrayTypes = (Array, )
end

@inline function auxiliary_state_initialization!(aux, x, y, z)
  @inbounds begin
    r = hypot(x, y, z)
    aux[1] = r
    aux[2] = x / r
    aux[3] = y / r
    aux[4] = z / r
  end
end

using Test
function run(mpicomm, topl, ArrayType, N, DFloat)
  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = DFloat,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = 0,
                           flux! = (x...) -> (),
                           numerical_flux! = (x...) -> (),
                           numerical_boundary_flux! = (x...) -> (),
                           auxiliary_state_length = 4,
                           auxiliary_state_initialization! = (x...) ->
                           auxiliary_state_initialization!(x...))

  exact_aux = copy(spacedisc.auxstate)

  DGBalanceLawDiscretizations.grad_auxiliary_state!(spacedisc, 1, (2,3,4))

  euclidean_distance(exact_aux, spacedisc.auxstate)
end

let
  MPI.Initialized() || MPI.Init()
  Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())

  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
    device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end

  numelem = (5, 5, 1)

  polynomialorder = 4

  base_Nhorz = 4
  base_Nvert = 2
  Rinner = 1//2
  Router = 1

  expected_result = [2.6990557059226476e-04;
                     1.5939109758106752e-05;
                     1.0009810808262489e-06;
                     6.2636204797003758e-08]

  lvls = integration_testing ? length(expected_result) : 1

  for ArrayType in ArrayTypes
    for DFloat in (Float64,) #Float32)
      err = zeros(DFloat, lvls)
      for l = 1:lvls
        @info (ArrayType, DFloat, "sphere", l)
        Nhorz = 2^(l-1) * base_Nhorz
        Nvert = 2^(l-1) * base_Nvert
        Rrange = range(DFloat(Rinner); length=Nvert+1, stop=Router)
        topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
        err[l] = run(mpicomm, topl, ArrayType, polynomialorder, DFloat)
        @test expected_result[l] ≈ err[l]
      end
      if integration_testing
        @info begin
          msg = ""
          for l = 1:lvls-1
            rate = log2(err[l]) - log2(err[l+1])
            msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
          end
          msg
        end
      end
    end
  end
end

isinteractive() || MPI.Finalize()

nothing

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
using Base
using Base.Test
using CRF

v = Features(100)
w = Features(100)

# Check if adding features works like it should (you sould be able to access
# the value of global variables in the append! macro)
g = true
for i = 1:10, j = 1:10
    @append! v ((g) & ((i < 7) | (j > 4)))
end

g = false
for i = 1:10, j = 1:10
    @append! w ((g) & ((i < 7) | (j > 4)))
end

@test (v.i == 101) & (sum(v.x) == 84)
@test (w.i == 101) & (sum(w.x) == 0)

empty!(v)
empty!(w)
@test (v.i == 1) & (sum(v.x) == 0)
@test (w.i == 1) & (sum(w.x) == 0)


text:
function get_grids1()

    r = ( 1/β )
    w = 1/(Lₛ * r * 𝜋)
    t = (1 - Lₛ) * ϐ / ( Lₛ * w)
    loan_max =  (w * (1 - t) - ϐ) / ((r - 1) * w * (1 - t))

    loan_grid_points = Int(32)

    b = 0.25
    a = ( - log( loan_max ) + log( b )) / ( ( loan_grid_points - 1 ) * log( 1 - δᶠ))

    function GridGen(x::Int64,a_::Float64=a, b_::Float64=b)

        return (b_ * ( (1-δᶠ )^(-a_*x) ) )

    end

    loan_grid = Vector{Float64}([0.0])
    push!(loan_grid, [ GridGen(x) for x in 1:(loan_grid_points-1) ]...)

    depo_grid_points = Int(32)
    depo_grid = Vector{Float64}([0.0])
    push!(depo_grid, [ GridGen(x) for x in 1:(depo_grid_points-1) ]...)

    return(loan_grid, depo_grid)

end

function get_grids2(x::Int)

    r = ( 1/β )
    w = 1/(Lₛ * r * 𝜋)
    t = (1 - Lₛ) * ϐ / ( Lₛ * w)
    loan_max =  (w * (1 - t) - ϐ) / ((r - 1) ) #* w * (1 - t))
    depo_max = loan_max
    loan_grid_points = Int(x)

    # "distance between two neighbouring grid points. Used in utility function"
    𝔤 = loan_max / (loan_grid_points-1)
    loan_grid = 0:𝔤:loan_max

    # loan_grid = (0:(loan_grid_points - 1))*(loan_max/(loan_grid_points - 1))
    depo_grid_points = Int(x)

    𝔤 = depo_max / (depo_grid_points-1)
    depo_grid = 0:𝔤:loan_max

    return(loan_grid, depo_grid)

end

function get_grids3()

    r = ( 1/β )
    w = 1/(Lₛ * r * 𝜋)
    t = (1 - Lₛ) * ϐ / ( Lₛ * w)
    loan_max = 1.16 * (w * (1 - t) - ϐ) / ((r - 1) * w * (1 - t))
    depo_max = loan_max
    loan_min =  68.0
    depo_min = loan_min
    loan_grid_points = Int(32)

    loan_grid = [0,collect( range(loan_min, loan_max, length = (loan_grid_points) ))... ]

    # loan_grid = (0:(loan_grid_points - 1))*(loan_max/(loan_grid_points - 1))
    depo_grid_points = Int(32)

    depo_grid = [0,collect( range(loan_min, depo_max, length = (depo_grid_points) ))... ]

    return(loan_grid, depo_grid)

end

loan_grid, depo_grid = get_grids2(32)
#

## Set starting prices (vector Real):

## Define Grid. Run if have any questions (uncomment first):
# ?help GridDef
## make use of build in function from HANS_utilities:

grid_def = make_dimmensions( ["d", "l","λᵢ"], [ depo_grid , loan_grid, Int.([1, 2]) ])


text:
"""
    resample(f::Function, x::AbstractUncertainValue, n::Int, args...; kwargs...)

Draw an `n`-element sample from `x` according to its furnishing distribution, then
call `f(x_draw, args...; kwargs...)` on the length-`n` draws.
"""
function resample(f::Function, x::T, n::Int, args...; kwargs...) where T <: AbstractUncertainValue

    draw = resample(x, n)

    f(draw, args...; kwargs...)
end

"""
    resample(f::Function,
        x::AbstractUncertainValue, y::AbstractUncertainValue,
        n::Int, args...; kwargs...)

Draw an `n`-element sample from `x` according to its furnishing distribution, and
draw an `n`-element sample from `y` according to its furnishing distribution. Then,
call `f(x_draw, y_draw, args...; kwargs...)` on the length-`n` draws.
"""
function resample(f::Function,
        x::T1, y::T2,
        n::Int, args...; kwargs...) where {T1 <: AbstractUncertainValue, T2 <: AbstractUncertainValue}

    draw_x = float.(resample(x, n))
    draw_y = float.(resample(y, n))

    f(draw_x, draw_y, args...; kwargs...)
end

"""
    resample(f::Function, x::UVAL_COLLECTION_TYPES, n::Int, args...; kwargs...)

Resample the elements of `x` according to their furnishing uncertain values, yielding
a length-`l` realisation of `x` if `length(x) = l`. The elements of `x` are resampled
independently, assuming no sequential dependence between the elements. Then, call
`f(x, args...; kwargs...)` on the length-`l` sample with the given `args` and `kwargs`.
This process is repeated `n` times, yielding a length-`n` distribution of evaluations
of `f`.
"""
function resample(f::Function, x::T, n::Int, args...; kwargs...) where {
        T <: UVAL_COLLECTION_TYPES}

    draws = resample(x, n)
    [f(draw, args...; kwargs...) for draw in draws]
end

"""
    resample(f::Function, x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int, args...; kwargs...)

Resample the elements of `x` according to their furnishing uncertain values, yielding
a length-`l` realisation of `x` if `length(x) = l`. Then, do the same for `y`.
The elements of `x` and `y` are resampled independently, assuming no sequential dependence
between the elements of neither `x` nor `y`.

Then, call `f(x_draw, y_draw, args..., kwargs...)` on the length-`l` samples
`x_draw` and `y_draw`. This process is repeated `n` times, yielding a length-`n`
distribution of evaluations of `f`.
"""
function resample(f::Function, x::TX, y::TY, n::Int, args...; kwargs...) where {
        TX <: UVAL_COLLECTION_TYPES,
        TY <: UVAL_COLLECTION_TYPES}

    f_vals = zeros(Float64, n)

    for i = 1:n
        draw_x = float.(resample(x))
        draw_y = float.(resample(y))
        f_vals[i] = f(draw_x, draw_y, args...; kwargs...)
    end

    return f_vals
end

export resample


text:
export GradientBanditLearner

using Flux: softmax, onehot

"""
    GradientBanditLearner(;approximator::A, optimizer::O, baseline::B)
"""
Base.@kwdef mutable struct GradientBanditLearner{A,O,B} <: AbstractLearner
    approximator::A
    optimizer::O
    baseline::B
end

(learner::GradientBanditLearner)(s::Int) = s |> learner.approximator |> softmax
(learner::GradientBanditLearner)(obs) = learner(get_state(obs))

RLBase.update!(learner::GradientBanditLearner, ::Nothing) = nothing

function RLBase.update!(learner::GradientBanditLearner, experience::Tuple)
    s, a, r = experience
    probs = s |> learner.approximator |> softmax
    r̄ = learner.baseline isa Number ? learner.baseline : learner.baseline(r)
    errors = (r - r̄) .* (onehot(a, 1:length(probs)) .- probs)
    update!(learner.approximator, s => apply!(learner.optimizer, s, errors))
end

function extract_experience(t::AbstractTrajectory, ::GradientBanditLearner)
    if length(t) > 0
        (
            get_trace(t, :state)[end],
            get_trace(t, :action)[end],
            get_trace(t, :reward)[end],
        )
    else
        nothing
    end
end


text:
function trace_interpolation(d::Array{Tv}, order::Ti) where {Tv<:AbstractFloat, Ti<:Int64}

    Ndims = ndims(d)
    Ndims == 1 ? nx = 1 :  nx = prod(size(d)[2:Ndims])
    dims  = size(d)
    N     = dims[1]
    Npad  = order*N
    nf    = 2*Npad
    nw    = convert(Int,floor(nf/2)) + 1
    di    = zeros(eltype(d),nf,dims[2:end]...)

    for k = 1:nx

        dd = zeros(eltype(d),nf)
        dd[1:order:Npad] = d[1:1:end,k]
        D = fft(dd)
        nyq = convert(Int,floor(nw/order)) + 1
        for iw = nyq : nw
            D[iw] *= 0
        end

        # symmetries
        for iw=nw+1:nf
            D[iw] = conj(D[nf-iw+2])
        end
        di[:,k] = real(ifft(D,1))
    end

    return  reshape(di[1:Npad,:]*order,Npad,dims[2:end]...)
end


text:
"""
## Convex hull

Calculation f convex hull of a sequence of points represented through an Float64
Array of N rows and 2 columns
"""
function convexhull{T <: Real}(points::Array{T,2})
    n = size(points,1)
    # Sort the points by x-coordinate, resulting in a sequence p 1 , . . . , p n
    p = sortrows(points)
    # Put the points p_1 and p_2 in a list L_upper , with p_1 as the first point
    upper = [p[1,:], p[2,:]]

    for i in 3:n
        while size(upper,1) >= 2 && direction(upper[end-1], upper[end], p[i,:]) <= 0
            pop!(upper)
        end

        push!(upper,p[i,:])
    end

    lower = [p[end,:], p[end-1,:]]

    for i in n-2:-1:1
        while size(lower,1) >= 2 && direction(lower[end-1], lower[end], p[i,:]) <= 0
            pop!(lower)
        end

        push!(lower,p[i,:])
    end

    append!(upper, lower[2:end-1])
end

function convexhull{T <: Real}(points::Array{Point2D{T}, 1})
    n = size(points,1)
    # Sort the points by x-coordinate, resulting in a sequence p 1 , . . . , p n
    p = sort(points)
    # Put the points p_1 and p_2 in a list L_upper , with p_1 as the first point
    upper = [p[1], p[2]]

    for i in 3:n
        while size(upper,1) >= 2 && direction(upper[end-1], upper[end], p[i]) <= 0
            pop!(upper)
        end

        push!(upper,p[i])
    end

    lower = [p[end], p[end-1]]

    for i in n-2:-1:1
        while size(lower,1) >= 2 && direction(lower[end-1], lower[end], p[i]) <= 0
            pop!(lower)
        end

        push!(lower,p[i])
    end

    append!(upper, lower[2:end-1])
end

max_stars_repo_name: helgith/GeometricStructures

max_issues_repo_name: helgith/GeometricStructures

max_forks_repo_name: helgith/GeometricStructures


text:
using SimradRaw
using SimradEK60TestData
filename = EK60_SAMPLE # or some other RAW file
datagrams = SimradRaw.load(filename)


text:
"""
    AndersonMooreAlg(h, qcols, neq)

Solve a linear perfect foresight model using the julia eig
function to find the invariant subspace associated with the big
roots.  This procedure will fail if the companion matrix is
defective and does not have a linearly independent set of
eigenvectors associated with the big roots.

  Input arguments:

    h         Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
    neq       Number of equations.
    nlag      Number of lags.
    nlead     Number of leads.
    condn     Zero tolerance used as a condition number test
              by numericShift and reducedForm.
    upper     Inclusive upper bound for the modulus of roots
              allowed in the reduced form.

  Output arguments:

    b         Reduced form coefficient matrix (neq,neq*nlag).
    rts       Roots returned by eig.
    ia        Dimension of companion matrix (number of non-trivial
              elements in rts).
    nexact    Number of exact shiftRights.
    nnumeric  Number of numeric shiftRights.
    lgroots   Number of roots greater in modulus than upper.
    AMAcode   Return code: see function AMAerr.
"""
function AndersonMooreAlg(hh::Array{Float64,2}, neq::Int64, nlag::Int64, nlead::Int64, anEpsi::Float64, upper::Float64)

    if(nlag < 1 || nlead < 1)
        error("AMA_eig: model must have at least one lag and one lead.")
    end

    # Initialization.
    nexact   = 0
    nnumeric = 0
    lgroots  = 0
    iq       = 0
    AMAcode  = 0
    bb       = 0
    qrows    = neq * nlead
    qcols    = neq * (nlag + nlead)
    bcols    = neq * nlag
    qq       = zeros(qrows, qcols)
    rts      = zeros(qcols, 1)

    # Compute the auxiliary initial conditions and store them in q.

    (hh, qq, iq, nexact) = exactShift!(hh, qq, iq, qrows, qcols, neq)
    if (iq > qrows)
        AMAcode = 61
        return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
    end

    (hh, qq, iq, nnumeric) = numericShift!(hh, qq, iq, qrows, qcols, neq, anEpsi)
    if (iq > qrows)
        AMAcode = 62
        return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
    end

    #  Build the companion matrix.  Compute the stability conditions, and
    #  combine them with the auxiliary initial conditions in q.

    (aa, ia, js) = buildA!(hh, qcols, neq)

    if (ia != 0)
        for element in aa
            if isnan(element) || isinf(element)
                display("A is NAN or INF")
                AMAcode = 63
                return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
            end
        end
        (ww, rts, lgroots) = eigenSys!(aa, upper, min(size(js, 1), qrows - iq + 1))

        qq = augmentQ!(qq, ww, js, iq, qrows)
    end

    test = nexact + nnumeric + lgroots
    if (test > qrows)
        AMAcode = 3
    elseif (test < qrows)
        AMAcode = 4
    end

    # If the right-hand block of q is invertible, compute the reduced form.

    if(AMAcode == 0)
        (nonsing,bb) = reducedForm(qq, qrows, qcols, bcols, neq, anEpsi)
        if ( nonsing && AMAcode==0)
            AMAcode =  1
        elseif (!nonsing && AMAcode==0)
            AMAcode =  5
        elseif (!nonsing && AMAcode==3)
            AMAcode = 35
        elseif (!nonsing && AMAcode==4)
            AMAcode = 45
        end
    end

    return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
 # (bbJulia,rtsJulia,iaJulia,nexJulia,nnumJulia,lgrtsJulia,AMAcodeJulia) = AMAalg(hh,neq,nlag,nlead,anEpsi,1+anEpsi)
end # AndersonMooreAlg


text:
using Random

Random.rand(rng::AbstractRNG, ::Random.SamplerType{Infinite}) = Infinite(rand(rng, Bool))


text:
# Reconstruct a String from its Genome Path

function calc_string_spelled_by_genome_path(kmers)
    text = collect(kmers[1])
    for kmer in kmers[2:end]
        push!(text, kmer[end])
    end
    join(text)
end

function main()
    data = readlines()
    result = calc_string_spelled_by_genome_path(data)
    println(result)
end

main()


text:
using LinearAlgebra
using SparseArrays
using Plots

"This code solves the 2D Laplace's equation with nonzero Dirichlet and Neumann
boundary conditions using a finite difference method."

m = 100 # number of points
x = LinRange(-1,1,m+2)
y = LinRange(-1,1,m+2)
xint = x[2:end-1]
yint = y[2:end-1]
h = x[2]-x[1]

# ordering of unknowns
# u14 u24 u34 u44
# u13 u23 u33 u43
# u12 u22 u32 u42
# u11 u21 u31 u41

id(i,j,m) = i + (j-1)*m # return global index from i,j

function FD_matrix_2D(h,m)
    A = spzeros(m*m,m*m)
    for i = 1:m, j = 1:m # loop thru x,y indices
        A[id(i,j,m),id(i,j,m)] = 4/h^2
        if i > 1 # avoids leftmost line of nodes
            A[id(i,j,m),id(i-1,j,m)] = -1/h^2 # x-derivative
        end
        if i < m # avoids rightmost line of nodes
            A[id(i,j,m),id(i+1,j,m)] = -1/h^2
        end
        if j > 1
            A[id(i,j,m),id(i,j-1,m)] = -1/h^2 # y-derivative
        end
        if j < m
            A[id(i,j,m),id(i,j+1,m)] = -1/h^2
        end
    end
    return A
end

A = FD_matrix_2D(h,m)

f(x,y) = exp(-10*((x-.25)^2 + (y-.75)^2))
f(x,y) = 0.0
b = vec(f.(xint,yint'))

# impose non-zero Dirichlet BCs
u_left(y) = 0.0
u_right(y) = 1.0 # 0 if y in [-1,0], 1 if y in [0,1]
u_bottom(x) = 0.0
u_top(x) = -1.0

for i = 1:m, j = 1:m
    if i==1
        b[id(i,j,m)] += u_left(yint[j])/h^2
    end
    if i==m
        b[id(i,j,m)] += u_right(yint[j])/h^2
    end
    if j==1
        b[id(i,j,m)] += u_bottom(xint[i])/h^2
    end
    if j==m
        b[id(i,j,m)] += u_top(xint[i])/h^2
    end
end

# solve Au = f(xi,yj)
u = A\b
contourf(xint,yint,u,clims=(-1,1)) #vec(reshape(u,m,m)'))


text:
# System module.
module System

import PyCall: pyimport, pycall

export LIB_SKL_AVAILABLE,
       LIB_CRT_AVAILABLE

function check_py_dep(package::String)
  is_available = true
  try
    pyimport(package)
  catch
    is_available = false
  end
  return is_available
end

function check_r_dep(package::String)
  is_available = true
  try
    pycall(pyimport("rpy2.robjects.packages")["importr"], Any, package)
  catch
    is_available = false
  end
  return is_available
end

# Check system for python dependencies.
LIB_SKL_AVAILABLE = check_py_dep("sklearn")
LIB_CRT_AVAILABLE = check_r_dep("caret")

end # module


text:
using .StencilBuilder

const ifactor = 1_000_000

mutable struct Timeloop
    # Float numbers.
    start_time::Float64
    end_time::Float64
    save_time::Float64
    check_time::Float64
    dt::Float64
    time::Float64

    # Integer values.
    istart_time::Int64
    iend_time::Int64
    isave_time::Int64
    icheck_time::Int64
    idt::Int64
    itime::Int64

    rkstep::Int64
end

function Timeloop(d::Dict)
    start_time = d["start_time"]
    end_time = d["end_time"]
    save_time = d["save_time"]
    check_time = d["check_time"]
    dt = d["dt"]
    time = start_time

    istart_time = round(Int64, start_time * ifactor)
    iend_time = round(Int64, end_time * ifactor)
    isave_time = round(Int64, save_time * ifactor)
    icheck_time = round(Int64, check_time * ifactor)
    idt = round(Int64, dt * ifactor)
    itime = istart_time

    rkstep = 1

    Timeloop(
        start_time, end_time, save_time, check_time, dt, time,
        istart_time, iend_time, isave_time, icheck_time, idt, itime,
        rkstep)
end

function integrate_time_kernel!(
    a, at,
    rkstep, dt,
    is, ie, js, je, ks, ke)

    c_a = [0, -5//9, -153//128];
    c_b = [1//3, 15//16, 8//15];

    rkstep_next = (rkstep + 1) == 4 ? 1 : rkstep + 1

    c_b_step = c_b[rkstep] * dt
    c_a_step = c_a[rkstep_next]

    @fast3d begin
        @fd (a, at) a += c_b_step*at
        @fd (a, at) at *= c_a_step
    end
end

function integrate_time!(
    f::Fields, g::Grid, t::Timeloop)

    # Make sure the dt is in Float32 if the array is.
    dt = convert(typeof(f.u[1]), t.dt)

    integrate_time_kernel!(
        f.u, f.u_tend,
        t.rkstep, dt,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)

    integrate_time_kernel!(
        f.v, f.v_tend,
        t.rkstep, dt,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)

    integrate_time_kernel!(
        f.w, f.w_tend,
        t.rkstep, dt,
        g.is, g.ie, g.js, g.je, g.ks, g.keh)

    integrate_time_kernel!(
        f.s, f.s_tend,
        t.rkstep, dt,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)
end

function step_time!(t::Timeloop)
    t.rkstep += 1

    if t.rkstep > 3
        t.rkstep = 1
        t.itime += t.idt
        t.time = convert(Float64, t.itime) / ifactor
    end
end

function get_sub_dt(t::Timeloop)
    c_b = [1//3, 15//16, 8//15];
    return c_b[t.rkstep]*t.dt;
end

function is_in_progress(t::Timeloop)
    return t.itime < t.iend_time
end

function is_save_time(t::Timeloop)
    return t.itime % t.isave_time == 0 && t.rkstep == 1 && t.itime != t.istart_time
end

function is_check_time(t::Timeloop)
    return t.itime % t.icheck_time == 0
end


text:
module StaticTools

    if !(v"1.8.0-DEV" < VERSION < v"1.9.0-DEV")
        @warn "StaticTools + StaticCompiler integration is best supported on Julia 1.8.\n Consider switching to a different Julia version"
    end

    # External dependencies
    using ManualMemory: MemoryBuffer, load, store!

    # Declare some types we'll use later
    struct FILE end # Plain struct to denote and dispatch on file pointers

    # Manual memory allocation
    include("mallocarray.jl")
    include("staticrng.jl")

    # String handling
    include("abstractstaticstring.jl")  # Shared string infrastructure
    include("unescape.jl")      # You don't want to know
    include("staticstring.jl")  # StaticCompiler-safe stack-allocated strings
    include("mallocstring.jl")  # StaticCompiler-safe heap-allocated strings

    # Union of things that don't need GC.@protect
    const AbstractMallocdMemory = Union{MallocString, MallocArray}

    # Here there be `llvmcall`s
    include("llvmio.jl")        # Best way to print things? LLVM IR obviously!
    include("llvmlibc.jl")      # strtod, strtol, parse, etc...

    # higher-level printing, parsing, etc.
    include("printformats.jl")
    include("parsedlm.jl")

    # Types
    export StaticString, MallocString, StringView, AbstractStaticString
    export MallocArray, MallocMatrix, MallocVector
    # Macros
    export @c_str, @m_str, @mm_str
    # Functions
    export malloc, calloc, free, memset!, memcpy!, memcmp                       # Memory management
    export stdinp, stdoutp, stderrp                                             # File pointers
    export fopen, fclose, frewind, fseek, SEEK_SET, SEEK_CUR, SEEK_END          # File open, close, seek
    export usleep                                                               # Other libc utility functions
    export newline, putchar, getchar, getc, puts, gets!, printf                 # Char & String IO
    export unsafe_mallocstring, strlen                                          # String management
    export parsedlm                                                             # File parsing
    export static_rng, xoshiro256✴✴, Xoshiro256✴✴, splitmix64, SplitMix64       # Random number generation
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b141a11c-6857-11eb-00c9-0d2d762cd167
begin
	using Test
	using BenchmarkTools
	using PlutoUI
end

# ╔═╡ 9b574d24-6855-11eb-1af8-d1c2ccf2f3b0
md"""
### Alphametcis aka CryptArithmetics

[Alphametics](https://en.wikipedia.org/wiki/Verbal_arithmetic) is a type of cryptarithm (or  cryptarithmetic) in which a set of words is written down in the form of a long addition sum or some other mathematical problem. The objective is to replace the letters of the alphabet with decimal digits to make a valid arithmetic sum.

For this kata, your objective is to write a function that accepts an alphametic equation in the form of a single-line string and returns a valid arithmetic equation in the form of a single-line string.

"""

# ╔═╡ e420b680-6855-11eb-2464-8b09ef6b45ae
md"""
**Test Examples**

INPUT: `"SEND + MORE = MONEY"`
SOLUT: `"9567 + 1085 = 10652"`
```
  S E N D
+ M O R E
---------
M O N E Y
```

INPUT: `"ELEVEN + NINE + FIVE + FIVE = THIRTY"`
SOLUT: `"797275 + 5057 + 4027 + 4027 = 810386"`
```
  E L E V E N
+     N I N E
+     F I V E
+     F I V E
-------------
  T H I R T Y
```

Some puzzles may have multiple valid solutions; your function only needs to return one:
BIG + CAT = LION
403 + 679 = 1082
326 + 954 = 1280
304 + 758 = 1062
...etc.

Technical Details
- All alphabetic letters in the input will be uppercase
- Each unique letter may only be assigned to one unique digit (bijective mapping)
- As a corollary to the above, there will be a maximum of 10 unique letters in any given test
- No leading zeroes
- The equations will only deal with addition with multiple summands on the left side
(summands are LHS, ex: SEND, MORE ...) and one term on the right side
- The number of summands will range between 2 and 7, inclusive
- The length of each summand will range from 2 to 8 characters, inclusive
- All test cases will be valid and will have one or more possible solutions
- Optimize your code -- a naive, brute-force algorithm may time out before the first test completes
"""

# ╔═╡ 93fbf81c-6856-11eb-0a43-bfaa709a504a
md"""
**Resolution**

```
  S E N D
+ M O R E
---------
M O N E Y
```

*  $M is 1$, obviously (it cannot be 0)
*  $D + E = Y + 10.c_1$
*  $c_1 + N + R = E + 10.c_2$
*  $c_2 + E + O = N + 10.c_3$
*  $c_3 + S + M = O + 10.c_4$
*  $c_4 = M  \equiv 1 = M$

Hence $\color{blue}{M = 1}$

System can be rewritten as follows:
* (i) $\space D + E = Y + 10.c_1$
* (ii) $\space c_1 + N + R = E + 10.c_2$
* (iii) $\space c_2 + E + O = N + 10.c_3$
* (iv) $\space c_3 + S + 1 = O + 10 \equiv c_3 + S = O + 9$

From (iv):
$\space c_3 + S = O + 9 \Rightarrow (c_3=0 \land S=9, O=0) \lor (c_3=1 \Rightarrow S = O + 8 \Rightarrow S \ge 8, O \lt 2) => O=1 \Rightarrow \varnothing$ ($M=1$ already).

Hence $\color{blue}{O = 0}$

From (iii) $r_2 + E + 0 = N + 10.r_3 \Rightarrow r_2 + E = N + 10.r_3$
$r_2 + E = N + 10.r_3$, if $r_3 = 1 \Rightarrow r_2 + E = N + 10$, 2 cases:
*  $r_2 = 0, \space 0 + E = N + 10 \equiv E = N + 10$, this is impossible
*  $r_2 = 1, \space 1 + E = N + 10 \equiv E = N + 9 => N=0 \land E=9$, this is impossible, because $O=0$

Hence $r_3 = 0 \land \color{blue}{S=9}$

From (iii): $r_2 + E = N$, 2 cases:
*  $r_2 = 0, \space E = N \space \Rightarrow impossible$
*  $r_2 = 1, \space 1 + E = N$

Hence $r_2 = 1$

From (ii): $r_1 + N + R = E + 10 \equiv r_1 + N + R = N - 1 + 10 \equiv r_1 + R = 9$, given $S = 9$:
$r_1 = 1 \land \color{blue}{R = 8}$

From (i), $D + E = Y + 10$ and $1 + E = N$:
* we have $Y \ge 1$
* we know maximum $D + E = 13$ (as 9, 8 are already taken, so D, E would be (6, 7) or (7, 6).
* But $E = 7 \Rightarrow N = 8$, impossible
* And $E = 6 \Rightarrow N = 7 \land D = 7$, impossible

Thus:
* E is 5, 4, ~3, 2, 1~ and
* N is 6, 5, ~4, 3, 2~ and
* D 7, 6 and
* Y = (7 + 5 - 10 =) 2, ~(6 + 4 -10) 0~

Hence $\color{blue}{Y = 2}$ (only possibility) $\Rightarrow \color{blue}{E =5} \land \color{blue}{N=6}$

"""

# ╔═╡ 9c4e7976-6857-11eb-3197-7da4e0930506
md"""
Some more examples:
* `  D A Y S  +     T O O =   S H O R T`
* `    T W O  +   D A Y S =     M O R E`
* `  W A I T  +     A L L =   G I F T S`
* `  N I N E  +   F I V E =   W I V E S`
* `  E V E R  + S I N C E = D A R W I N`
* `    U S A  +   U S S R =   P E A C E`
* `P O I N T  +   Z E R O = E N E R G Y`
* `      G O  +       T O =       O U T`

* ` E A T  +  E A T  +  E A T = B E E T`
"""

# ╔═╡ bce4f3cc-6857-11eb-19f4-dd998acaeafe
begin
	const ALL_DIGITS = [d for d in 0:9]

	_input = "SEND + MORE = MONEY"
	_s, _result = map(strip, split(_input, '='))
end

# ╔═╡ bee7b3f8-6857-11eb-2ccd-e7202568155d
_summands = map(strip, split(_s, "+"))

# ╔═╡ bfcee1b0-6857-11eb-36b6-e1bb043a08e8
function words_result(str::String)
    s, res = map(strip, split(str, '='))
    summands = map(strip, split(s, "+"))

    (result=split(res, ""), summands=map(s -> split(s, ""), summands))
end

# ╔═╡ c14aade4-6857-11eb-21ea-51b53292b1af
with_terminal() do
	println(words_result("SEND + MORE = MONEY"))
end

# ╔═╡ c1f15cf0-6857-11eb-2fb8-d95ee8623319
function wr_tuple_rev(wr_tuple)
    """
    ```jldoctest
    julia> wr_tuple = (result = ["M", "O", "N", "E", "Y"], summands = [["S", "E", "N", "D"], ["M", "O", "R", "E"]])
    julia> wr_tuple_rev(wr_tuple)
    (result = ["Y", "E", "N", "O", "M"], summands = [["D", "N", "E", "S"], ["E", "R", "O", "M"])
    ```
    """
    res = reverse(wr_tuple.result)
    sums = map(reverse, wr_tuple.summands)

	(result=res, summands=sums)
end

# ╔═╡ c20ce150-6857-11eb-0c50-bdb013339921
with_terminal() do
	wr_tuple = (result = ["M", "O", "N", "E", "Y"],
		summands = [["S", "E", "N", "D"], ["M", "O", "R", "E"]])

	println(wr_tuple_rev(wr_tuple))
end

# ╔═╡ 48cb2262-6858-11eb-3dc4-fbfb36ff2ab3
begin
	wr_tuple1 = (result = ["M", "O", "N", "E", "Y"],
		summands = [["S", "E", "N", "D"], ["M", "O", "R", "E"]])
	@test wr_tuple_rev(wr_tuple1) == (result = ["Y", "E", "N", "O", "M"],
		summands = [["D", "N", "E", "S"], ["E", "R", "O", "M"]])

	wr_tuple2 = (result = ["T", "H", "I", "R", "T", "Y"],
		summands = [["E", "L", "E", "V", "E", "N"], ["N", "I", "N", "E"], ["F", "I", "V", "E"], ["F", "I", "V", "E"]])

	@test wr_tuple_rev(wr_tuple2) == (result = ["Y", "T", "R", "I", "H", "T"],
		summands = [["N", "E", "V", "E", "L", "E"], ["E", "N", "I", "N"], ["E", "V", "I", "F"], ["E", "V", "I", "F"]])

end

# ╔═╡ 8d0d5940-6858-11eb-1155-dbe742a87729
function all_letters(str::String)
    split(replace(str, ['+', ' ', '='] => ""),
		"") |>
		sort |>
		unique
end

# ╔═╡ e2b38092-6858-11eb-1305-115e60d9ce26
with_terminal() do
	println(all_letters("SEND + MORE = MONEY"))
end

# ╔═╡ 123426e4-6859-11eb-1857-150829245eca
begin
	@test all_letters("SEND + MORE = MONEY") == [
		"D", "E", "M", "N", "O", "R", "S", "Y"
	]

	@test all_letters("ELEVEN + NINE + FIVE + FIVE = THIRTY") == [
		"E", "F", "H", "I", "L", "N", "R", "T", "V", "Y"
	]
end

# ╔═╡ 3882e254-6859-11eb-19d4-bfb9921902f2
mutable struct AExpr
    lhs::String
    rhs::String
end

# ╔═╡ 3f8cfd78-6859-11eb-07c2-8f02f60fc0a3
Base.show(io::IO, aexpr::AExpr) = print(io, "lhs: $(aexpr.lhs) = rhs: $(aexpr.rhs)")

# ╔═╡ 45a9e11c-6859-11eb-3266-cb2f4081e7de
function gen_constraints(wr_tuple_rev)::Tuple{Array{AExpr, 1}, Array{String, 1}}
    len_res = length(wr_tuple_rev.result)
    constraints = []
    ix, aexpr = 1, AExpr("", "")
    symbols = Vector{String}([])

    while true
        ## lhs
        for ary in wr_tuple_rev.summands
            if ix ≤ length(ary)
                ary[ix] ∈ symbols || push!(symbols, ary[ix])
                aexpr.lhs = isempty(aexpr.lhs) ? "$(ary[ix])" : "$(aexpr.lhs) + $(ary[ix])"
            end
        end

        ## rhs
        wr_tuple_rev.result[ix] ∈ symbols || push!(symbols, wr_tuple_rev.result[ix])
        aexpr.rhs = ix == length(wr_tuple_rev.result) ? "$(wr_tuple_rev.result[ix])" :
        	"$(wr_tuple_rev.result[ix]) + 10 * c_$(ix)"
        push!(constraints, aexpr)

        ## next
        ix += 1
        ix > len_res && break

        ## next carry
        aexpr = AExpr("c_$(ix - 1)", "")
        push!(symbols, "c_$(ix - 1)")
    end

    (constraints, symbols)
end

# ╔═╡ 10d6c6a2-685a-11eb-269d-3f92f1a94fe1
begin
	wr_tuple_rev1 = wr_tuple_rev(wr_tuple1)

	constraint_ary, symbol_ary = gen_constraints(wr_tuple_rev1)
	with_terminal() do
		println("constraints: $(constraint_ary)")
		println("symbols: $(symbol_ary)")
	end
end

# ╔═╡ 8b681600-685a-11eb-0bfd-bbaf19e1a22b
md"""
Remarks:
-  $\forall i, c_i ∈ [0, 1]$
-  $S \neq 0 \land M \neq 0$
"""

# ╔═╡ b8855382-685a-11eb-2bd4-7f7d3557d836
begin
	const ETYPES = (:carry, :var)  # Symbol Expression type either :carry or :var

	with_terminal() do
		println(typeof(ETYPES))
	end
end

# ╔═╡ e2e282da-685a-11eb-32b2-116fce22d9d3
mutable struct Elem
    symbol::String
    value::Union{Nothing, Int}
    pvalues::Set{Int}
    etype::Symbol
end

# ╔═╡ e71d56b8-685a-11eb-367c-c37ad541bdf4
begin
	hsymbols = Dict{String, Elem}()

	for k ∈ symbol_ary
    	cond = match(r"\Ac_\d+", k) != nothing
    	hsymbols[k] = Elem(
        	k,
        	nothing,
        	cond ? Set([0, 1]) : Set(collect(0:9)),
        	cond ? ETYPES[1] : ETYPES[2] # ETYPES][1] == :carry, ETYPES][2] == :var
    	)
end

	with_terminal() do
		println(hsymbols)
	end
end

# ╔═╡ 0f3277dc-685b-11eb-2bc5-d7cadea01388
begin
const NT{T} = NamedTuple{(:lhs, :rhs), T} where T <: Tuple

function gen_constraints_as_expr(ary::Vector{AExpr})::Vector{NT}
	map(aexpr -> (lhs=Meta.parse(aexpr.lhs), rhs=Meta.parse(aexpr.rhs)),
			ary)
end

end

# ╔═╡ 9fb53268-685b-11eb-18cc-5d103022febf
const AEXPR_LST = gen_constraints_as_expr(constraint_ary)

# ╔═╡ 0c9579e2-685c-11eb-0985-b323f500064c
md"""
#### Evaluation
"""

# ╔═╡ 1816ffd4-685c-11eb-218c-3985c1e8fd40
function eval_aexpr(aexpr_lst::Array{NamedTuple{(:lhs, :rhs),T} where T<:Tuple,1}, d_env)
    # set the env
    [@eval $(Symbol("$(k)")) = ($v) for (k, v) in d_env]

    for aexpr in aexpr_lst
      try
      	lhs = eval(aexpr.lhs)
        rhs = eval(aexpr.rhs)

        if (lhs != rhs)
          println("eval: $(aexpr.lhs)[==$(lhs)] != $(aexpr.rhs)[==$(rhs)] - false / stop")
          return false
        else
          println("eval: $(aexpr.lhs)[==$(lhs)] == $(aexpr.rhs)[==$(rhs)] - true / continue")
        end

      catch ex
      	println("[!] lhs: $(aexpr.lhs) / rhs: $(aexpr.rhs) - not eval-uable")
        continue
      end
    end

    true
end

# ╔═╡ 686aed92-685c-11eb-295b-81d603e5eedc
with_terminal() do
	eval_aexpr(AEXPR_LST,
		merge(hsymbols, Dict("D" => 9, "E" => 2, "Y" => 3, "c_1" => 0)))
end

# ╔═╡ 9a573124-685c-11eb-107a-05702bbfc2f2
with_terminal() do
	eval_aexpr(AEXPR_LST,
		merge(hsymbols, Dict("c_4" => 1, "M" => 1)))
end

# ╔═╡ a0855ff0-685c-11eb-3b70-b7c6d257b1bf
with_terminal() do
	eval_aexpr(AEXPR_LST,
		merge(hsymbols, Dict("c_4" => 1, "M" => 2, "D" => 3, "E" => 1, "Y" => 5, "c_1" => 0)))
end

# ╔═╡ bd6bf02a-685c-11eb-1b01-d9dfda6e8f91
with_terminal() do
	eval_aexpr(AEXPR_LST,
		merge(hsymbols, Dict("c_4" => 1, "M" => 2, "D" => 4, "E" => 1, "Y" => 5, "c_1" => 0)))
end

# ╔═╡ db3fe086-685c-11eb-0a28-e7529365d7b8
with_terminal() do
	eval_aexpr(AEXPR_LST,
		merge(hsymbols, Dict("c_4" => 1, "M" => 1, "D" => 5, "E" => 2, "Y" => 7, "c_1" => 0)))
end

# ╔═╡ f288a028-685c-11eb-07f2-3de80f49b384
md"""
#### System

Possible Strategy:
* Start with equation with less constraints (in term of length and number of constraints
* LHS (1 symb) = RHS (1 symb) and LHS $\neq$ RHS $\Rightarrow$ solution = LHS $\cap$ RHS

Constraints:
* Leftmost symbols cannot be 0
* All Symbols have different values

"""

# ╔═╡ ff22edf2-685c-11eb-226e-c9470f8902e9
with_terminal() do
	println(AEXPR_LST)
end

# ╔═╡ 2304a31e-685d-11eb-2c68-037fefdb684f
with_terminal() do
	t = AEXPR_LST[1]

	println(cat(t.lhs.args[2:end], t.rhs.args[2:end]; dims=1))
end

# ╔═╡ 38b95e48-685d-11eb-3c4a-8fd1e28d00a6
with_terminal() do
	t = AEXPR_LST[1]

	println(vcat(t.lhs.args[2:end], t.rhs.args[2:end]))
end

# ╔═╡ 659a12a2-685d-11eb-1c32-6b53b431ac9a
function find_eqn_with_less_contraints(aexpr_lst)
    """
    Find equation in given aexpr_lst, which has the less element(s)
    returns first equation if all equations have same number of elements (or members)
    """
    ary, nmin = fill(0, 3), 10_000_000

    for (ix, t) in enumerate(aexpr_lst)
        _ary = fill(ix, 3)
        for (jx, ths) in enumerate(t)
            _ary[jx + 1] = if typeof(ths) == Expr
				length(ths.args[2:end])
			elseif typeof(ths) == Symbol
				1
			end
        end
        n = _ary[2] + _ary[end]
        if n < nmin
            nmin = n
            for i in 1:3
                ary[i] = _ary[i]
            end
        end
    end
    ary
end

# ╔═╡ b3775770-685d-11eb-3090-dbbc1689f1d0
with_terminal() do
	ary = find_eqn_with_less_contraints(AEXPR_LST)
	println(ary)
end

# ╔═╡ e5ae84de-685d-11eb-2377-918cd4b53af4
function non_leading_0!(hsymbols, wr_tuple)
    sym = wr_tuple.result[1]
    hsymbols[sym].pvalues = Set(collect(1:9))

    for s in wr_tuple.summands
        hsymbols[s[1]].pvalues = Set(collect(1:9))
    end

    hsymbols
end

# ╔═╡ f2fb7c8c-685d-11eb-0def-01b7643c6547
with_terminal() do
	wr_tuple = (result = ["M", "O", "N", "E", "Y"],
		summands = [["S", "E", "N", "D"], ["M", "O", "R", "E"]])

	println(non_leading_0!(hsymbols, wr_tuple))
end

# ╔═╡ 4ced4446-685e-11eb-3431-c3de701b63a6
function lhs_rhs_singleton!(hsymb, aexpr_lst)
    ix, n_l, n_r = find_eqn_with_less_contraints(aexpr_lst)

    if n_l == n_r && n_r == 1
		## we have 1 equation with singleton element on both side, simplify!
        lhs_sym = aexpr_lst[ix].lhs
        rhs_sym = aexpr_lst[ix].rhs
        set_ = intersect(
			hsymb[String(lhs_sym)].pvalues,
			hsymb[String(rhs_sym)].pvalues
		)
        hsym[String(lhs_sym)].pvalues = set_
        hsym[String(rhs_sym)].pvalues = set_

        if length(set_) == 1
            hsymb[String(lhs_sym)].value = first(set_)
            hsymb[String(rhs_sym)].value = first(set_)

			## this value needs to be remove from all pvalues of all
			## remaining variables
            for sym ∈ keys(hsymb)
                hsymb[sym].etype == :carry && continue # ignore carries
                hsymb[sym].pvalues = setdiff(hsymb[sym].pvalues, set)
            end
        end
    end
end

# ╔═╡ dcc9489e-685e-11eb-315e-e34875afae1e
with_terminal() do
	println(eval_aexpr(AEXPR_LST,
		Dict(k => hsymbols[k].value for k in keys(hsymbols))
	))
end

# ╔═╡ fb37bad6-685e-11eb-0314-73160a66b806
md"""
*TBC...*
"""

# ╔═╡ Cell order:
# ╟─9b574d24-6855-11eb-1af8-d1c2ccf2f3b0
# ╟─e420b680-6855-11eb-2464-8b09ef6b45ae
# ╠═93fbf81c-6856-11eb-0a43-bfaa709a504a
# ╟─9c4e7976-6857-11eb-3197-7da4e0930506
# ╠═b141a11c-6857-11eb-00c9-0d2d762cd167
# ╠═bce4f3cc-6857-11eb-19f4-dd998acaeafe
# ╠═bee7b3f8-6857-11eb-2ccd-e7202568155d
# ╠═bfcee1b0-6857-11eb-36b6-e1bb043a08e8
# ╠═c14aade4-6857-11eb-21ea-51b53292b1af
# ╠═c1f15cf0-6857-11eb-2fb8-d95ee8623319
# ╠═c20ce150-6857-11eb-0c50-bdb013339921
# ╠═48cb2262-6858-11eb-3dc4-fbfb36ff2ab3
# ╠═8d0d5940-6858-11eb-1155-dbe742a87729
# ╠═e2b38092-6858-11eb-1305-115e60d9ce26
# ╠═123426e4-6859-11eb-1857-150829245eca
# ╠═3882e254-6859-11eb-19d4-bfb9921902f2
# ╠═3f8cfd78-6859-11eb-07c2-8f02f60fc0a3
# ╠═45a9e11c-6859-11eb-3266-cb2f4081e7de
# ╠═10d6c6a2-685a-11eb-269d-3f92f1a94fe1
# ╟─8b681600-685a-11eb-0bfd-bbaf19e1a22b
# ╠═b8855382-685a-11eb-2bd4-7f7d3557d836
# ╠═e2e282da-685a-11eb-32b2-116fce22d9d3
# ╠═e71d56b8-685a-11eb-367c-c37ad541bdf4
# ╠═0f3277dc-685b-11eb-2bc5-d7cadea01388
# ╠═9fb53268-685b-11eb-18cc-5d103022febf
# ╟─0c9579e2-685c-11eb-0985-b323f500064c
# ╠═1816ffd4-685c-11eb-218c-3985c1e8fd40
# ╠═686aed92-685c-11eb-295b-81d603e5eedc
# ╠═9a573124-685c-11eb-107a-05702bbfc2f2
# ╠═a0855ff0-685c-11eb-3b70-b7c6d257b1bf
# ╠═bd6bf02a-685c-11eb-1b01-d9dfda6e8f91
# ╠═db3fe086-685c-11eb-0a28-e7529365d7b8
# ╟─f288a028-685c-11eb-07f2-3de80f49b384
# ╠═ff22edf2-685c-11eb-226e-c9470f8902e9
# ╠═2304a31e-685d-11eb-2c68-037fefdb684f
# ╠═38b95e48-685d-11eb-3c4a-8fd1e28d00a6
# ╠═659a12a2-685d-11eb-1c32-6b53b431ac9a
# ╠═b3775770-685d-11eb-3090-dbbc1689f1d0
# ╠═e5ae84de-685d-11eb-2377-918cd4b53af4
# ╠═f2fb7c8c-685d-11eb-0def-01b7643c6547
# ╠═4ced4446-685e-11eb-3431-c3de701b63a6
# ╠═dcc9489e-685e-11eb-315e-e34875afae1e
# ╟─fb37bad6-685e-11eb-0314-73160a66b806

max_stars_repo_path: Julia Kata/07_alphametics_solver.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Julia Kata/07_alphametics_solver.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Julia Kata/07_alphametics_solver.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
"""
Reference: The Annotated Transformer (http://nlp.seas.harvard.edu/2018/04/03/attention.html)
"""

using Flux
using Flux: onecold, gradient
import Flux.Optimise: update!

using WordTokenizers

using Transformers
using Transformers.Basic

include("./0-data.jl")

const vocab = Vocabulary(labels, unksym)
const embed = gpu(Embed(512, length(vocab)))

embedding(x) = embed(x, inv(sqrt(512)))

const encoder = gpu(Stack(
    @nntopo(e → pe:(e, pe) → x → x → $N),
    PositionEmbedding(512),
    (e, pe) -> e .+ pe,
    Dropout(0.1),
    [Transformer(512, 8, 64, 2048) for i = 1:N]...
))

const decoder = gpu(Stack(
    @nntopo((e, m, mask):e → pe:(e, pe) → t → (t:(t, m, mask) → t:(t, m, mask)) → $N:t → c),
    PositionEmbedding(512),
    (e, pe) -> e .+ pe,
    Dropout(0.1),
    [TransformerDecoder(512, 8, 64, 2048) for i = 1:N]...,
    Positionwise(Dense(512, length(labels)), logsoftmax)
))

const ps = params(embed, encoder, decoder)
const opt = ADAM(lr)

function smooth(et)
    global Smooth
    sm = fill!(similar(et, Float32), Smooth/size(embed, 2))
    p = sm .* (1 .+ -et)
    label = p .+ et .* (1 - convert(Float32, Smooth))
    label
end

function loss(src, trg, src_mask, trg_mask)
    lab = onehot(vocab, trg)

    src = embedding(src)
    trg = embedding(trg)

    if src_mask === nothing || trg_mask === nothing
        mask = nothing
    else
        mask = getmask(src_mask, trg_mask)
    end

    enc = encoder(src)
    dec = decoder(trg, enc, mask)

    #label smoothing
    label = smooth(lab)[:, 2:end, :]

    if mask == nothing
        loss = logkldivergence(label, dec[:, 1:end-1, :])
    else
        loss = logkldivergence(label, dec[:, 1:end-1, :], trg_mask[:, 1:end-1, :])
    end
end

function translate(x)
    ix = todevice(vocab(mkline(x)))
    seq = [startsym]

    src = embedding(ix)
    enc = encoder(src)

    len = length(ix)
    for i = 1:2len
        trg = embedding(todevice(vocab(seq)))
        dec = decoder(trg, enc, nothing)
        #move back to gpu due to argmax wrong result on CuArrays
        ntok = onecold(collect(dec), labels)
        push!(seq, ntok[end])
        ntok[end] == endsym && break
    end
    seq
end

train!()


text:
# poor man's info dict for testing
info_dict(MM::Type{<:Model}) =
    Dict(trait => eval(:($trait))(MM) for trait in M.MODEL_TRAITS)

@mlj_model mutable struct FooRegressor <: Deterministic
    a::Int = 0::(_ ≥ 0)
    b
end

metadata_pkg(FooRegressor,
    name="FooRegressorPkg",
    uuid="10745b16-79ce-11e8-11f9-7d13ad32a3b2",
    url="http://existentialcomics.com/",
    julia=true,
    license="MIT",
    is_wrapper=false
    )
metadata_model(FooRegressor,
    input=Table(Continuous),
    target=AbstractVector{Continuous},
    descr="La di da")

@testset "metadata" begin
    setfull()
    M.implemented_methods(::FI, M::Type{<:MLJType}) =
        getfield.(methodswith(M), :name)
    infos = info_dict(FooRegressor)

    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:output_scitype] == Unknown
    @test infos[:target_scitype] == AbstractVector{Continuous}
    @test infos[:is_pure_julia]
    @test infos[:package_name] == "FooRegressorPkg"
    @test infos[:package_license] == "MIT"
    @test infos[:load_path] == "MLJModels.FooRegressorPkg_.FooRegressor"
    @test infos[:package_uuid] == "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
    @test infos[:package_url] == "http://existentialcomics.com/"
    @test !infos[:is_wrapper]
    @test !infos[:supports_weights]
    @test !infos[:supports_online]
    @test startswith(infos[:docstring], "La di da")
    @test infos[:name] == "FooRegressor"
    @test infos[:is_supervised]
    @test infos[:prediction_type] == :deterministic
    @test infos[:implemented_methods] == [:clean!]
    @test infos[:hyperparameters] == (:a, :b)
    @test infos[:hyperparameter_types] == ("Int64", "Any")
    @test infos[:hyperparameter_ranges] == (nothing, nothing)
end


text:
# These are some macros (and supporting functions) to make it easier to define rules.

"""
    @scalar_rule(f(x₁, x₂, ...),
                 @setup(statement₁, statement₂, ...),
                 (∂f₁_∂x₁, ∂f₁_∂x₂, ...),
                 (∂f₂_∂x₁, ∂f₂_∂x₂, ...),
                 ...)

A convenience macro that generates simple scalar forward or reverse rules using
the provided partial derivatives. Specifically, generates the corresponding
methods for `frule` and `rrule`:

    function ChainRulesCore.frule(::typeof(f), x₁::Number, x₂::Number, ...)
        Ω = f(x₁, x₂, ...)
        \$(statement₁, statement₂, ...)
        return Ω, (Rule((Δx₁, Δx₂, ...) -> ∂f₁_∂x₁ * Δx₁ + ∂f₁_∂x₂ * Δx₂ + ...),
                   Rule((Δx₁, Δx₂, ...) -> ∂f₂_∂x₁ * Δx₁ + ∂f₂_∂x₂ * Δx₂ + ...),
                   ...)
    end

    function ChainRulesCore.rrule(::typeof(f), x₁::Number, x₂::Number, ...)
        Ω = f(x₁, x₂, ...)
        \$(statement₁, statement₂, ...)
        return Ω, (Rule((ΔΩ₁, ΔΩ₂, ...) -> ∂f₁_∂x₁ * ΔΩ₁ + ∂f₂_∂x₁ * ΔΩ₂ + ...),
                   Rule((ΔΩ₁, ΔΩ₂, ...) -> ∂f₁_∂x₂ * ΔΩ₁ + ∂f₂_∂x₂ * ΔΩ₂ + ...),
                   ...)
    end

If no type constraints in `f(x₁, x₂, ...)` within the call to `@scalar_rule` are
provided, each parameter in the resulting `frule`/`rrule` definition is given a
type constraint of `Number`.
Constraints may also be explicitly be provided to override the `Number` constraint,
e.g. `f(x₁::Complex, x₂)`, which will constrain `x₁` to `Complex` and `x₂` to
`Number`.

Note that the result of `f(x₁, x₂, ...)` is automatically bound to `Ω`. This
allows the primal result to be conveniently referenced (as `Ω`) within the
derivative/setup expressions.

Note that the `@setup` argument can be elided if no setup code is need. In other
words:

    @scalar_rule(f(x₁, x₂, ...),
                 (∂f₁_∂x₁, ∂f₁_∂x₂, ...),
                 (∂f₂_∂x₁, ∂f₂_∂x₂, ...),
                 ...)

is equivalent to:

    @scalar_rule(f(x₁, x₂, ...),
                 @setup(nothing),
                 (∂f₁_∂x₁, ∂f₁_∂x₂, ...),
                 (∂f₂_∂x₁, ∂f₂_∂x₂, ...),
                 ...)

For examples, see ChainRulesCore' `rules` directory.

See also: [`frule`](@ref), [`rrule`](@ref), [`AbstractRule`](@ref)
"""
macro scalar_rule(call, maybe_setup, partials...)
    if Meta.isexpr(maybe_setup, :macrocall) && maybe_setup.args[1] == Symbol("@setup")
        setup_stmts = map(esc, maybe_setup.args[3:end])
    else
        setup_stmts = (nothing,)
        partials = (maybe_setup, partials...)
    end
    @assert Meta.isexpr(call, :call)
    f = esc(call.args[1])
    # Annotate all arguments in the signature as scalars
    inputs = map(call.args[2:end]) do arg
        esc(Meta.isexpr(arg, :(::)) ? arg : Expr(:(::), arg, :Number))
    end
    # Remove annotations and escape names for the call
    for (i, arg) in enumerate(call.args)
        if Meta.isexpr(arg, :(::))
            call.args[i] = esc(first(arg.args))
        else
            call.args[i] = esc(arg)
        end
    end
    if all(Meta.isexpr(partial, :tuple) for partial in partials)
        input_rep = :(first(promote($(inputs...))))  # stand-in with the right type for an input
        forward_rules = Any[rule_from_partials(input_rep, partial.args...) for partial in partials]
        reverse_rules = Any[]
        for i in 1:length(inputs)
            reverse_partials = [partial.args[i] for partial in partials]
            push!(reverse_rules, rule_from_partials(inputs[i], reverse_partials...))
        end
    else
        @assert length(inputs) == 1 && all(!Meta.isexpr(partial, :tuple) for partial in partials)
        forward_rules = Any[rule_from_partials(inputs[1], partial) for partial in partials]
        reverse_rules = Any[rule_from_partials(inputs[1], partials...)]
    end
    forward_rules = length(forward_rules) == 1 ? forward_rules[1] : Expr(:tuple, forward_rules...)
    reverse_rules = length(reverse_rules) == 1 ? reverse_rules[1] : Expr(:tuple, reverse_rules...)
    return quote
        function ChainRulesCore.frule(::typeof($f), $(inputs...))
            $(esc(:Ω)) = $call
            $(setup_stmts...)
            return $(esc(:Ω)), $forward_rules
        end
        function ChainRulesCore.rrule(::typeof($f), $(inputs...))
            $(esc(:Ω)) = $call
            $(setup_stmts...)
            return $(esc(:Ω)), $reverse_rules
        end
    end
end

function rule_from_partials(input_arg, ∂s...)
    wirtinger_indices = findall(x -> Meta.isexpr(x, :call) && x.args[1] === :Wirtinger,  ∂s)
    ∂s = map(esc, ∂s)
    Δs = [Symbol(string(:Δ, i)) for i in 1:length(∂s)]
    Δs_tuple = Expr(:tuple, Δs...)
    if isempty(wirtinger_indices)
        ∂_mul_Δs = [:(@thunk($(∂s[i])) * $(Δs[i])) for i in 1:length(∂s)]
        return :(Rule($Δs_tuple -> +($(∂_mul_Δs...))))
    else
        ∂_mul_Δs_primal = Any[]
        ∂_mul_Δs_conjugate = Any[]
        ∂_wirtinger_defs = Any[]
        for i in 1:length(∂s)
            if i in wirtinger_indices
                Δi = Δs[i]
                ∂i = Symbol(string(:∂, i))
                push!(∂_wirtinger_defs, :($∂i = $(∂s[i])))
                ∂f∂i_mul_Δ = :(wirtinger_primal($∂i) * wirtinger_primal($Δi))
                ∂f∂ī_mul_Δ̄ = :(conj(wirtinger_conjugate($∂i)) * wirtinger_conjugate($Δi))
                ∂f̄∂i_mul_Δ = :(wirtinger_conjugate($∂i) * wirtinger_primal($Δi))
                ∂f̄∂ī_mul_Δ̄ = :(conj(wirtinger_primal($∂i)) * wirtinger_conjugate($Δi))
                push!(∂_mul_Δs_primal, :($∂f∂i_mul_Δ + $∂f∂ī_mul_Δ̄))
                push!(∂_mul_Δs_conjugate, :($∂f̄∂i_mul_Δ + $∂f̄∂ī_mul_Δ̄))
            else
                ∂_mul_Δ = :(@thunk($(∂s[i])) * $(Δs[i]))
                push!(∂_mul_Δs_primal, ∂_mul_Δ)
                push!(∂_mul_Δs_conjugate, ∂_mul_Δ)
            end
        end
        primal_rule = :(Rule($Δs_tuple -> +($(∂_mul_Δs_primal...))))
        conjugate_rule = :(Rule($Δs_tuple -> +($(∂_mul_Δs_conjugate...))))
        return quote
            $(∂_wirtinger_defs...)
            AbstractRule(typeof($input_arg), $primal_rule, $conjugate_rule)
        end
    end
end


text:
using JuMP
using Clp

m = Model(Clp.Optimizer)

@variable(m, x₁ ≥ 0)
@variable(m, x₂ ≥ 0)

@objective(m, Max, 3000x₁  + 5000x₂)
@constraint(m, 3x₁ + 2x₂ ≤ 18)
@constraint(m, x₁ ≤ 4)
@constraint(m, x₂ ≤ 6)

print(m)

optimize!(m)

status = termination_status(m)

println("Solution status: ", status)

println("Objective value: ", objective_value(m))
println("x₁ = ", value(x₁))
println("x₂ = ", value(x₂))