text:
# === PRIVATE FUNCTIONS THAT ARE NOT EXPORTED =============================================================== #
# =========================================================================================================== #

# === PUBLIC FUNCTIONS THAT ARE EXPORTED ==================================================================== #
"""
    extract_minimum_complete_data_set(dataFrame::DataFrame, colNameArray::Array{String,1})::VLResult
"""
function extract_minimum_complete_data_set(dataFrame::DataFrame, colNameArray::Array{String,1})::VLResult

    # initialize -
    rows_to_keep_array = Array{Int64,1}()

    try

        # ok: grab the cols of interest -
        df_tmp = dataFrame[!,colNameArray]

        # logic: we need to go through each row, if we econter *any* missing value, then we need stop and throwthat row out
        (number_of_rows, number_of_cols) = size(df_tmp)
        for row_index = 1:number_of_rows

            # tmp -
            tmp_array = Array{Any,1}()
            for col_index = 1:number_of_cols
                push!(tmp_array,df_tmp[row_index,col_index])
            end

            # ok, so let's walk through -
            bit_array = ismissing.(tmp_array)

            # do we have *any* 1's?
            idx_one = findall(x->x==1,bit_array)
            if (length(idx_one)==0)
                push!(rows_to_keep_array,row_index)
            end
        end

        # ok: return filtered array -
        df_complete = dataFrame[rows_to_keep_array,colNameArray]

        # return -
        return VLResult(df_complete)
    catch error
        return VLResult(error)
    end

end

"""
load_study_data_set(filePath::String;
    removeMissingColumn::Bool = true, missingValueSentinal::String = "#N/A")::VLResult
"""
function load_study_data_set(filePath::String;
    removeMissingColumn::Bool = true, missingValueSentinal::String = "#N/A")::VLResult

    # initialize -
    col_to_keep_array = Array{Int64,1}()
    row_to_keep_array = Array{Int64,1}()

    try

        # check: is filePath ok?
        check_result = is_file_path_ok(filePath)
        if (isnothing(check_result.value) == false)
            return check_result
        end

        # ok, if we get here we can load the file -
        df = CSV.read(filePath, DataFrame)

        # what are the original names?
        col_name_array = names(df)

        # go trhough the cols, are there #N/A
        df_missing = DataFrame(replace(convert(Matrix, df), "$(missingValueSentinal)"=>missing))

        # update the col headers, to set back the original names -
        rename!(df_missing, col_name_array)

        # if removeMissingColumn = true, then lets go through each col of the data set.
        # if *all* the values in the col are missing, then remove that col -
        if (removeMissingColumn == false)
            return VLResult(df_missing)
        end

        # process cols -
        (number_of_rows,number_of_cols) = size(df_missing)
        for col_index = 1:number_of_cols

            # grab the col -
            data_col = df_missing[!,col_index]

            # ok, so let's through -
            bit_array = ismissing.(data_col)
            idx_one = findall(x->x==1,bit_array)

            # if the length(idx_one) => the number of rows, then all rows are missing
            if (length(idx_one) != number_of_rows)
                push!(col_to_keep_array, col_index)
            end
        end

        # grab the no-missing cols -
        df_populated_cols = df_missing[!,col_to_keep_array]

        # ok, last - let's remove rows that are all missing -
        (number_of_rows,number_of_cols) = size(df_populated_cols)
        for row_index = 1:number_of_rows

            # grab the row -
            data_row = df_populated_cols[row_index,:]

            # check -
            for col_index = 1:number_of_cols
                data_row_value = data_row[col_index]
                if (ismissing(data_row_value) == false)
                    push!(row_to_keep_array, row_index)
                    break
                end
            end
        end

        # grab the non-missing rows -
        df_populated_row_and_cols = df_populated_cols[row_to_keep_array,:]

        # return -
        return VLResult(df_populated_row_and_cols)
    catch error
        return VLResult(error)
    end
end
# =========================================================================================================== #

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit


text:
using MOOSE

include("ADCoupledConvection.jl")

# Create the Mesh
mesh = buildSquare(0, 1, 0, 1, 20, 20)

# Create the System to hold the equations
# By using Dual here the Jacobian will be computed
# using automatic differentiation.
# The "8" is for the maximum number of DoFs on
# any one element.
# Since we're using Lagrange with Quads with Two
# variables... that's how we get to 8
diffusion_system = System{Dual{8,Float64}}(mesh)

# Add some variables
u = addVariable!(diffusion_system, "u")
v = addVariable!(diffusion_system, "v")

# Apply the Laplacian operator to the "u" variable
addKernel!(diffusion_system, Diffusion(u))

# Apply the Coupled Convection operator
addKernel!(diffusion_system, ADCoupledConvection(u, v))

# Apply the Laplacian operator to the "v" variable
addKernel!(diffusion_system, Diffusion(v))

# u = 0 on the Left
addBC!(diffusion_system, DirichletBC(u, [4], 0.0))

# u = 1 on the Right
addBC!(diffusion_system, DirichletBC(u, [2], 1.0))

# v = 0 on the Left
addBC!(diffusion_system, DirichletBC(v, [4], 0.0))

# v = 1 on the Right
addBC!(diffusion_system, DirichletBC(v, [2], 1.0))

# Initialize the system of equations
initialize!(diffusion_system)

# Create a solver and solve
solver = JuliaDenseNonlinearImplicitSolver(diffusion_system)
solve!(solver, nl_max_its=5)

# Output
out = VTKOutput()
output(out, solver, "ad_coupled_convection_out")

ext: jl


text:
## Experimental syntax for indexing within vectors:
using Test

## use # to include
tpl = mt"{{#:vec}}{{#.[1]}}{{.}}{{/.[1]}}{{/:vec}}"
out = Mustache.render(tpl, vec=["A1","B2","C3"])
@test out == "A1"

## use ^ to exclude
tpl = mt"{{#:vec}}{{.}}{{^.[end]}}, {{/.[end]}}{{/:vec}}"
out = Mustache.render(tpl, vec=["A1","B2","C3"])
@test out == "A1, B2, C3"


text:
### Variables for Thermal Generation ####

"""
This function add the variables for power generation output to the model
"""
function activepowervariables(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64) where {A <: JumpExpressionMatrix, T <: PowerSystems.ThermalGen}

    on_set = [d.name for d in devices if d.available == true]

    t = 1:time_periods

    p_th = @variable(m, p_th[on_set,t], start = 0.0) # Power output of generators

    return p_th
end

"""
This function add the variables for power generation output to the model
"""
function reactivepowervariables(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64) where {A <: JumpExpressionMatrix, T <: PowerSystems.ThermalGen}

    on_set = [d.name for d in devices if d.available == true]

    t = 1:time_periods

    q_th = @variable(m, q_th[on_set,t], start = 0.0) # Power output of generators

    return q_th
end

"""
This function add the variables for power generation commitment to the model
"""
function commitmentvariables(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64) where T <: PowerSystems.ThermalGen

    on_set = [d.name for d in devices if d.available == true]

    t = 1:time_periods

    @variable(m, on_th[on_set,t], Bin) # Power output of generators
    @variable(m, start_th[on_set,t], Bin) # Power output of generators
    @variable(m, stop_th[on_set,t], Bin) # Power output of generators

end


text:
# arithmétique de virgule fixe
# precision q = 65536 i.e. 16 bits pour la partie décimale

function fp_add(x, y)
    x + y
end
function fp_sub(x, y)
    x - y
end
function fp_mul(x, y)
    t = x * y
    div(t + div(65536,2), 65536)
end
function fp_div(x, y)
    t = x * 65536
    div(t + div(y,2), y)
end
function fp_of_int(x)
    x * 65536
end

function iter(n, a, b, xn, yn)
    if n == 100
        return 1
    end
    xn2 = fp_mul(xn, xn)
    yn2 = fp_mul(yn, yn)
    if fp_add(xn2, yn2) > fp_of_int(4)
	0
    else
	iter(n+1, a, b, fp_add(fp_sub(xn2, yn2), a),
	     fp_add(fp_mul(fp_of_int(2), fp_mul(xn, yn)), b))
    end
end

function inside(x, y)
    iter(0, x, y, fp_of_int(0), fp_of_int(0)) > 0
end

function run(steps)
    xmin = fp_of_int(-2)
    xmax = fp_of_int(1)
    deltax = fp_div(fp_sub(xmax, xmin), fp_of_int(2*steps))
    ymin = fp_of_int(-1)
    ymax = fp_of_int(1)
    deltay = fp_div(fp_sub(ymax, ymin), fp_of_int(steps))
    for i = 0 : steps-1
	y = fp_add(ymin, fp_mul(fp_of_int(i), deltay))
	for j = 0 : 2*steps-1
	    x = fp_add(xmin, fp_mul(fp_of_int(j), deltax))
	    if inside(x, y)
		print("0")
	    else
		print("1")
	    end
	end
	print("\n")
    end
end

run(30)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
#=
Comparing two numbers written in index form like 211 and 37 is not difficult, as any calculator would confirm that 211 = 2048 < 37 = 2187.

However, confirming that 632382518061 > 519432525806 would be much more difficult, as both numbers contain over three million digits.

Using base_exp.txt (right click and 'Save Link/Target As...'), a 22K text file containing one thousand lines with a base/exponent pair on each line, determine which line number has the greatest numerical value.

NOTE: The first two lines in the file represent the numbers in the example given above.
=#
using DelimitedFiles
function calc()
  combined = readdlm(dirname(@__FILE__()) * "/../Resources/p99.txt", ',', Int64)
  bases = vec(combined[1:1000])
  exponents = vec(combined[1001:2000])
  max, line = 0,0
  for i in 1:1000
    result = exponents[i]*log10(bases[i])
    if result > max
      max, line = result, i
    end
  end
  line
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using HDF5
using JLD
using DataFrames

filename = joinpath(tempdir(), "mydata.jld")

df = DataFrame({[2:6], pi*[1:5]})
df2 = @DataFrame(a => [1:5], b => pi * [1:5])

file = jldopen(filename, "w")
write(file, "df", df)
write(file, "df2", df2)
close(file)

file = jldopen(filename, "r")
x = read(file, "df")
y = read(file, "df2")
close(file)

using Base.Test
@test all(df .== x)
@test all(df2 .== y)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export numberinput

"""
    numberinput(;kwargs...)

A NumberInput component.
NumberInput component
Keyword arguments:
- `allowEmpty` (Bool; optional): `true` to allow empty string.
- `style` (Dict; optional): style of the input
- `ariaLabel` (String; optional): Provide a description that would be used to best describe the use case of the NumberInput component
- `className` (String; optional): Specify an optional className to be applied to the wrapper node
- `defaultValue` (Real | String; optional): Optional starting value for uncontrolled state
- `disabled` (Bool; optional): Specify if the control should be disabled, or not
- `helperText` (String; optional): Provide text that is used alongside the control label for additional help
- `hideLabel` (Bool; optional): Specify whether you want the underlying label to be visually hidden
- `iconDescription` (String; optional): Provide a description for up/down icons that can be read by screen readers
- `id` (String; required): Specify a custom `id` for the input
- `invalid` (Bool; optional): Specify if the currently value is invalid.
- `invalidText` (String; optional): Message which is displayed if the value is invalid.
- `label` (String; optional): Generic `label` that will be used as the textual representation of what
this field is for
- `light` (Bool; optional): `true` to use the light version.
- `max` (Real; optional): The maximum value.
- `min` (Real; optional): The minimum value.
- `readOnly` (Bool; optional): Specify if the component should be read-only
- `size` (a value equal to: 'sm', 'xl'; optional): Specify the size of the Number Input. Currently supports either `sm` or `xl` as an option.
- `step` (Real; optional): Specify how much the valus should increase/decrease upon clicking on up/down button
- `value` (Real | String; optional): Specify the value of the input
- `warn` (Bool; optional): Specify whether the control is currently in warning state
- `warnText` (String; optional): Provide the text that is displayed when the control is in warning state
"""
function numberinput(; kwargs...)
        available_props = Symbol[:allowEmpty, :style, :ariaLabel, :className, :defaultValue, :disabled, :helperText, :hideLabel, :iconDescription, :id, :invalid, :invalidText, :label, :light, :max, :min, :readOnly, :size, :step, :value, :warn, :warnText]
        wild_props = Symbol[]
        return Component("numberinput", "NumberInput", "dash_carbon_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
using Flux
using AlphaGo
using AlphaGo: select_leaf, incorporate_results!, child_U, inject_noise!,
                Q, N, child_Q, initialize_game!, tree_search!, MCTSRules, extract_data
using AlphaGo: pass_move!, play_move!, PlayerMove, BLACK, WHITE, PlayerMove,
                to_kgs, score, from_kgs

include("test_utils.jl")

struct DummyNet
  fake_priors
  fake_value

  function DummyNet(env; fake_priors = nothing, fake_value = 0)
    if fake_priors == nothing
      fake_priors = ones(env.action_space) / (env.action_space)
    end
    new(fake_priors, fake_value)
  end
end

(dn::DummyNet)(position::Position) = param(dn.fake_priors),
                                          param(dn.fake_value)

function (dn::DummyNet)(positions = nothing)
  if positions == nothing || length(positions) == 0
    throw(ArgumentError("No positions passed!"))
  end
  len = length(positions)
  param(repeat(dn.fake_priors, outer = [1, len])),
          param(repeat([dn.fake_value], len))
end

env = GoEnv(9)

ALMOST_DONE_BOARD = load_board("""
                              .XO.XO.OO
                              X.XXOOOO.
                              XXXXXOOOO
                              XXXXXOOOO
                              .XXXXOOO.
                              XXXXXOOOO
                              .XXXXOOO.
                              XXXXXOOOO
                              XXXXOOOOO
                              """, env)

SEND_TWO_RETURN_ONE = Position(env,
    board = ALMOST_DONE_BOARD,
    n = 70,
    komi = 2.5,
    caps = (1, 4),
    ko = nothing,
    recent = [PlayerMove(BLACK, (1, 2)),
    PlayerMove(WHITE, (1, 9))],
    to_play = BLACK
    );

function initialize_basic_player(env)
  player = MCTSPlayer(env, DummyNet(env))
  initialize_game!(player)
  first_node = select_leaf(player.root)
  p, v = player.network(player.root.position)
  incorporate_results!(first_node, p.data, v.data, player.root)
  return player
end

function initialize_almost_done_player(env)
  probs = ones(env.action_space) * 0.001
  probs[3:5] .= 0.2  # some legal moves along the top.
  probs[end] = 0.2  # passing is also ok
  net = DummyNet(env; fake_priors = probs)
  player = MCTSPlayer(env, net)
  # root position is white to play with no history == white passed.
  initialize_game!(player, SEND_TWO_RETURN_ONE)
  return player
end

@testset "MCTSPlayer" begin
  # Tromp taylor means black can win if we hit the move limit.
  TT_FTW_BOARD = load_board("""
                            .XXOOOOOO
                            X.XOO...O
                            .XXOO...O
                            X.XOO...O
                            .XXOO..OO
                            X.XOOOOOO
                            .XXOOOOOO
                            X.XXXXXXX
                            XXXXXXXXX
                            """, env)

  @testset "inject_noise" begin
    player = initialize_basic_player(env)
    sum_priors = sum(player.root.child_prior)
    # dummyNet should return normalized priors.
    @test sum_priors ≈ 1
    u = child_U(player.root)
    @test all(u .== u[1])

    inject_noise!(player.root)
    new_sum_priors = sum(player.root.child_prior)
    # priors should still be normalized after injecting noise
    @test sum_priors ≈ new_sum_priors

    # With dirichelet noise, majority of density should be in one node.
    max_p = maximum(player.root.child_prior)
    @test max_p > 3 / (env.action_space)
  end

  @testset "pick_moves" begin
    player = initialize_basic_player(env)
    root = player.root
    root.child_N[to_flat((3, 1), player.root.position)] = 10
    root.child_N[to_flat((2, 1), player.root.position)] = 5
    root.child_N[to_flat((4, 1), player.root.position)] = 1

    root.position.n = env.action_space  # move 81, or 361, or... Endgame.

    # Assert we're picking deterministically
    @test root.position.n > player.τ_threshold
    move = pick_move(player)
    @test move == (3, 1)

    # But if we're in the early part of the game, pick randomly
    root.position.n = 3
    @test player.root.position.n ≤ player.τ_threshold

    #TODO: complete this test
    #with mock.patch('random.random', lambda: .5)
    #move = pick_move(player)
    #@test move == (3, 1)

    #with mock.patch('random.random', lambda: .99):
    #move = pick_move(player)
    #@test move == (4, 1)
  end

 @testset "dont_pass_if_losing" begin
    player = initialize_almost_done_player(env)

    # check -- white is losing.
    @test score(player.root.position) == -0.5

    for i = 1:20
      tree_search!(player)
    end

    # uncomment to debug this test
    #println(describe(player.root))

    # Search should converge on D9 as only winning move.
    flattened = to_flat(from_kgs("D9", env), player.root.position)
    best_move = findmax(player.root.child_N)[2]
    @test best_move == flattened
    # D9 should have a positive value
    @test Q(player.root.children[flattened]) > 0
    @test N(player.root) ≥ 20
    # passing should be ineffective.
    @test child_Q(player.root)[end] < 0
    # no virtual losses should be pending
    @test assertNoPendingVirtualLosses(player.root)
    # uncomment to debug this test
    #println(describe(player.root))
  end

  @testset "parallel_tree_search" begin
    player = initialize_almost_done_player(env)
    # check -- white is losing.
    @assert score(player.root.position) == -0.5
    # initialize the tree so that the root node has populated children.
    tree_search!(player, 1)
    # virtual losses should enable multiple searches to happen simultaneously
    # without throwing an error...
    for i = 1:6
      tree_search!(player, 10)
    end
    # uncomment to debug this test
    # print(player.root.describe())

    # Search should converge on D9 as only winning move.
    flattened = to_flat(from_kgs("D9", env), player.root.position)
    best_moves = findall(x -> x .== maximum(player.root.child_N), player.root.child_N)
    @test flattened ∈ best_moves
    # D9 should have a positive value
    @test Q(player.root.children[flattened]) > 0
    @test N(player.root) ≥ 20
    # passing should be ineffective.
    child_Q(player.root)[end] < 0
    # no virtual losses should be pending
    @test assertNoPendingVirtualLosses(player.root)
  end

  @testset "ridiculously_parallel_tree_search" begin
    player = initialize_almost_done_player(env)
    # Test that an almost complete game
    # will tree search with # parallelism > # legal moves.
    for i = 1:10
      tree_search!(player, 50)
    end
    @test assertNoPendingVirtualLosses(player.root)
  end

  @testset "long_game_tree_search" begin
    mcts_rules = MCTSRules(env)
    player = MCTSPlayer(env, DummyNet(env))
    endgame = Position(env,
        board = TT_FTW_BOARD,
        n = mcts_rules.max_game_length - 2,
        komi = 2.5,
        ko = nothing,
        recent = [PlayerMove(BLACK, (1, 2)),
                PlayerMove(WHITE, (1, 9))],
        to_play = BLACK
    )
    initialize_game!(player, endgame)

    # Test that MCTS can deduce that B wins because of TT-scoring
    # triggered by move limit.
    for i = 1:10
      tree_search!(player, 8)
    end
    @test assertNoPendingVirtualLosses(player.root)
    @test Q(player.root) > 0
  end

  @testset "cold_start_parallel_tree_search" begin
    # Test that parallel tree search doesn't trip on an empty tree
    player = MCTSPlayer(env, DummyNet(env, fake_value = 0.17))
    initialize_game!(player)
    @test N(player.root) == 0
    @test !player.root.is_expanded
    tree_search!(player, 4)
    @test assertNoPendingVirtualLosses(player.root)
    # Even though the root gets selected 4 times by tree search, its
    # final visit count should just be 1.
    N(player.root) == 1
    # 0.085 = average(0, 0.17), since 0 is the prior on the root.
    @test Q(player.root) ≈ 0.085
  end

  @testset "tree_search_failsafe" begin
    # Test that the failsafe works correctly. It can trigger if the MCTS
    # repeatedly visits a finished game state.
    probs = ones(env.N * env.N + 1) * 0.001
    probs[end] = 1  # Make the dummy net always want to pass
    player = MCTSPlayer(env, DummyNet(env, fake_priors = probs))
    pass_position = pass_move!(Position(env))
    initialize_game!(player, pass_position)
    tree_search!(player, 1)
    @test assertNoPendingVirtualLosses(player.root)
  end

  @testset "only_check_game_end_once" begin
    # When presented with a situation where the last move was a pass,
    # and we have to decide whether to pass, it should be the first thing
    # we check, but not more than that.

    white_passed_pos = pass_move!(
                        play_move!(
                          play_move!(
                            play_move!(
                              Position(env), (4,4) # b plays
                                ), (4,5)  # w plays
                              ), (5,4)  # b plays
                            ) # w passes - if B passes too, B would lose by komi.
                          )

    player = MCTSPlayer(env, DummyNet(env))
    initialize_game!(player, white_passed_pos)
    # initialize the root
    for i = 1:15
      tree_search!(player)
    end
    # explore a child - should be a pass move.
    #tree_search!(player, 10)
    pass_move = env.N * env.N + 1
    @test N(player.root.children[pass_move]) == 1
    @test player.root.child_N[pass_move] == 1
    tree_search!(player)
    # check that we didn't visit the pass node any more times.
    @test player.root.child_N[pass_move] == 1
  end

  @testset "extract_data_normal_end" begin
    player = MCTSPlayer(env, DummyNet(env))
    initialize_game!(player)
    tree_search!(player)
    play_move!(player, nothing)
    tree_search!(player)
    play_move!(player, nothing)
    @test is_done(player.root)
    set_result!(player, AlphaGo.result(player.root.position), false)

    positions, πs, results = extract_data(player)
    @test length(positions) == length(πs) == length(results) == 2
    position, pi, result = positions[1], πs[1], results[1]
    # White wins by komi
    @test result == WHITE
    @test player.result_string == "W+$(player.root.position.komi)"
  end

  @testset "extract_data_resign_end" begin
    player = MCTSPlayer(env, DummyNet(env))
    initialize_game!(player)
    tree_search!(player)
    play_move!(player, (1, 1))
    tree_search!(player)
    play_move!(player, nothing)
    tree_search!(player)
    # Black is winning on the board
    @test AlphaGo.result(player.root.position) == BLACK
    # But if Black resigns
    set_result!(player, WHITE, true)

    data = extract_data(player)
    position, pi, result = data[1], data[2], data[3]
    # Result should say White is the winner
    @test result[1] == WHITE
    @test player.result_string == "W+R"
  end
end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    m = 15
    n = 5
    data = rand(m,n)
    # some fake ticks
    alphabet = 'A':'E'
    yticks = string.(collect(alphabet))
    k = 4
    itr = Iterators.product(ntuple(_ -> alphabet, k)...)
    xticks = []
    for word in Base.Generator(join, itr)
        push!(xticks, word)
        if length(xticks) == m
            break
        end
    end
    xticks = string.(xticks)

    fig = Figure(resolution = (1200, 600), fontsize = 20)
    ax = Axis(fig)
    hmap = heatmap!(ax, data, colormap = :plasma)
    for i in 1:15, j in 1:5
        if data[i,j] < 0.15
            txtcolor = :white
        else
            txtcolor = :black
        end
        text!(ax, "$(round(data[i,j], digits = 2))", position = (i,j),
            color = txtcolor, align = (:center, :center))
    end

    cbar = Colorbar(fig, hmap, label = "values", width = 15,
                ticksize=15, tickalign = 1, height = Relative(3.55/4))
    ax.xticks = (1:m, xticks)
    ax.yticks = (1:n, yticks)
    ax.xticklabelrotation = π/3
    #ax.xticklabelsize = 12
    #ax.xticklabelalign = (:center, :center)
    #ax.xticklabelpad = 20
    fig[1, 1] = ax
    fig[1, 2] = cbar
    colgap!(fig.layout, 7)
    fig
    save(joinpath(@__DIR__, "output", "heatmapText.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using DifferentialEquations, Plots, Statistics

function Diffusion(μ0,μ1,κ,α,σ,ρ,u0,tspan)
    f = function (du,u,p,t)
        du[1] = μ0 + μ1*(u[2]-α)/σ # drift in log prices
        du[2] = κ*(α-u[2]) # mean reversion in shocks
    end
    g = function (du,u,p,t)
        du[1] = exp(u[2]/2.0)
        du[2] = σ
    end
    Γ = [1.0 ρ;ρ 1.0] # Covariance Matrix
    noise = CorrelatedWienerProcess!(Γ,tspan[1],zeros(2),zeros(2))
    sde_f = SDEFunction{true}(f,g)
    SDEProblem(sde_f,g,u0,tspan,noise=noise)
end

function dgp()
# trading days, closing, etc

# assume trading period is 1/3 of day (8 hours)
# but that latent price evolves continuously
# observed return is daily difference of log price at
# closing time
TradingDays = 1000 # total days in sample
Days = Int(ceil(TradingDays*7/5))+3  # calendar days
MinPerDay = 1440 # minutes per day
MinPerTic = 5 # minutes between tics, lower for better accuracy
tics = Int(MinPerDay/MinPerTic) # number of tics in day
dt = 1/tics # divisions per day
closing = Int(floor(tics/3)) # closing tic: closing happens after 1/3 of day
# parameters
μ0 = 0.0
μ1 = 0.0
κ = 0.1
α = 0.15
σ = 0.15
ρ = -0.7
λ0 = 0.005 # jump rate per day
λ1 = 4.0 # scaling factor st. dev. of jumps
u0 = [0.0;α]
prob = Diffusion(μ0, μ1, κ, α, σ, ρ, u0, (0.0,Days))

## jump in price
rate(u,p,t) = λ0
# jump is normal with st. dev. equal to λ1 times current st. dev.
affect1!(integrator) = (integrator.u[1] = integrator.u[1].+randn(size(integrator.u[1])).*λ1.*exp(integrator.u[2]./2.0))
#affect1!(integrator) = (integrator.u[1] = λ1*integrator.u[1])
jump = ConstantRateJump(rate,affect1!)
jump_prob = JumpProblem(prob,Direct(), jump)
sol = solve(jump_prob,SRIW1(), dt=dt, adaptive=false)
# get lnP at each tic
lnP = [sol(t)[1] for t in 0:dt:(Days-1)]
# get log price at end of trading days
z = zeros(TradingDays+1)
RV = zeros(TradingDays+1)
global DayofWeek = 0 # counter for day of week
global TradingDay = 0 # counter for trading days
global Day = 0
while TradingDay<=TradingDays
    # set day of week, and record if it's a trading day
    Day +=1
    DayofWeek +=1
    if DayofWeek<6
        TradingDay +=1 # advance trading day
        # compute realized measures
        for tic = 2:closing
            RV[TradingDay]+=(lnP[Day*tics+tic]-lnP[Day*tics+tic-1])^2.0
        end
        z[TradingDay] = lnP[Day*tics+closing]
    end
    if DayofWeek==7 # restart the week if Sunday
        DayofWeek = 0
    end
end
rets = z[2:end]-z[1:end-1]
RV = RV[2:end]

# checking the jump rate
t = sol.t
jump = t[2:end].==t[1:end-1] # find where time does not change
jump[1]=0 # there is alway a "jump" at beginning, drop it
println("number of jumps: ", sum(jump), " in ", TradingDays, " days")
println("empirical jump rate (?): ", sum(jump)/TradingDays)
println("nominal jump rate: ", λ0)
return z, jump, sol, rets, RV;
end
z, jump, sol, rets, RV = dgp();
t = sol.t;
t = t[2:end]
u = sol(t[jump]).u
jumpvals = [u[i][1] for i in 1:sum(jump)]
plot(rets, label="returns");
vline!(t[jump].*1000/1400, label="jumps")

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
#=
rect.jl
=#

export rect

using Test: @test, @inferred

"""
    rect(x::Real)
Unit width rect function. Potential problem? Bring up with fess.
"""
function rect(x::Real)
	T = promote_type(eltype(x), Float32)
    T(abs(x) < 0.5)
end

"""
    rect(:test)
self test
"""
function rect(test::Symbol)
	test != :test && throw("test $test")
	@test (@inferred rect(Int(3))) === zero(Float32)
	@test (@inferred rect(Int(0))) === one(Float32)
	@test (@inferred rect(0.25)) === one(Float64)
	@test (@inferred rect(Float16(0.25))) === one(Float32)
	true
end


text:
@benchgroup "mapwindow" begin
    img1d = randn(1000)
    img2d = randn(30,30)
    img3d = randn(10,11,12)
    @bench "cheap f, tiny window" mapwindow(first, img1d, (1,))
    @bench "extrema" mapwindow(extrema, img2d, (5,5))
    @bench "median!" mapwindow(median!, img2d, (5,5))
    @bench "mean, small window" mapwindow(mean, img1d, (3,))
    @bench "mean, large window" mapwindow(mean, img3d, (5,5,5))
    @bench "expensive f" mapwindow(x -> quantile(vec(x), 0.7), img3d, (3,3,3))
end


text:
__precompile__(false)
#
#
#
# A Shadowsocks implementation written in Julia
# Create by John Xiong <imgk@mail.ustc.edu.cn>
#
#
#
#
#
module Shadowsocks

include("Common.jl")
include("Server.jl")
include("Client.jl")
include("Obfs.jl")

using .Common: @terminal, SSConfig, readConfigFile, parseURI
using .Server: runServer
using .Client: runClient

import Base.run

function run(configs::Array{SSConfig}) # configure multi servers at client side
    runClient(configs)
end

function run(ch::Channel{Dict{String, Any}}) # configure multi servers at server side
    servers = Dict{String, Condition}()

    while true
        configs = take!(ch)

        key = keys(servers) # add server
        for (k, v) in configs
            if !(k in key)
                servers[k] = Condition()
                @async runServer(v, servers[k])
            end
        end

        key = keys(configs) # disable server
        for (k, ~) in servers
            if !(k in key)
                notify(servers[k])
                delete!(servers, k)
            end
        end
    end
end

end # module


text:
# try-catch
try
    throw(DomainError("그냥 에러"))
catch e
    println("예외 던져짐: $e")
end

# try-finally
try
    throw(DomainError("그냥 에러"))
finally
    println("종료됨")
end

max_stars_repo_name: myyrakle/Julia-Example

max_issues_repo_name: myyrakle/Julia-Example

max_forks_repo_name: myyrakle/Julia-Example


text:
function simulate!(X, mcProcess::finiteActivityProcess, rfCurve::AbstractZeroRateCurve, mcBaseData::MonteCarloConfiguration{type1, type2, type3, type4, type5}, T::numb) where {finiteActivityProcess <: FiniteActivityProcess, numb <: Number, type1 <: Number, type2 <: Number, type3 <: AbstractMonteCarloMethod, type4 <: BaseMode, type5 <: Random.AbstractRNG}
    r = rfCurve.r
    d = dividend(mcProcess)
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    λ = mcProcess.λ

    @assert T > 0.0

    ####Simulation
    ## Simulate
    # r-d-psi(-i)
    drift_RN = (r - d) - σ^2 / 2 - λ * compute_drift(mcProcess)
    simulate!(X, BrownianMotion(σ, drift_RN), mcBaseData, T)

    dt = T / Nstep
    PoissonRV = Poisson(λ * T)
    NJumps = rand(mcBaseData.rng, PoissonRV, Nsim)

    for (index_, njumps_) in enumerate(NJumps)
        # Simulate the times of jump
        @views tjumps = rand(mcBaseData.rng, njumps_) * T
        for tjump in tjumps
            # Add the jump size
            idx1 = ceil(Int64, tjump / dt) + 1
            #jump_size=compute_jump_size(mcProcess,mcBaseData);
            @views X[index_, idx1:end] .+= compute_jump_size(mcProcess, mcBaseData)
        end
    end
    ## Conclude
    S0 = mcProcess.underlying.S0
    #f(x)=S0*exp(x);
    X .= S0 .* exp.(X)
    #broadcast!(f,X,X)
    nothing
end


text:
function dataProvider()
    srand(100000)
    N=1000
    x=[1+rand(1,N) -1-rand(1,N)]
    y=[2*randn(1,N) 2*randn(1,N)]
    inputBase=[x;y]
    outputBase=label2vec([ones(1,N) 2*ones(1,N)],2)
    idx=randcycle(2*N)
    inputBase=inputBase[:,idx]
    outputBase=outputBase[:,idx]
  return inputBase,outputBase
end


text:
"""
`@style(object, args)` forces the displaying of `object` (even if it is not of a
  redirected type). `args` are pairs `parameter => string` for the styling of the
  display.
"""
macro style(obj, args...)
  _, opts = parseargs(args...)

  try
    obj = eval(current_module(), obj)
  catch
    error("cannot evaluate object $(obj)")
  end

  showmsg(obj, opts)
end


text:
using Documenter
using GeometricFlux

makedocs(
    sitename = "GeometricFlux",
    format = Documenter.HTML(
      canonical = "https://yuehhua.github.io/GeometricFlux.jl/stable"
    ),
    modules = [GeometricFlux],
    pages = ["Home" => "index.md",
             "Get started" => "start.md",
             "Basics" =>
               ["Building layers" => "basics/layers.md",
                "Graph passing" => "basics/passgraph.md"],
             "Abstractions" =>
               ["Message passing scheme" => "abstractions/msgpass.md",
                "Graph network block" => "abstractions/gn.md"],
             "Manual" =>
               ["Convolutional Layers" => "manual/conv.md",
                "Pooling Layers" => "manual/pool.md",
                "Models" => "manual/models.md",
                "Linear Algebra" => "manual/linalg.md",
                "Utilities" => "manual/utils.md"]
    ]
)

deploydocs(repo = "github.com/yuehhua/GeometricFlux.jl.git")


text:
using Documenter, Temporal

makedocs(modules=[Temporal],
         doctest=false,
         clean=true,
         format=Documenter.HTML(),
         authors="Jacob Amos",
         sitename="Temporal",
         pages=["Home" => "index.md",
                "Overview" => "overview.md",
                "Calculations" => ["Methods" => "calculation.md",
                                   "Aggregation" => "aggregation.md"],
                "Data" => ["I/O" => "io.md",
                           "Subsetting" => "indexing.md",
                           "Combining" => "combining.md"]])

deploydocs(deps=Deps.pip("mkdocs", "python-markdown-math"),
           repo="github.com/dysonance/Temporal.jl.git",
           devbranch="master",
           devurl="dev",
           versions=["stable" => "v^", "v#.#", "dev" => "dev"])


text:
module Database

const Q = 1.60217662e-19

const Rydberg = 13.6

const Material = Dict(
    "Si" => Dict(
        "latticeconstant" => 5.431e-10,
        "elasticconstant" => Dict(
            "C11" => 1660,
            "C12" => 640,
            "C44" => 796,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.21 * Q * Rydberg,
                4 => 0 * Q *  Rydberg,
                8 => 0.04 * Q * Rydberg,
                11 => 0.08 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.00 * Rydberg,
                4 => 0.00 * Rydberg,
                8 => 0.00 * Rydberg,
                11 => 0.00 * Rydberg,
            )
        )
    ),
    "Ge" => Dict(
        "latticeconstant" => 5.658e-10,
        "elasticconstant" => Dict(
            "C11" => 1260,
            "C12" => 440,
            "C44" => 677,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.230 * Q * Rydberg,
                4 => 0 * Q * Rydberg,
                8 => 0.01 * Q * Rydberg,
                11 => 0.06 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.00 * Rydberg,
                4 => 0.00 * Rydberg,
                8 => 0.00 * Rydberg,
                11 => 0.00 * Rydberg,
            )
        )
    ),
    "GaAs" => Dict(
        "latticeconstant" => 5.65325e-10,
        "elasticconstant" => Dict(
            "C11" => 1221,
            "C12" => 566,
            "C44" => 600,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.23 * Q * Rydberg,
                4 => 0.00 * Q * Rydberg,
                8 => 0.01 * Q * Rydberg,
                11 => 0.06 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.07 * Q * Rydberg,
                4 => 0.05 * Q * Rydberg,
                8 => 0.00 * Q * Rydberg,
                11 => 0.01 * Q * Rydberg,
            )
        ),

    ),
    "GaP" => Dict(
        "latticeconstant" => 5.4505e-10,
        "elasticconstant" => Dict(
            "C11" => 1405,
            "C12" => 620.3,
            "C44" => 703.3,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.22 * Q * Rydberg,
                4 => 0 * Q * Rydberg,
                8 => 0.03 * Q * Rydberg,
                11 => 0.07 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.12 * Q * Rydberg,
                4 => 0.07 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.02 * Q * Rydberg,
            )
        ),
    ),
    "InP" => Dict(
        "latticeconstant" => 5.8697e-10,
        "elasticconstant" => Dict(
            "C11" => 1011,
            "C12" => 561,
            "C44" => 456,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.23 * Q * Rydberg,
                4 => 0.00 * Q * Rydberg,
                8 => 0.01 * Q * Rydberg,
                11 => 0.06 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.07 * Q * Rydberg,
                4 => 0.05 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.01 * Q * Rydberg,
            )
        ),
    ),
    "InAs" => Dict(
        "latticeconstant" => 6.0583e-10,
        "elasticconstant" => Dict(
            "C11" => 832.9,
            "C12" => 452.6,
            "C44" => 395.9,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.22 * Q * Rydberg,
                4 => 0.00 * Q * Rydberg,
                8 => 0.00 * Q * Rydberg,
                11 => 0.05 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.08 * Q * Rydberg,
                4 => 0.05 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.03 * Q * Rydberg,
            )
        ),
    ),
    "InSb" => Dict(
        "latticeconstant" => 6.4794e-10,
        "elasticconstant" => Dict(
            "C11" => 684.7,
            "C12" => 373.5,
            "C44" => 311.1,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.20 * Q * Rydberg,
                4 => 0.00 * Q * Rydberg,
                8 => 0.00 * Q * Rydberg,
                11 => 0.04 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.06 * Q * Rydberg,
                4 => 0.05 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.01 * Q * Rydberg,
            )
        ),
    ),
    "AlSb" => Dict(
        "latticeconstant" => 6.1355e-10,
        "elasticconstant" => Dict(
            "C11" => 876.9,
            "C12" => 434.1,
            "C44" => 407.6,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.21 * Q * Rydberg,
                4 => 0 * Q * Rydberg,
                8 => 0.02 * Q * Rydberg,
                11 => 0.06 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.06 * Q * Rydberg,
                4 => 0.04 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.02 * Q * Rydberg,
            )
        ),
    ),
    "GaSb" => Dict(
        "latticeconstant" => 6.0959e-10,
        "elasticconstant" => Dict(
            "C11" => 884.2,
            "C12" => 402.6,
            "C44" => 432.2,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.22 * Q * Rydberg,
                4 => 0.00 * Q * Rydberg,
                8 => 0.00 * Q * Rydberg,
                11 => 0.05 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.06 * Q * Rydberg,
                4 => 0.05 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.01 * Q * Rydberg,
            )
        ),
    ),
    "ZnS" => Dict(
        "latticeconstant" => 0.0,
        "elasticconstant" => Dict(
            "C11" => 10,
            "C12" => 10,
            "C44" => 10,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.22 * Q * Rydberg,
                4 => 0 * Q * Rydberg,
                8 => 0.03 * Q * Rydberg,
                11 => 0.07 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.24 * Q * Rydberg,
                4 => 0.14 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.04 * Q * Rydberg,
            )
        ),
    ),
    "ZnSe" => Dict(
        "latticeconstant" => 0.0,
        "elasticconstant" => Dict(
            "C11" => 10,
            "C12" => 10,
            "C44" => 10,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.23 * Q * Rydberg,
                4 => 0 * Q * Rydberg,
                8 => 0.03 * Q * Rydberg,
                11 => 0.06 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.24 * Q * Rydberg,
                4 => 0.14 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.04 * Q * Rydberg,
            )
        ),
    ),
    "ZnTe" => Dict(
        "latticeconstant" => 0.0,
        "elasticconstant" => Dict(
            "C11" => 10,
            "C12" => 10,
            "C44" => 10,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.22 * Q * Rydberg,
                4 => 0 * Q * Rydberg,
                8 => 0.00 * Q * Rydberg,
                11 => 0.05 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.13 * Q * Rydberg,
                4 => 0.10 * Q * Rydberg,
                8 => 0.00 * Q * Rydberg,
                11 => 0.01 * Q * Rydberg,
            )
        ),
    ),
    "CdTe" => Dict(
        "latticeconstant" => 0.0,
        "elasticconstant" => Dict(
            "C11" => 10,
            "C12" => 10,
            "C44" => 10,
        ),
        "pseudopotential" => Dict(
            "Fs" => Dict(
                3 => -0.20 * Q * Rydberg,
                4 => 0.00 * Q * Rydberg,
                8 => 0.00 * Q * Rydberg,
                11 => 0.04 * Q * Rydberg,
            ),
            "Fa" => Dict(
                3 => 0.15 * Q * Rydberg,
                4 => 0.09 * Q * Rydberg,
                8 => 0 * Q * Rydberg,
                11 => 0.04 * Q * Rydberg,
            )
        ),
    )
)

end


text:
"""
General Interface for Pricing

		Rho=rho(mcProcess,rfCurve,mcBaseData,payoff_)

Where:\n
		mcProcess          = Process to be simulated.
		rfCurve  = Zero Rate Data.
		mcBaseData = Basic properties of MonteCarlo simulation
		payoff_ = Payoff(s) to be priced

		Rho     = Rho of the derivative

"""
function rho(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoff::AbstractPayoff, dr::Real = 1e-7)
    Price = pricer(mcProcess, rfCurve, mcConfig, abstractPayoff)
    rfCurve_1 = ZeroRate(rfCurve.r + dr)
    PriceUp = pricer(mcProcess, rfCurve_1, mcConfig, abstractPayoff)
    rho = (PriceUp - Price) / dr

    return rho
end

function rho(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoffs::Array{abstractPayoff_}, dr::Real = 1e-7) where {abstractPayoff_ <: AbstractPayoff}
    Prices = pricer(mcProcess, rfCurve, mcConfig, abstractPayoffs)
    rfCurve_1 = ZeroRate(rfCurve.r + dr)
    PricesUp = pricer(mcProcess, rfCurve_1, mcConfig, abstractPayoffs)
    rho = (PricesUp .- Prices) ./ dr

    return rho
end

export rho;


text:
using DiffOpt
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

using JuMP
using LinearAlgebra
import Ipopt
import OSQP
import Clp
import SCS

using DelimitedFiles
import MatrixOptInterface
const MatOI = MatrixOptInterface

@testset "Testing forward on trivial QP" begin
    # using example on https://osqp.org/docs/examples/setup-and-solve.html
    Q = [
        4.0 1.0
        1.0 2.0
    ]
    q = [1.0, 1.0]
    G = [
         1.0 1.0
         1.0 0.0
         0.0 1.0
        -1.0 -1.0
        -1.0 0.0
         0.0 -1.0
    ]
    h = [1, 0.7, 0.7, -1, 0, 0]

    model = JuMP.Model(() -> diff_optimizer(Ipopt.Optimizer))
    MOI.set(model, MOI.Silent(), true)
    @variable(model, x[1:2])
    @objective(model, Min, dot(Q * x, x) + dot(q, x))
    @constraint(model,
        G * x .<= h,
    )
    optimize!(model)

    @test JuMP.value.(x) ≈ [0.3, 0.7] atol=ATOL rtol=RTOL
end

@testset "Differentiating trivial QP 1" begin
    Q = [
        4.0 1.0
        1.0 2.0
    ]
    q = [1.0, 1.0]
    G = [1.0 1.0]
    h = [-1.0]

    model = JuMP.direct_model(diff_optimizer(Ipopt.Optimizer))
    MOI.set(model, MOI.Silent(), true)
    x = @variable(model, [1:2])
    @objective(model, Min, dot(Q * x, x) + dot(q, x))
    @constraint(model,
        G * x .<= h,
    )
    optimize!(model)

    @test JuMP.value.(x) ≈ [-0.25, -0.75] atol=ATOL rtol=RTOL

    grad_wrt_h = backward(JuMP.backend(model), ["h"], ones(2))[1]

    @test grad_wrt_h ≈ [1.0] atol=2ATOL rtol=RTOL
end

@testset "Differentiating QP with inequality and equality constraints" begin
    # refered from: https://www.mathworks.com/help/optim/ug/quadprog.html#d120e113424
    # Find equivalent qpth program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_4_py.ipynb

    Q = [
         1.0 -1.0 1.0;
        -1.0  2.0 -2.0;
         1.0 -2.0 4.0
    ]
    q = [2.0, -3.0, 1.0]
    G = [0.0 0.0 1.0
         0.0 1.0 0.0
         1.0 0.0 0.0
         0.0 0.0 -1.0
         0.0 -1.0 0.0
         -1.0 0.0 0.0
    ]
    h = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0,]
    A = [1.0 1.0 1.0;]
    b = [0.5]

    model = Model(() -> diff_optimizer(Ipopt.Optimizer))
    MOI.set(model, MOI.Silent(), true)
    @variable(model, x[1:3])
    @objective(model, Min, dot(Q * x, x) + dot(q, x))
    @constraint(model,
        G * x .<= h,
    )
    @constraint(model,
        A * x .== b,
    )
    optimize!(model)
    doptimizer = JuMP.backend(model).optimizer.model
    z = doptimizer.primal_optimal

    @test z ≈ [0.0, 0.5, 0.0] atol=ATOL rtol=RTOL

    grads = backward(doptimizer, ["Q","q","G","h","A","b"], ones(3))

    dl_dQ = grads[1]
    dl_dq = grads[2]
    dl_dG = grads[3]
    dl_dh = grads[4]
    dl_dA = grads[5]
    dl_db = grads[6]

    @test dl_dQ ≈ zeros(3,3)  atol=ATOL rtol=RTOL

    @test dl_dq ≈ zeros(3,1) atol=ATOL rtol=RTOL

    @test dl_dG ≈ zeros(6,3) atol=ATOL rtol=RTOL

    @test dl_dh ≈ zeros(6,1) atol=ATOL rtol=RTOL

    @test dl_dA ≈ [0.0 -0.5 0.0] atol=ATOL rtol=RTOL

    @test dl_db ≈ [1.0] atol=ATOL rtol=RTOL
end

# refered from https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contquadratic.jl#L3
# Find equivalent CVXPYLayers and QPTH code here:
#               https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_1_py.ipynb
@testset "Differentiating MOI examples 1" begin
    # homogeneous quadratic objective
    # Min x^2 + xy + y^2 + yz + z^2
    # st  x + 2y + 3z >= 4 (c1)
    #     x +  y      >= 1 (c2)
    #     x, y, z \in R

    model = JuMP.direct_model(diff_optimizer(OSQP.Optimizer))
    MOI.set(JuMP.backend(model).optimizer, MOI.RawParameter("eps_prim_inf"), 1e-7)
    MOI.set(JuMP.backend(model).optimizer, MOI.RawParameter("eps_dual_inf"), 1e-7)
    MOI.set(model, MOI.Silent(), true)
    @variables(model, begin
        x
        y
        z
    end)
    @constraint(model, c1,
        x + 2y + 3z >= 4
    )
    @constraint(model, c2,
        x +  y      >= 1
    )
    @objective(model, Min,
        x^2 + x*y + y^2 + y*z + z^2
    )

    optimize!(model)

    z = JuMP.backend(model).primal_optimal

    @test z ≈ [4/7, 3/7, 6/7] atol=ATOL rtol=RTOL

    # obtain gradients
    grads = backward(JuMP.backend(model), ["Q","q","G","h"], ones(3))

    dl_dQ = grads[1]
    dl_dq = grads[2]
    dl_dG = grads[3]
    dl_dh = grads[4]

    @test dl_dQ ≈ [-0.12244895  0.01530609 -0.11224488;
                    0.01530609  0.09183674  0.07653058;
                   -0.11224488  0.07653058 -0.06122449]  atol=ATOL rtol=RTOL

    @test dl_dq ≈ [-0.2142857;  0.21428567; -0.07142857] atol=ATOL rtol=RTOL

    @test_broken dl_dG ≈ [0.05102035   0.30612245  0.255102;
                   0.06122443   0.36734694  0.3061224] atol=ATOL rtol=RTOL

    @test dl_dh ≈ [-0.35714284; -0.4285714] atol=ATOL rtol=RTOL
end

@testset "Differentiating MOI examples 2 - non trivial backward pass vector" begin
    # non-homogeneous quadratic objective
    #    minimize 2 x^2 + y^2 + xy + x + y
    #       s.t.  x, y >= 0
    #             x + y = 1 (c1)

    model = JuMP.Model(() -> diff_optimizer(Ipopt.Optimizer))
    MOI.set(model, MOI.Silent(), true)
    @variable(model, x ≥ 0);
    @variable(model, y ≥ 0);
    @constraint(model, c1, x + y == 1)
    @objective(model, Min, 2 * x^2 + y^2 + x * y + x + y)

    optimize!(model)

    doptimizer = JuMP.backend(model).optimizer.model
    z = doptimizer.primal_optimal

    @test z ≈ [0.25, 0.75] atol=ATOL rtol=RTOL

    # obtain gradients
    dl_dz = [1.3, 0.5]   # choosing a non trivial backward pass vector
    grads = backward(doptimizer, ["Q", "q", "G", "h", "A", "b"], dl_dz)

    dl_dQ = grads[1]
    dl_dq = grads[2]
    dl_dG = grads[3]
    dl_dh = grads[4]
    dl_dA = grads[5]
    dl_db = grads[6]

    @test dl_dQ ≈ [-0.05   -0.05;
                   -0.05    0.15]  atol=ATOL rtol=RTOL

    @test dl_dq ≈ [-0.2; 0.2] atol=ATOL rtol=RTOL

    @test dl_dG ≈ [1e-8  1e-8; 1e-8 1e-8] atol=ATOL rtol=RTOL

    @test dl_dh ≈ [1e-8; 1e-8] atol=ATOL rtol=RTOL

    @test dl_dA ≈ [0.375 -1.075] atol=ATOL rtol=RTOL

    @test dl_db ≈ [0.7] atol=ATOL rtol=RTOL
end

@testset "Differentiating non trivial convex QP JuMP" begin
    nz = 10
    nineq_le = 25
    neq = 10

    # read matrices from files
    names = ["P", "q", "G", "h", "A", "b"]
    matrices = []

    for name in names
        push!(matrices, readdlm(joinpath(dirname(dirname(pathof(DiffOpt))), "test", "data", name * ".txt"), ' ', Float64, '\n'))
    end

    Q, q, G, h, A, b = matrices
    q = vec(q)
    h = vec(h)
    b = vec(b)

    model = JuMP.Model(() -> diff_optimizer(Ipopt.Optimizer))
    MOI.set(model, MOI.Silent(), true)

    @variable(model, x[1:nz])

    @objective(model, Min, dot(Q * x, x) + dot(q, x))
    @constraint(model, G * x .<= h)
    @constraint(model, A * x .== b)

    optimize!(model)

    doptimizer = JuMP.backend(model).optimizer.model

    # obtain gradients
    grads = backward(doptimizer, ["Q", "q", "G", "h", "A", "b"], ones(nz))  # using dl_dz=[1,1,1,1,1,....]

    # read gradients from files
    param_names = ["dP", "dq", "dG", "dh", "dA", "db"]
    grads_actual = []

    for name in param_names
        push!(grads_actual,
            readdlm(
                joinpath(@__DIR__, "data", "$(name).txt"),
                ' ', Float64, '\n',
            )
        )
    end

    grads_actual[2] = vec(grads_actual[2])
    grads_actual[4] = vec(grads_actual[4])
    grads_actual[6] = vec(grads_actual[6])

    # testing differences
    for i in 1:size(grads)[1]
        @test grads[i] ≈  grads_actual[i] atol=1e-2 rtol=1e-2
    end
end

@testset "Differentiating LP; checking gradients for non-active contraints" begin
    # Issue #40 from Gurobi.jl
    # min  x
    # s.t. x >= 0
    #      x >= 3

    model = direct_model(diff_optimizer(Clp.Optimizer))
    MOI.set(model, MOI.Silent(), true)

    @variable(model, x[1:1])

    @objective(model, Min, 1.1 * x[1])
    @constraint(model, x[1] ≥ 0)
    @constraint(model, x[1] ≥ 3)

    optimize!(model)

    # obtain gradients
    grads = backward(JuMP.backend(model), ["G", "h"], [1.0])

    @test grads[1] ≈ [0.0, 3.0] atol=ATOL rtol=RTOL
    @test grads[2] ≈ [0.0, -1.0] atol=ATOL rtol=RTOL

    model = direct_model(diff_optimizer(Clp.Optimizer))
    MOI.set(model, MOI.Silent(), true)

    @variable(model, x[1:1])

    @objective(model, Min, x[1])
    @constraint(model, x[1] ≥ 0)
    @constraint(model, x[1] ≥ 3)

    optimize!(model)

    # obtain gradients
    # fixes https://github.com/jump-dev/DiffOpt.jl/issues/82
    grads = backward(JuMP.backend(model), ["G", "h"], [1.0])
    @test grads[1] ≈ [0.0, 3.0] atol=ATOL rtol=RTOL
    @test grads[2] ≈ [0.0, -1.0] atol=ATOL rtol=RTOL
end

@testset "Differentiating LP; checking gradients for non-active contraints" begin
    # refered from - https://en.wikipedia.org/wiki/Simplex_algorithm#Example

    # max 2x + 3y + 4z
    # s.t. 3x+2y+z <= 10
    #      2x+5y+3z <= 15
    #      x,y,z >= 0

    model = direct_model(diff_optimizer(SCS.Optimizer))
    MOI.set(model, MOI.Silent(), true)
    @variable(model, v[1:3] ≥ 0)
    @objective(model, Min, dot([-2.0, -3.0, -4.0], v))
    (x, y, z) = v

    @constraint(model, 3x+2y+z <= 10)
    @constraint(model, 2x+5y+3z <= 15)
    optimize!(model)

    # obtain gradients
    grads = backward(JuMP.backend(model), ["Q", "q", "G", "h"], ones(3))  # using dl_dz=[1,1,1,1,1,....]

    @test grads[1] ≈ zeros(3,3) atol=ATOL rtol=RTOL
    @test grads[2] ≈ zeros(3) atol=ATOL rtol=RTOL
    @test grads[3] ≈ [0.0 0.0 0.0;
                    0.0 0.0 -5/3;
                    0.0 0.0 5/3;
                    0.0 0.0 -10/3;
                    0.0 0.0 0.0]   atol=ATOL rtol=RTOL
    @test grads[4] ≈ [0.0; 1/3; -1/3; 2/3; 0.0]   atol=ATOL rtol=RTOL
end

@testset "Differentiating LP with variable bounds" begin
    # max 2x + 3y + 4z
    # s.t. 3x+2y+z <= 10
    #      2x+5y+3z <= 15
    #      x ≤ 3
    #      0 ≤ y ≤ 2
    #      z ≥ 2
    #      x,y,z >= 0
    # variant of previous test with same solution

    model = direct_model(diff_optimizer(Clp.Optimizer))
    MOI.set(model, MOI.Silent(), true)
    @variable(model, v[1:3] ≥ 0)
    (x, y, z) = v
    @objective(model, Max, 2x + 3y + 4z)
    @constraint(model, 3x+2y+z <= 10)
    @constraint(model, 2x+5y+3z <= 15)
    @constraint(model, x ≤ 3)
    @constraint(model, 1.0 * y ≤ 2)
    @constraint(model, z ≥ 2)

    optimize!(model)

    # obtain gradients
    grads = backward(JuMP.backend(model), ["Q", "q", "G", "h"], ones(3))  # using dl_dz=[1,1,1]

    @test grads[1] ≈ zeros(3,3) atol=ATOL rtol=RTOL
    @test grads[2] ≈ zeros(3) atol=ATOL rtol=RTOL
    @test grads[3] ≈ [0.0 0.0 0.0
                      0.0 0.0 -5/3
                      0.0 0.0 0.0
                      0.0 0.0 0.0
                      0.0 0.0 0.0
                      0.0 0.0 5/3
                      0.0 0.0 -10/3
                      0.0 0.0 0.0
                      ]   atol=ATOL rtol=RTOL
    @test grads[4] ≈ [0.0, 1/3, 0.0, 0.0, 0.0, -1/3, 2/3, 0.0] atol=ATOL rtol=RTOL

end

@testset "Differentiating simple SOCP" begin
    # referred from _soc2test, https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L1355
    # find equivalent diffcp python program here: https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_socp_1_py.ipynb

    # Problem SOC2
    # min  x
    # s.t. y ≥ 1/√2
    #      x² + y² ≤ 1
    # in conic form:
    # min  x
    # s.t.  -1/√2 + y ∈ R₊
    #        1 - t ∈ {0}
    #      (t,x,y) ∈ SOC₃

    model = JuMP.direct_model(diff_optimizer(SCS.Optimizer))
    MOI.set(model, MOI.Silent(), true)
    @variable(model, x)
    @variable(model, y)
    @variable(model, t)

    @objective(model, Min, x)
    @constraint(model, c1, [-1/√2 + y] in MOI.Nonnegatives(1))
    @constraint(model, c2, [1 - t] in MOI.Zeros(1))
    @constraint(model, c3, [1.0 * t, x, y] in SecondOrderCone())

    optimize!(model)

    v = JuMP.backend(model).primal_optimal
    # slack variables
    s = collect(Iterators.flatten(JuMP.value.([c1, c2, c3])))
    y = collect(Iterators.flatten(JuMP.dual.([c1, c2, c3])))

    # these matrices are benchmarked with the output generated by diffcp
    # refer the python file mentioned above to get equivalent python source code
    @test v ≈ [-1/√2; 1/√2; 1.0] atol=ATOL rtol=RTOL
    @test s ≈ [0.0, 0.0, 1.0, -1/√2, 1/√2] atol=ATOL rtol=RTOL
    @test y ≈ [1, √2, √2, 1, -1] atol=ATOL rtol=RTOL

    dA = zeros(5, 3)
    dA[1:3, :] .= Matrix(1.0I, 3, 3)
    db = zeros(5)
    dc = zeros(3)

    @test_broken backward(JuMP.backend(model), dA, db, dc)

end


text:
import Base: mean

export mean, stdev

immutable Statistics <: AFImpl
    mean::Ptr{Void}
    stdev::Ptr{Void}

    function Statistics(ptr)
        new(
            Libdl.dlsym(ptr, :af_mean),
            Libdl.dlsym(ptr, :af_stdev)
        )
    end
end

@afCall_Arr_Arr_DimT(mean, statistics, mean, -1)

@afCall_Arr_Arr_DimT(stdev, statistics, stdev, -1)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
#=
By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

What is the 10 001st prime number?
=#

using Primes
function findnthprime(n)
  i = 2; foundprimes = 0; lastprime = 0
  while foundprimes < n
    if isprime(i)
      foundprimes += 1
      lastprime = i
    end
    i += 1
  end
  lastprime
end

@time println(findnthprime(10_001))


text:
using Mimi
using Distributions
using CSVFiles
using DataFrames
using CSV
using Random
using StatsBase
using Statistics

# get main_model file
include("../../src/main_model.jl")
include("../../src/mcs.jl")
include("../../src/compute_scc.jl")

if !@isdefined samplesize
    samplesize = 50000
end

for scenario in ["RCP1.9 & SSP1", "RCP2.6 & SSP1", "RCP4.5 & SSP2", "RCP8.5 & SSP5"]
    model = "PAGE-ICE"
    # define model, default settings: getpage(NDCs scenario, permafrost, no sea-ice, no page09damages)
    m = getpage(scenario, true, true)
    # run model
    run(m)

    # open up Explorer UI, for visual exploration of the variables
    # explore(m)

    # do general monte carlo simulation and save the output
    do_monte_carlo_runs(samplesize, scenario, joinpath(@__DIR__, "../../output", scenario, model))

    # get the social cost of carbon for the Monte Carlo simulations and save the output
    sccs = compute_scc_mcs(m, samplesize, year=2020)
    # store results in DataFrame
    df = DataFrame(Any[fill(model, samplesize), fill(scenario, samplesize), sccs], [:ModelName, :ScenarioName, :SCC])
    # write out to csv
    DIR = joinpath(@__DIR__, "..", "..", "output")
    CSV.write(joinpath(DIR, string(model, "_", scenario, "_", "scc.csv")), df)
end


text:
using PkgBenchmark
include("pprinthelper.jl")

if length(ARGS) == 1
    result = PkgBenchmark.readresults(ARGS[1])
elseif length(ARGS) == 0
    res_dir = joinpath(@__DIR__, "results")
    path = maximum(readdir(res_dir))
    path = joinpath(res_dir, path)
    result = PkgBenchmark.readresults(path)
end

displayresult(result)


text:
# Execution time; to send data, do in the TERMINAL  $echo "message" | socat - UNIX-CONNECT:/tmp/tmpfifo

if nworkers() < 2
	newproc = addprocs(1)
	@show workers()
end

# Make module available (it is assumed we are in j4pr root directory)
@everywhere include("src/exp/socket_servers.jl")

# Assign worker, unix socket
pid = workers()[1]			# which worker
sockname = "/tmp/tmpsocket"		# which unix socket
cc=RemoteChannel(()->Channel(1)) 	# notification channel
#srv = socket_servers.start_sync_sock_server
srv = socket_servers.start_async_sock_server

# start server and wait for notification that the worker is no longer running
@async remote_do(srv, pid, sockname, cc)
wait(cc)


text:
"""
myIVBound(y, x)

An implementation of the monotone IV bounds of Manski and Pepper (2000).
"""
struct myIVBound
    px
    Eyx
    unq_x

    function myIVBound(y, x)
        unq_x = sort(unique(x))
        px = zeros(length(unq_x))
        Eyx = zeros(length(unq_x))
        for u in 1:length(unq_x)
            # Estimate P(x = u) and E[y|x = u]
            indx = x .== unq_x[u]
            px[u] = sum(indx) / n
            Eyx[u] = mean(y[indx])
        end
        new(px, Eyx, unq_x)
    end #MYIVBOUND
end #MYIVBOUND

# Methods ======================================================================
"""
coef(fit::myIVBound, s, t)

Function to calculate the upper bound via Monotone Instrumental Variables.
"""
function coef(fit::myIVBound, s, t)
    Eyx_px = fit.Eyx .* fit.px

    from_t = fit.unq_x .> t
    until_s = fit.unq_x .< s

    ub = sum(Eyx_px[from_t]) + (fit.Eyx[fit.unq_x .== t] * sum(fit.px[.!from_t]))[1]
    lb = sum(Eyx_px[until_s]) + (fit.Eyx[fit.unq_x .== s] * sum(fit.px[.!until_s]))[1]

    return ub - lb
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 82ecf9c7-bc7c-4827-a025-920609acfb80
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
    using NativeSVG
	using Plots
	using Ipopt
	using JuMP
end

# ╔═╡ 444a0740-0191-11eb-382f-dd522542f7c4
md"# Interior Point Methods"

# ╔═╡ db7adcc0-0191-11eb-07b4-5b3fbfc45f72
md"""## Introduction

In the 1980s it was discovered that many large linear programs could
be solved efficiently by using formulations and algorithms from nonlinear
programming. One characteristic of these methods was that they required
all iterates to satisfy the inequality constraints in the problem
strictly, so they became known as _interior-point methods_. In
the 1990s, a subclass of interior-point methods known as _primal-dual
methods_ had distinguished themselves as the most efficient practical
approaches, and proved to be strong competitors to the _simplex
method_ on large problems. Recently, it has been shown that
interior-point methods are as successful for nonlinear optimization
as for linear programming. General primal-dual interior point methods
are the focus of this chapter"""

# ╔═╡ e2fe62a0-0191-11eb-0a1d-51d009d5ab5c
md"""## Interior-point Quadratic Programming

For simplicity, we restrict our attention to convex quadratic programs,
which we write as follows:
```math
\begin{aligned}
\min_{\vec{x}}\, & f\left(\vec{x}\right)\overset{\vartriangle}{=}\frac{1}{2} \vec{x}^\mathsf{T}Q\vec{x}- \vec{c}^\mathsf{T}\vec{x}\\
\textrm{subject to}\, & \begin{cases}
A_{\textrm{eq}}\vec{x}=\vec{b}_{\textrm{eq}}\,,\\
A_{\textrm{in}}\vec{x}\leq\vec{b}_{\textrm{in}}\,,
\end{cases}
\end{aligned}
```
where ``Q`` is a symmetric and positive semidefinite ``n\times n`` matrix,
``\vec{c}\in\mathbb{R}^{n}``, ``A_{\textrm{eq}}`` is a ``m\times n`` matrix,
``\vec{b}_{\textrm{eq}}\in\mathbb{R}^{m}``, ``A_{\textrm{in}}`` is a ``p\times n``
matrix and ``\vec{b}_{\textrm{in}}\in\mathbb{R}^{p}``. Rewriting the KKT
conditions in this notation, we obtain
```math
\begin{aligned}
\vec{\mu} & \geq\vec{0}\,,\\
Q\vec{x}+ A_{\textrm{eq}}^\mathsf{T}\vec{\lambda}+ A_{\textrm{in}}^\mathsf{T}\vec{\mu}-\vec{c} & =\vec{0}\,,\\
\left(A_{\textrm{in}}\vec{x}-\vec{b}_{\textrm{in}}\right)_{j}\mu_{j} & =0\,,\quad j=1,\dots,p\,,\\
A_{\textrm{eq}}\vec{x} & =\vec{b}_{\textrm{eq}}\,,\\
A_{\textrm{in}}\vec{x} & \leq\vec{b}_{\textrm{in}}\,.
\end{aligned}
```
By introducing the slack vector ``\vec{y}\geq\vec{0}``, we can rewrite
these conditions as
```math
\begin{aligned}
\left(\vec{\mu},\vec{y}\right) & \geq\vec{0}\,,\\
Q\vec{x}+ A_{\textrm{eq}}^\mathsf{T}\vec{\lambda}+ A_{\textrm{in}}^\mathsf{T}\vec{\mu}-\vec{c} & =\vec{0}\,,\\
y_{j}\mu_{j} & =0\,,\quad j=1,\dots,p\,,\\
A_{\textrm{eq}}\vec{x}-\vec{b}_{\textrm{eq}} & =\vec{0}\,,\\
A_{\textrm{in}}\vec{x}-\vec{b}_{\textrm{in}}+\vec{y} & =\vec{0}\,.
\end{aligned}
```
Since we assume that ``Q`` is positive semidefinite, these KKT conditions
are not only necessary but also sufficient, so we can solve the convex
quadratic program by finding solutions of this system.

Primal-dual methods generate iterates that satisfy the bounds strictly; that is, ``\vec{y}>0``
and ``\vec{\mu}>0``. This property is the origin of the term interior-point.
By respecting these bounds, the methods avoid spurious solutions,
points that satisfy the system but not the bounds. Spurious solutions
abound, and do not provide useful information about real solutions,
so it makes sense to exclude them altogether. Given a current iterate
``\left(\vec{x}^{\left(k\right)},\vec{y}^{\left(k\right)},\vec{\lambda}^{\left(k\right)},\vec{\mu}^{\left(k\right)}\right)``
that satisfies ``\left(\vec{\mu}^{\left(k\right)},\vec{y}^{\left(k\right)}\right)>0``,
we can define a _complementary measure_
```math
\nu_{k}=\frac{ \left(\vec{y}^{\left(k\right)}\right)^\mathsf{T}\vec{\mu}^{\left(k\right)}}{p}\,.
```
This measure gives an indication of the desirability of the couple
``\left(\vec{\mu}^{\left(k\right)},\vec{y}^{\left(k\right)}\right)``.

We derive a path-following, primal-dual method by considering the
perturbed KKT conditions by
```math
\vec{F}\left(\vec{x}^{\left(k+1\right)},\vec{y}^{\left(k+1\right)},\vec{\lambda}^{\left(k+1\right)},\vec{\mu}^{\left(k+1\right)},\sigma_{k}\nu_{k}\right)=\begin{pmatrix}Q\vec{x}^{\left(k+1\right)}+ A_{\textrm{eq}}^\mathsf{T}\vec{\lambda}^{\left(k+1\right)}+ A_{\textrm{in}}^\mathsf{T}\vec{\mu}^{\left(k+1\right)}-\vec{c}\\
Y_{k+1}M_{k+1}\vec{1}-\sigma_{k}\nu_{k}\vec{1}\\
A_{\textrm{eq}}\vec{x}^{\left(k+1\right)}-\vec{b}_{\textrm{eq}}\\
A_{\textrm{in}}\vec{x}^{\left(k+1\right)}-\vec{b}_{\textrm{in}}+\vec{y}^{\left(k+1\right)}
\end{pmatrix}=\vec{0}\,,
```
where
```math
Y_{k+1}=\begin{pmatrix}y_{1}^{\left(k+1\right)} & 0 & \cdots & 0\\
0 & y_{2}^{\left(k+1\right)} & \ddots & 0\\
\vdots & \ddots & \ddots & 0\\
0 & 0 & 0 & y_{p}^{\left(k+1\right)}
\end{pmatrix}\,,\quad M_{k+1}=\begin{pmatrix}\mu_{1}^{\left(k+1\right)} & 0 & \cdots & 0\\
0 & \mu_{2}^{\left(k+1\right)} & \ddots & 0\\
\vdots & \ddots & \ddots & 0\\
0 & 0 & 0 & \mu_{p}^{\left(k+1\right)}
\end{pmatrix}\,,
```
and ``\sigma\in\left[0,1\right]`` is the reduction factor that we wish
to achieve in the complementary measure on one step. We call ``\sigma``
the _centering parameter_. The solution of this system for all
positive values of ``\sigma`` and ``\nu`` define the _central path_,
which is a trajectory that leads to the solution of the quadratic
program as ``\sigma\nu`` tends to zero.

By fixing ``\sigma_{k}`` and applying Newton's method to the system,
we obtain the linear system
```math
\begin{pmatrix}Q & 0 &  A_{\textrm{eq}}^\mathsf{T} &  A_{\textrm{in}}^\mathsf{T}\\
0 & M_{k} & 0 & Y_{k}\\
A_{\textrm{eq}} & 0 & 0 & 0\\
A_{\textrm{in}} & I & 0 & 0
\end{pmatrix}\begin{pmatrix}\vec{d}_{\vec{x}}^{\left(k\right)}\\
\vec{d}_{\vec{y}}^{\left(k\right)}\\
\vec{d}_{\vec{\lambda}}^{\left(k\right)}\\
\vec{d}_{\vec{\mu}}^{\left(k\right)}
\end{pmatrix}=-\begin{pmatrix}Q\vec{x}^{\left(k\right)}+ A_{\textrm{eq}}^\mathsf{T}\vec{\lambda}^{\left(k\right)}+ A_{\textrm{in}}^\mathsf{T}\vec{\mu}^{\left(k\right)}-\vec{c}\\
Y_{k}M_{k}\vec{1}-\sigma_{k}\nu_{k}\vec{1}\\
A_{\textrm{eq}}\vec{x}^{\left(k\right)}-\vec{b}_{\textrm{eq}}\\
A_{\textrm{in}}\vec{x}^{\left(k\right)}-\vec{b}_{\textrm{in}}+\vec{y}^{\left(k\right)}
\end{pmatrix}\,.
```
We obtain the next iterate by setting
```math
\begin{pmatrix}\vec{x}^{\left(k+1\right)}\\
\vec{y}^{\left(k+1\right)}\\
\vec{\lambda}^{\left(k+1\right)}\\
\vec{\mu}^{\left(k+1\right)}
\end{pmatrix}=\begin{pmatrix}\vec{x}^{\left(k\right)}\\
\vec{y}^{\left(k\right)}\\
\vec{\lambda}^{\left(k\right)}\\
\vec{\mu}^{\left(k\right)}
\end{pmatrix}+\alpha_{k}\begin{pmatrix}\vec{d}_{\vec{x}}^{\left(k\right)}\\
\vec{d}_{\vec{y}}^{\left(k\right)}\\
\vec{d}_{\vec{\lambda}}^{\left(k\right)}\\
\vec{d}_{\vec{\mu}}^{\left(k\right)}
\end{pmatrix}\,,
```
where ``\alpha_{k}`` is chosen to retain the bounds ``\left(\vec{\mu}^{\left(k+1\right)},\vec{y}^{\left(k+1\right)}\right)>0``
and possibly to satisfy various other conditions.

The choices of centering parameter ``\sigma_{k}`` and step-length ``\alpha_{k}``
are crucial for the performance of the method. Techniques for controlling
these parameters, directly and indirectly, give rise to a wide variety
of methods with diverse properties. One option is to use equal step
length for the primal and dual updates, and to set ``\alpha_{k}=\min\left\{ \alpha_{k}^{\textrm{pri}},\alpha_{k}^{\textrm{dual}}\right\} ``,
where
```math
\begin{aligned}
\alpha_{k}^{\textrm{pri}} & =\max\left\{ \alpha\in\left\{ 0,1\right\} :\vec{y}^{\left(k\right)}+\alpha\vec{d}_{\vec{y}}^{\left(k\right)}\geq\left(1-\tau\right)\vec{y}^{\left(k\right)}\right\} \,,\\
\alpha_{k}^{\textrm{dual}} & =\max\left\{ \alpha\in\left\{ 0,1\right\} :\vec{\mu}^{\left(k\right)}+\alpha\vec{d}_{\vec{\mu}}^{\left(k\right)}\geq\left(1-\tau\right)\vec{\mu}^{\left(k\right)}\right\} \,,
\end{aligned}
```
the parameter ``\tau\in\left]0,1\right[`` controls how far we back
off from the maximum step for which the conditions ``\vec{y}^{\left(k\right)}+\alpha\vec{d}_{\vec{y}}^{\left(k\right)}\geq\vec{0}``
and ``\vec{\mu}^{\left(k\right)}+\alpha\vec{d}_{\vec{\mu}}^{\left(k\right)}\geq\vec{0}``
are satisfied. A typical value of ``\tau=0.995`` and we can choose
``\tau_{k}`` to approach ``1`` as the iterates approach the solution,
to accelerate the convergence.

The most popular interior-point method for convex QP is based on Mehrotra's
predictor-corrector. First we compute an affine scaling step ``\left(\vec{d}_{\vec{x},\textrm{aff}},\vec{d}_{\vec{y},\textrm{aff}},\vec{d}_{\vec{\lambda},\textrm{aff}},\vec{d}_{\vec{\mu},\textrm{aff}}\right)``
by setting ``\sigma_{k}=0``. We improve upon this step by computing
a corrector step. Next, we compute the centering parameter ``\sigma_{k}``
using following heuristic
```math
\sigma_{k}=\left(\frac{\nu_{\textrm{aff}}}{\nu_{k}}\right)^{3}\,,
```
where ``\nu_{\textrm{aff}}=\frac{ \left(\vec{y}_{\textrm{aff}}\right)^\mathsf{T}\left(\vec{\mu}_{\textrm{aff}}\right)}{p}``.
The total step is obtained by solving the following system
```math
\begin{pmatrix}Q & 0 &  A_{\textrm{eq}}^\mathsf{T} &  A_{\textrm{in}}^\mathsf{T}\\
0 & M_{k} & 0 & Y_{k}\\
A_{\textrm{eq}} & 0 & 0 & 0\\
A_{\textrm{in}} & I & 0 & 0
\end{pmatrix}\begin{pmatrix}\vec{d}_{\vec{x}}^{\left(k\right)}\\
\vec{d}_{\vec{y}}^{\left(k\right)}\\
\vec{d}_{\vec{\lambda}}^{\left(k\right)}\\
\vec{d}_{\vec{\mu}}^{\left(k\right)}
\end{pmatrix}=-\begin{pmatrix}Q\vec{x}^{\left(k\right)}+ A_{\textrm{eq}}^\mathsf{T}\vec{\lambda}^{\left(k\right)}+ A_{\textrm{in}}^\mathsf{T}\vec{\mu}^{\left(k\right)}-\vec{c}\\
Y_{k}M_{k}\vec{1}+\Delta Y_{\textrm{aff}}\Delta M_{\textrm{aff}}\vec{1}-\sigma_{k}\nu_{k}\vec{1}\\
A_{\textrm{eq}}\vec{x}^{\left(k\right)}-\vec{b}_{\textrm{eq}}\\
A_{\textrm{in}}\vec{x}^{\left(k\right)}-\vec{b}_{\textrm{in}}+\vec{y}^{\left(k\right)}
\end{pmatrix}\,,
```
where
```math
\Delta Y_{\textrm{aff}}=Y_{\textrm{aff}}-Y_{k}\,,\quad\Delta M_{\textrm{aff}}=M_{\textrm{aff}}-M_{k}\,.
```"""

# ╔═╡ 578662d0-0192-11eb-027a-7d04d8c5524d
md"""## Julia

The package `Ipopt` integrated in `JuMP` solves a quadratic problem with an interior point method.
"""

# ╔═╡ 705fed7e-0192-11eb-3f1f-9d13ebbb4fb9
let
	model = Model(with_optimizer(Ipopt.Optimizer))
	@variable(model, 0 <= x, start=2)
	@variable(model, 0 <= y, start=0)
	@NLobjective(model, Min, (x-1)^2+(y-2.5)^2)
	@constraint(model, con1, -x+2y <= 2)
	@constraint(model, con2,  x+2y <= 6)
	@constraint(model, con3,  x-2y <= 2)
	optimize!(model)
    value(x), value(y)
end

# ╔═╡ 947daf90-0192-11eb-38f5-d921e71228e3
md"""## General interior-point method

We consider the general constrained problem
```math
\begin{aligned}
\min\, & f\left(\vec{x}\right)\\
\textrm{subject to} & \begin{cases}
\vec{h}\left(\vec{x}\right)=\vec{0}\\
\vec{g}\left(\vec{x}\right)\leq\vec{0}
\end{cases}
\end{aligned}
```
where $f: \mathbb{R}^{n}\rightarrow \mathbb{R}$, $\vec{h}: \mathbb{R}^{n}\rightarrow \mathbb{R}^{m}$,
$m\leq n$, and $\vec{g}: \mathbb{R}^{n}\rightarrow \mathbb{R}^{p}$. The general
interior-point method can be viewed as a direct extension of interior-point
methods for quadratic programming.

The KKT conditions for the nonlinear optimization problem can be written
as
```math
\begin{aligned}
\left(\vec{\mu},\vec{y}\right) & \geq\vec{0}\,,\\
\nabla f\left(\vec{x}\right)+ \mathsf{J}^\mathsf{T}\vec{h}\left(\vec{x}\right)\vec{\lambda}+ \mathsf{J}^\mathsf{T}\vec{g}\left(\vec{x}\right)\vec{\mu} & =\vec{0}\,,\\
YM\vec{1} -\sigma\nu\vec{1}&=0\,,\\
\vec{h}\left(\vec{x}\right) & =\vec{0}\,,\\
\vec{g}\left(\vec{x}\right)+\vec{y} & =\vec{0}\,,
\end{aligned}
```
with $\nu=0$. This first and third equation from the system introduce
into the problem the combinatorial aspects of determining the optimal
active set. We circumvent this difficulty be letting $\nu$ be strict
positive, thus forcing the variables $\vec{\mu}$ and $\vec{y}$ to
take positive values. We solve approximately the perturbed KKT conditions
for a sequence of positive parameters $\left\{ \sigma_{k}\nu_{k}\right\} $
that converges to zero, while maintaining $\left(\vec{\mu}^{\left(k\right)},\vec{y}^{\left(k\right)}\right)\geq\vec{0}$.

Applying Newton's method to the nonlinear system, we obtain
```math
\begin{gathered}
\begin{bmatrix}\mathsf{H} \mathcal{L}\left(\vec{x}^{\left(k\right)},\vec{y}^{\left(k\right)},\vec{\lambda}^{\left(k\right)},\vec{\mu}^{\left(k\right)}\right) & 0 &  \mathsf{J}^\mathsf{T}\vec{h}\left(\vec{x}^{\left(k\right)}\right) &  \mathsf{J}^\mathsf{T}\vec{g}\left(\vec{x}^{\left(k\right)}\right)\\
0 & M_{k} & 0 & Y_{k}\\
\mathsf{J}\vec{h}\left(\vec{x}^{\left(k\right)}\right) & 0 & 0 & 0\\
\mathsf{J}\vec{g}\left(\vec{x}^{\left(k\right)}\right) & I & 0 & 0
\end{bmatrix}\begin{bmatrix}\vec{d}_{\vec{x}}^{\left(k\right)}\\
\vec{d}_{\vec{y}}^{\left(k\right)}\\
\vec{d}_{\vec{\lambda}}^{\left(k\right)}\\
\vec{d}_{\vec{\mu}}^{\left(k\right)}
\end{bmatrix}=\\\qquad\qquad\qquad\qquad-\begin{bmatrix}\nabla f\left(\vec{x}^{\left(k\right)}\right)+ \mathsf{J}^\mathsf{T}\vec{h}\left(\vec{x}^{\left(k\right)}\right)\vec{\lambda}^{\left(k\right)}+ \mathsf{J}^\mathsf{T}\vec{g}\left(\vec{x}^{\left(k\right)}\right)\vec{\mu}^{\left(k\right)}\\
Y_{k}M_{k}-\sigma_{k}\nu_{k}\vec{1}\\
\vec{h}\left(\vec{x}^{\left(k\right)}\right)\\
\vec{g}\left(\vec{x}^{\left(k\right)}\right)+\vec{y}^{\left(k\right)}
\end{bmatrix}\,.
\end{gathered}
```
Numerical experience has shown, however, that using different step-lengths
in the primal and dual variables often leads to better convergence.
After the step $k$ has been computed and the maximum step lengths have been
determined, we perform a backtracking line search that computes the
step-lengths.
```math
\alpha_{k,\textrm{back}}^{\textrm{pri}}\in\left]0,\alpha_{k}^{\textrm{pri}}\right]\,,\quad\alpha_{k,\textrm{back}}^{\textrm{dual}}\in\left]0,\alpha_{k}^{\textrm{dual}}\right]\,,
```
providing sufficient decrease.

The centering parameter $\sigma_{k}$ can be calculated by a predictor-corrector
strategy as in the quadratic programming case. """

# ╔═╡ 96b90420-0193-11eb-3f8f-b1fdc9e66341
let
	m = Model(with_optimizer(Ipopt.Optimizer))

	a1 = 2
	b1 = 0
	a2 = -1
	b2 = 1

	@variable(m, x1, start=1.234)
	@variable(m, x2 >= 0, start=5.678)

	@NLobjective(m, Min, sqrt(x2))
	@NLconstraint(m, x2 >= (a1*x1+b1)^3)
	@NLconstraint(m, x2 >= (a2*x1+b2)^3)

	optimize!(m)
	value(x1), value(x2)
end

# ╔═╡ Cell order:
# ╟─82ecf9c7-bc7c-4827-a025-920609acfb80
# ╟─444a0740-0191-11eb-382f-dd522542f7c4
# ╟─db7adcc0-0191-11eb-07b4-5b3fbfc45f72
# ╟─e2fe62a0-0191-11eb-0a1d-51d009d5ab5c
# ╟─578662d0-0192-11eb-027a-7d04d8c5524d
# ╠═705fed7e-0192-11eb-3f1f-9d13ebbb4fb9
# ╟─947daf90-0192-11eb-38f5-d921e71228e3
# ╠═96b90420-0193-11eb-3f8f-b1fdc9e66341


text:
###
#   Simple implmentation of the Metropolis Hasting algorithm in julia
#
#            Institute for Risk and Uncertainty, Uni of Liverpool
#
#                       Authors: Ander Gray, Adolphus Lye
#
#                       Email: Ander.Gray@liverpool.ac.uk,
#                              Adolphus.Lye@liverpool.ac.uk
#
#
#
#       W. K. Hastings (1970). Monte Carlo sampling methods using
#       Markov chains and their applications. Biometrika,
#       57(1), 97-109. doi:10.1093/biomet/57.1.97
#
###

###
#   Add description of function and inputs
###

function MHsample(Target :: Function, Prop, start :: Vector{<:Real}, Nsamples :: Integer, burnin :: Integer = 50, thin ::Integer = 3; islogged :: Bool = true)

    dims = length(start)                                    # Dimensions of input/ prior

    PropRnd = mu -> rand(Prop(mu))                          # Generates a sample from the proposal given mean mu

    chain = zeros( Nsamples * thin + burnin, dims)
    chain[1,:] = start
    accRate = 0

    islogged ? evalDen =  x -> exp(Target(x)) : evalDen = x -> Target(x)

    for i = 2:( Nsamples*thin +burnin)

        next = PropRnd(chain[i-1,:])              # Draw candidate

        targDen = evalDen(next)                  # Target Density at next sample
        targPrevious = evalDen(chain[i-1,:])     # Target Density at current sample

        propDen = pdf(Prop(chain[i-1,:]), next)         # Proposal at next centred at current
        propPrevious = pdf(Prop(next), chain[i-1,:])    # Propsoal at current centred at next

        α = targDen/targPrevious * propPrevious/propDen     # General formula for acceptance probability

        accepted = α >= rand()

        if accepted
            chain[i,:] = next
            accRate = accRate +1
        else
            chain[i,:] = chain[i-1, :]
        end

    end
    accRate = accRate/(Nsamples*thin+burnin)
    return chain[burnin+1:thin:end,:], accRate
end

##
#   Symetric proposal and logged target
##
function MHsampleSimple(Target :: Function, PropRnd, start :: Vector{<:Real}, Nsamples :: Integer, burnin :: Integer = 50, thin ::Integer = 3)

    dims = length(start)                                    # Dimensions of input/ prior
    chain = zeros( Nsamples * thin + burnin, dims)
    chain[1,:] = start[:]
    accRate = 0

    for i = 2:( Nsamples*thin+burnin)

        next = PropRnd(chain[i-1,:])                 # Draw candidate

        targDen = Target(next)[1]                    # Target Density at next sample
        targPrevious = Target(chain[i-1,:])[1]       # Target Density at current sample

        α =  min(0, targDen - targPrevious)

        accepted = α >= log(rand())

        if accepted
            chain[i,:] = next
            accRate = accRate +1
        else
            chain[i,:] = chain[i-1, :]
        end

    end
    accRate = accRate/(Nsamples*thin+burnin)
    return chain[burnin+1:thin:end,:], accRate
end


text:
using Test
using PackedReadWrite

struct Example
    float::Float32
    byte::UInt8
    int::Int32
end

@reflect Example

struct ExampleNoReflect
    float::Float32
    byte::UInt8
    int::Int32
end

struct ExamplePrivate
    ExamplePrivate() = new(1, 2, 3)
    float::Float32
    byte::UInt8
    int::Int32
end

@reflect ExamplePrivate

module SomeOtherModule
    struct OtherStruct
        byte1::UInt8
        byte2::UInt8
    end
end

@reflect SomeOtherModule.OtherStruct

@testset "Test read." begin

    example = Example(10,20,30)
    io = IOBuffer()
    write(io, example.float)
    write(io, example.byte)
    write(io, example.int)
    seekstart(io)
    example_read = read_struct(io, Example)
    @test example_read == example

    example = ExampleNoReflect(10,20,30)
    io = IOBuffer()
    write(io, example.float)
    write(io, example.byte)
    write(io, example.int)
    seekstart(io)
    example_read = read_struct(io, ExampleNoReflect)
    @test example_read == example

    io = IOBuffer()
    other = SomeOtherModule.OtherStruct(1,2)
    write(io, other.byte1)
    write(io, other.byte2)
    seekstart(io)
    other_read = read_struct(io, SomeOtherModule.OtherStruct)
    @test other_read == other

end

@testset "Test read errors." begin

    io = IOBuffer(rand(UInt8, 100))

    @test_throws Exception read_struct(io, ExamplePrivate)
    @test_throws Exception read_struct(io, Int64)

end

@testset "Test write." begin

    example = Example(10,20,30)
    io = IOBuffer()
    bytes_written = write_struct(io, example)
    @test bytes_written == 9
    @test bytes_written == position(io)
    seekstart(io)
    @test read(io, Float32) == 10
    @test read(io, UInt8) == 20
    @test read(io, Int32) == 30

    example = ExampleNoReflect(10,20,30)
    io = IOBuffer()
    bytes_written = write_struct(io, example)
    @test bytes_written == 9
    @test bytes_written == position(io)
    seekstart(io)
    @test read(io, Float32) == 10
    @test read(io, UInt8) == 20
    @test read(io, Int32) == 30

    io = IOBuffer()
    example_private = ExamplePrivate()
    write_struct(io, example_private)
    seekstart(io)
    @test read(io, Float32) == 1
    @test read(io, UInt8) == 2
    @test read(io, Int32) == 3

    other = SomeOtherModule.OtherStruct(1,2)
    io = IOBuffer()
    @test write_struct(io, other) == 2
    seekstart(io)
    @test read(io, UInt8) == 1
    @test read(io, UInt8) == 2

end

@testset "Test write error." begin
    io = IOBuffer()
    @test_throws Exception write_struct(io, 42)
end

struct ExampleParam{T}
    item::T
end

@reflect ExampleParam{UInt16}

@testset "Test ExampleParam" begin

    io = IOBuffer()
    @test write_struct(io, ExampleParam{UInt16}(42)) == 2
    seekstart(io)
    @test read_struct(io, ExampleParam{UInt16}).item == 42

end


text:
using OrdinaryDiffEq
using JLD2
#using Plots

function f(x,θ,t)
    α,β,γ,δ = θ
    x1,x2=x
    dx1 = α*x1 - β*x1*x2
    dx2 = δ*x1*x2 - γ*x2
    [dx1,dx2]
end

θ0 = [0.1,0.2,0.3,0.2]
x0 = [1.0,1.0]
tspan = (0., 100.)
prob = ODEProblem(f,x0,tspan,θ0)
sol = solve(prob,Tsit5())
#p = plot(sol)

dt = 0.2
ts = (tspan[1]):dt:(tspan[2])
us = reduce(hcat, sol(ts).u)

data = Dict(:u=>us, :t=>collect(ts))
jldsave("lotkadata.jld2", u=us, t=collect(ts))


text:
#=

  Euler #14 in Julia.

  Problem 14
  """
  The following iterative sequence is defined for the set of positive integers:

  n n/2 (n is even)
  n 3n + 1 (n is odd)

  Using the rule above and starting with 13, we generate the following
  sequence:
  13 40 20 10 5 16 8 4 2 1

  It can be seen that this sequence (starting at 13 and finishing at 1)
  contains
  10 terms. Although it has not been proved yet (Collatz Problem), it is
  thought that all starting numbers finish at 1.

  Which starting number, under one million, produces the longest chain?

  NOTE: Once the chain starts the terms are allowed to go above one million.)
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using Memoization
include("Euler.jl")

# It's much slower with memoizing!
@memoize function hailstone_memo(n)
    if n % 2 == 0
        return round(Int,n / 2)
    else
        return 3*n+1;
    end
end

function hailstone(n)
    if n % 2 == 0
        return round(Int,n / 2)
    else
        return 3*n+1;
    end
end

#
# Using 1..1_000_000: 0.68331034s
# Using 3..2..1_000_000: 0.36079831s
function euler14a()
    maxN = 0
    maxLen = 0
    for n in 1:1_000_000
    # for n in 3:2:1_000_000
        m = n
        alen = 0
        while m > 1
            m = hailstone(m)
            alen += 1
        end
        alen += 1
        if alen > maxLen
            maxN = n
            maxLen = alen
        end
    end
    return maxN
end

# Cache the lengths.
# 1..999_999: 1.66648847s
# 3..2..999_999: 1.09878097s
function euler14b()
    hash = Dict()
    maxN = 0
    maxLen = 0
    for n in 2:1_000_000
    # for n in 3:2:1_000_000
        m = n
        mlen = 1
        while m > 1
            if haskey(hash,m)
                mlen = hash[m]+mlen-1
                m = 1
            else
                m = hailstone(m)
                mlen += 1
            end
        end
        if !haskey(hash,n)
            hash[n] = mlen
        end
        if mlen > maxLen
            maxN = n
            maxLen = mlen
        end
    end
    return maxN;
end

# Using an array instead:
# 0.04143852s
function euler14c()
    limit = 1_000_000
    hash = zeros(limit)
    maxN = 0
    maxLen = 0
    for n = 2:1_000_000
        m = n
        mlen = 1
        while m > 1
            if m <= limit && hash[m] > 0
                mlen = hash[m]+mlen-1
                m = 1
            else
                m = hailstone(m)
                mlen += 1
            end
        end
        if hash[n] == 0
            hash[n] = mlen
        end
        if mlen > maxLen
            maxN = n
            maxLen = mlen
        end
    end
    return maxN;
end

# run_euler(euler14a)
# run_euler(euler14b)
run_euler(euler14c)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module Blocks
  ##############################################################################
  # Exports
  ##############################################################################

  export Block
  export ==, tostring

  ##############################################################################
  # Implementation
  ##############################################################################

  """
  Represents a unique block of data within a file.
  """
  type Block
    filename::AbstractString
    number::Int
  end

  function ==(a::Block, b::Block)
    a.filename == b.filename && a.number == b.number
  end

  function !=(a::Block, b::Block)
    !(a == b)
  end

  function Base.hash(block::Block)
    hash(tostring(block))
  end

  # Pretty print
  function Base.show(io::IO, block::Block)
    tostring(block)
  end

  function tostring(block::Block)
    "Block[filename: $(block.filename), number: $(block.number)]"
  end
end


text:
__precompile__()

module GradOptim

using Compat
using Compat.LinearAlgebra.BLAS

export Adam, update!, iRProp

"""
    Adam(;lr=0.001, gclip=0, beta1=0.9, beta2=0.999, eps=1e-8)
    update!(w,g,p::Adam)

Container for parameters of the Adam optimization algorithm used by
[`update!`](@ref).

Adam is one of the methods that compute the adaptive learning rate. It
stores accumulated gradients (first moment) and the sum of the squared
of gradients (second).  It scales the first and second moment as a
function of time. Here is the update formulas:

    m = beta1 * m + (1 - beta1) * g
    v = beta2 * v + (1 - beta2) * g .* g
    mhat = m ./ (1 - beta1 ^ t)
    vhat = v ./ (1 - beta2 ^ t)
    w = w - (lr / (sqrt(vhat) + eps)) * mhat

where `w` is the weight, `g` is the gradient of the objective function
w.r.t `w`, `lr` is the learning rate, `m` is an array with the same
size and type of `w` and holds the accumulated gradients. `v` is an
array with the same size and type of `w` and holds the sum of the
squares of the gradients. `eps` is a small constant to prevent a zero
denominator. `beta1` and `beta2` are the parameters to calculate bias
corrected first and second moments. `t` is the update count.

If `vecnorm(g) > gclip > 0`, `g` is scaled so that its norm is equal
to `gclip`.  If `gclip==0` no scaling takes place.

Reference: [Kingma, D. P., & Ba,
J. L. (2015)](https://arxiv.org/abs/1412.6980). Adam: a Method for
Stochastic Optimization. International Conference on Learning
Representations, 1–13.

"""
mutable struct Adam
    lr::AbstractFloat
    gclip::AbstractFloat
    beta1::AbstractFloat
    beta2::AbstractFloat
    eps::AbstractFloat
    t::Int
    fstm
    scndm
end

Adam(; lr=0.001, gclip=0, beta1=0.9, beta2=0.999, eps=1e-8)=Adam(lr, gclip, beta1, beta2, eps, 0, nothing, nothing)

function update!(w, g, p::Adam)
    gclip!(g, p.gclip)
    if p.fstm===nothing; p.fstm=zeros(w); p.scndm=zeros(w); end
    p.t += 1
    BLAS.scale!(p.beta1, p.fstm)
    BLAS.axpy!(1-p.beta1, g, p.fstm)
    BLAS.scale!(p.beta2, p.scndm)
    BLAS.axpy!(1-p.beta2, g .* g, p.scndm)
    fstm_corrected = p.fstm / (1 - p.beta1 ^ p.t)
    scndm_corrected = p.scndm / (1 - p.beta2 ^ p.t)
    BLAS.axpy!(-p.lr, (fstm_corrected ./ (sqrt.(scndm_corrected) + p.eps)), w)
end

function gclip!(g, gclip)
    if gclip == 0
        g
    else
        gnorm = vecnorm(g)
        if gnorm <= gclip
            g
        else
            BLAS.scale!(gclip/gnorm, g)
        end
    end
end

mutable struct iRProp{T}
    eta_pos::T
    eta_neg::T
    dmax::T
    dmin::T
    d0::T
    gpre::Vector{T}
    Δ::Vector{T}
    iRProp(eta_pos::T, eta_neg::T; d0=1.5e-3π, dmin=1e-6π, dmax=6e-3π) where T = new{T}(eta_pos, eta_neg, dmax, dmin, d0)
end

function update!(w, g, p::iRProp{T}) where T
    N = length(g)
    if !isdefined(p, :gpre)
        p.gpre = zeros(T, N)
        p.Δ = p.d0*ones(T, N)
    end
    sgn = sign.(g.*p.gpre)
    posmask = sgn .> 0
    negmask = sgn .< 0
    p.Δ[posmask] = min.(p.eta_pos.*p.Δ[posmask], p.dmax)
    p.Δ[negmask] = max.(p.eta_neg.*p.Δ[negmask], p.dmin)

    p.gpre[:] = g
    p.gpre[negmask] = 0

    w[:] .-= sign.(p.gpre).*p.Δ
end

end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
##################################################################
#
#   Maximal cliques of undirected graph
#   Derived from Graphs.jl: https://github.com/julialang/Graphs.jl
#
##################################################################

"""
    maximal_cliques(g)

Return a vector of vectors representing the node indices in each of the maximal
cliques found in the undirected graph `g`.

```jldoctest
julia> using LightGraphs
julia> g = Graph(3)
julia> add_edge!(g, 1, 2)
julia> add_edge!(g, 2, 3)
julia> maximal_cliques(g)
2-element Array{Array{Int64,N},1}:
 [2,3]
 [2,1]
```
"""
function maximal_cliques end
@traitfn function maximal_cliques(g::::(!IsDirected))
    T = eltype(g)
    # Cache nbrs and find first pivot (highest degree)
    maxconn = -1
    nnbrs = Vector{Set{T}}()
    for n in vertices(g)
        push!(nnbrs, Set{T}())
    end
    pivotnbrs = Set{T}() # handle empty graph
    pivotdonenbrs = Set{T}()  # initialize

    for n in vertices(g)
        nbrs = Set{T}()
        union!(nbrs, out_neighbors(g, n))
        delete!(nbrs, n) # ignore edges between n and itself
        conn = length(nbrs)
        if conn > maxconn
            pivotnbrs = nbrs
            nnbrs[n] = pivotnbrs
            maxconn = conn
        else
            nnbrs[n] = nbrs
        end
    end

    # Initial setup
    cand = Set{T}(vertices(g))
    # union!(cand, keys(nnbrs))
    smallcand = setdiff(cand, pivotnbrs)
    done = Set{T}()
    stack = Vector{Tuple{Set{T}, Set{T}, Set{T}}}()
    clique_so_far = Vector{T}()
    cliques = Vector{Array{T}}()

    # Start main loop
    while !isempty(smallcand) || !isempty(stack)
        if !isempty(smallcand) # Any vertices left to check?
            n = pop!(smallcand)
        else
            # back out clique_so_far
            cand, done, smallcand = pop!(stack)
            pop!(clique_so_far)
            continue
        end
        # Add next node to clique
        push!(clique_so_far, n)
        delete!(cand, n)
        push!(done, n)
        nn = nnbrs[n]
        new_cand = intersect(cand, nn)
        new_done = intersect(done, nn)
        # check if we have more to search
        if isempty(new_cand)
            if isempty(new_done)
                # Found a clique!
                push!(cliques, collect(clique_so_far))
            end
            pop!(clique_so_far)
            continue
        end
        # Shortcut--only one node left!
        if isempty(new_done) && length(new_cand) == 1
            push!(cliques, cat(1, clique_so_far, collect(new_cand)))
            pop!(clique_so_far)
            continue
        end
        # find pivot node (max connected in cand)
        # look in done vertices first
        numb_cand = length(new_cand)
        maxconndone = -1
        for n in new_done
            cn = intersect(new_cand, nnbrs[n])
            conn = length(cn)
            if conn > maxconndone
                pivotdonenbrs = cn
                maxconndone = conn
                if maxconndone == numb_cand
                    break
                end
            end
        end
        # Shortcut--this part of tree already searched
        if maxconndone == numb_cand
            pop!(clique_so_far)
            continue
        end
        # still finding pivot node
        # look in cand vertices second
        maxconn = -1
        for n in new_cand
            cn = intersect(new_cand, nnbrs[n])
            conn = length(cn)
            if conn > maxconn
                pivotnbrs = cn
                maxconn = conn
                if maxconn == numb_cand - 1
                    break
                end
            end
        end
        # pivot node is max connected in cand from done or cand
        if maxconndone > maxconn
            pivotnbrs = pivotdonenbrs
        end
        # save search status for later backout
        push!(stack, (cand, done, smallcand))
        cand = new_cand
        done = new_done
        smallcand = setdiff(cand, pivotnbrs)
    end
    return cliques
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using Plots, SpectralDistances, DSP
m = examplemodels(2)

plot(m[1])
assignmentplot(m...)
flowplot(OptimalTransportRootDistance(domain=Continuous()), m...)

y = randn(10000)
plot(spectrogram(y))
plot(periodogram(y))
plot(randn(ComplexF64, 10), Continuous())
plot(ContinuousRoots(randn(ComplexF64, 3)))


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(setindex!),RegCol{1},FNum{Int64},TableIndex{3}})   # time: 0.0012955
end


text:
abstract type TemperatureProfile{FT} end

struct DecayingTemperatureProfile{FT} <: TemperatureProfile{FT}
    "Virtual temperature at surface (K)"
    T_virt_surf::FT
    "Minimum virtual temperature at the top of the atmosphere (K)"
    T_min_ref::FT
    "Height scale over which virtual temperature drops (m)"
    H_t::FT
    function DecayingTemperatureProfile{FT}(
        params::NamedTuple,
        _T_virt_surf::FT,
        _T_min_ref::FT,
        H_t::FT = FT(params.R_d) * _T_virt_surf / FT(params.g),
    ) where {FT}
        return new{FT}(_T_virt_surf, _T_min_ref, H_t)
    end
end

function (profile::DecayingTemperatureProfile)(params::NamedTuple, z::FT) where {FT}
    R_d = params.R_d
    g = params.g
    p₀ = params.pₒ

    # Scale height for surface temperature
    H_sfc = R_d * profile.T_virt_surf / g
    H_t = profile.H_t
    z′ = z / H_t
    tanh_z′ = tanh(z′)

    ΔTv = profile.T_virt_surf - profile.T_min_ref
    Tv = profile.T_virt_surf - ΔTv * tanh_z′

    ΔTv′ = ΔTv / profile.T_virt_surf
    p = -H_t * (z′ + ΔTv′ * (log(1 - ΔTv′ * tanh_z′) - log(1 + tanh_z′) + z′))
    p /= H_sfc * (1 - ΔTv′^2)
    p = p₀ * exp(p)
    return (Tv, p)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module BERL

using Cambrian
using NeuroEvolution
using CartesianGeneticProgramming
export atari_envs, gym_envs, environments, algorithms
atari_envs=[]
gym_envs=[]

algs = ["CGP", "NEAT", "HyperNEAT"]
envs = ["xor", "iris", "gym", "pybullet"]#, "atari_ram"]

include("environments/BERL_env.jl")
include("step.jl")

environments = Dict()
algorithms = Dict()

for e in envs
    include("environments/$e.jl")
end

for a in algs
    include("algorithms/$a/_setup.jl")
end

include("core.jl")
include("berl_run.jl")
include("plotting.jl")

end # module


text:
using Random
using SubspacePreservingSparsification: separated_min_max

v = rand(MersenneTwister(0), Float32, 20)
@test separated_min_max(v, 0.2f0) == (0.009370327f0, 0.14352024f0, 0.23794222f0, 0.9680481f0)
@test separated_min_max([1, 2, 3], 2) == (1, 2, 2, 3)


text:
@kwdef mutable struct DMultHandler <: AbstractHandler
    id = 0
end

@register(DMultHandler, *, [DBit, DBit] => [Number])

function (handler::DMultHandler)(netlist::Netlist,
                                 inputs::Vector{Variable},
                                 outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("Deterministic multiply does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN d_mult$(handler.id)
        determ_mult d_mult$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END d_mult$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end

@kwdef mutable struct DMultFxpHandler <: AbstractHandler
    id = 0
end

@register(DMultFxpHandler, *, begin
    [DBit, Number] => [Number]
    [Number, DBit] => [Number]
end)

function (handler::DMultFxpHandler)(netlist::Netlist,
                                    inputs::Vector{Variable},
                                    outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("Deterministic multiply does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN d_mult_fxp$(handler.id)
        determ_mult_fxp #(
                .BIT_WIDTH(BIT_WIDTH)
            ) d_mult_fxp$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END d_mult_fxp$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end


text:
using ThreeBodyDecay
using Test

@testset "Border and Flat Dalitz" begin
    ms = ThreeBodyMasses(6.0,1.1,1.3,1.5)
    fDP = flatDalitzPlotSample(ms)
    b = border31(ms)
    # make better test!
    #   It would be great to check if all points in `dDP` are inside of `b`
    @test true
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay