text:
module DoseResponse

    using LsqFit
    import LsqFit: LsqFitResult
    using StatsBase
    using GLM
    using Distributions
    using LinearAlgebra

    struct DoseResponseModel{Y <: AbstractVector, X <: AbstractVector, P,R,J,W <: AbstractArray} <: StatisticalModel
        result::LsqFitResult{P,R,J,W}
        response::Y
        dose::X
    end

    function DoseResponseModel(response::AbstractVector{T}, dose::AbstractVector{T}, theta0::AbstractVector{T}) where T <: Real
        fit = curve_fit(generalized_logistic_5PL!, jacobian_5PL!, dose, response, theta0, inplace=true)
        DoseResponseModel(fit, response, dose)
    end
    function DoseResponseModel(response::AbstractVector{T}, dose::AbstractVector{T}) where T <: Real
        theta0 = zeros(eltype(response), 5)
        initial_values!(theta0, response, dose)
        fit = curve_fit(generalized_logistic_5PL!, jacobian_5PL!, dose, response, theta0, inplace=true)
        DoseResponseModel(fit, response, dose)
    end
    function initial_values!(theta0::AbstractVector{T}, response::AbstractVector{T}, dose::AbstractVector{T}) where T <: Real
        ymin, ymax = minimum(response), maximum(response)
        theta0[1] = ymin - 0.001 * (ymax-ymin)
        theta0[2] = ymax + 0.001 * (ymax-ymin)

        ytrans = log.((theta0[2] .- response) ./ (response .- theta0[1]))
        idx = findall(dose .> 0)
        m = lm([ones(length(idx)) log.(x[idx])], ytrans[idx])
        theta0[3] = -coef(m)[2]
        theta0[4] = -coef(m)[1]/(-theta0[3])
        theta0[5] = 0.0
    end

    StatsBase.coef(m::DoseResponseModel) = m.result.param
    StatsBase.coefnames(m::DoseResponseModel) = ["A","K","B","M","logv"]
    StatsBase.confint(m::DoseResponseModel) = confidence_interval(m.result)
    StatsBase.dof(m::DoseResponseModel) = length(m.result.param) + 1
    StatsBase.nobs(m::DoseResponseModel) = length(m.result.resid)
    StatsBase.rss(m::DoseResponseModel) = sum(abs2, m.result.resid)
    StatsBase.weights(m::DoseResponseModel) = m.result.wt
    StatsBase.residuals(m::DoseResponseModel) = m.result.resid
    StatsBase.dof_residual(m::DoseResponseModel) = nobs(m) - length(coef(m))
    mse(m::DoseResponseModel) = rss(m)/dof_residual(m)

    function StatsBase.stderror(m::DoseResponseModel; rtol::Real=NaN, atol::Real=0)
        # computes standard error of estimates from
        #   fit   : a LsqFitResult from a curve_fit()
        #   atol  : absolute tolerance for approximate comparisson to 0.0 in negativity check
        #   rtol  : relative tolerance for approximate comparisson to 0.0 in negativity check
        covar = LsqFit.estimate_covar(m.result)
        # then the standard errors are given by the sqrt of the diagonal
        vars = diag(covar)
        vratio = minimum(vars)/maximum(vars)
        if !isapprox(vratio, 0.0, atol=atol, rtol=isnan(rtol) ? Base.rtoldefault(vratio, 0.0, 0) : rtol) && vratio < 0.0
            error("Covariance matrix is negative for atol=$atol and rtol=$rtol")
        end
        return sqrt.(abs.(vars))
    end

    function StatsBase.coeftable(m::DoseResponseModel; level::T=0.95) where T <: Real
        cc = coef(m)
        se = stderror(m)
        tt = cc ./ se
        p = ccdf.(Ref(FDist(1, dof_residual(m))), abs2.(tt))
        ci = se*quantile(TDist(dof_residual(m)), (1-level)/2)
        levstr = isinteger(level*100) ? string(Integer(level*100)) : string(level*100)
        CoefTable(hcat(cc,se,tt,p,cc+ci,cc-ci),
                ["Coef.","Std. Error","t","Pr(>|t|)","Lower $levstr%","Upper $levstr%"],
                coefnames(m), 4, 3)
    end

    StatsBase.loglikelihood(m::DoseResponseModel) = -nobs(m)/2 * (log(2π) + log(rss(m)) - log(nobs(m))+1)

    StatsBase.predict(m::DoseResponseModel) = generalized_logistic_5PL(m.dose, m.result.param)
    StatsBase.predict(m::DoseResponseModel{Y,X,P,R,J,W}, newX::X) where {Y<:AbstractVector,X<:AbstractVector,P,R,J,W<:AbstractArray} = generalized_logistic_5PL(newX, m.result.param)
    StatsBase.mss(m::DoseResponseModel) = sum(abs2, predict(m) .- mean(m.response))
    StatsBase.r2(m::DoseResponseModel) = 1 - rss(m)/sum(abs2,m.response .- mean(m.response))
    StatsBase.islinear(m::DoseResponseModel) = false

    ## 4PL
    function generalized_logistic_4PL!(F, x,theta)
        F .= ifelse.(isfinite.(log.(x)), theta[1] .+ (theta[2] - theta[1]) ./ (1.0 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))), theta[1] + (theta[2] - theta[1]) * ((theta[3] < 0) && ((theta[2] - theta[1]) > 0)))
    end
    function generalized_logistic_4PL(x,theta)
        ifelse.(isfinite.(log.(x)), theta[1] .+ (theta[2] - theta[1]) ./ (1.0 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))), theta[1] + (theta[2] - theta[1]) * ((theta[3] < 0) && ((theta[2] - theta[1]) > 0)))
    end

    function jacobian_4PL!(J, x, theta)
        J[:,1] .= 1 .- (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4])))
        J[:,2] .= (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4])))
        J[:,3] .= ifelse.(isfinite.(log.(x)), (theta[2] - theta[1]) .* (log.(x) .- theta[4]) .* exp.(-theta[3] .* (log.(x) .- theta[4])) ./ (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4])) ), 0.0)
        J[:,4] .=  -(theta[2] - theta[1]) .* theta[3] .* exp.(-theta[3] .* (log.(x) .- theta[4])) ./ (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4])) )
    end

    ## 5PL
    function generalized_logistic_5PL!(F, x, theta)
        F .= ifelse.(isfinite.(log.(x)), theta[1] .+ (theta[2] - theta[1]) ./ (1.0 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) .^ (1 / exp(theta[5])), theta[1] + (theta[2] - theta[1]) * ((theta[3] < 0) && ((theta[2] - theta[1]) > 0)))
    end
    function generalized_logistic_5PL(x, theta)
        ifelse.(isfinite.(log.(x)), theta[1] .+ (theta[2] - theta[1]) ./ (1.0 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) .^ (1 / exp(theta[5])), theta[1] + (theta[2] - theta[1]) * ((theta[3] < 0) && ((theta[2] - theta[1]) > 0)))
    end

    function jacobian_5PL!(J, x, theta)
        J[:,1] .= 1 .- (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) .^ (-1 / exp(theta[5]))
        J[:,2] .= (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) .^ (-1 / exp(theta[5]))
        J[:,3] .= ifelse.(isfinite.(log.(x)), (theta[2] - theta[1]) .* (log.(x) .- theta[4]) .* exp.(-theta[3] .* (log.(x) .- theta[4])) ./ (exp(theta[5]) .* (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) .^ (1/exp(theta[5]) + 1)), 0.0)
        J[:,4] .=  -(theta[2] - theta[1]) .* theta[3] .* exp.(-theta[3] .* (log.(x) .- theta[4])) ./ (exp(theta[5]) .* (1 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) .^ (1/exp(theta[5]) + 1))
        J[:,5] .= (theta[2] - theta[1]) .* log.(1 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) ./ (exp(theta[5]) .*(1 .+ exp.(-theta[3] .* (log.(x) .- theta[4]))) .^ (1/exp(theta[5])))
    end

end # module


text:
@testset "matrix difference" begin
    modinds = ([1, 2], [1, 2])
    addinds = ([3, 4], [3])
    reminds = ([5, 6], [5])
    modvals = [1, 2, 3, 4]
    addvals = [1, 2]
    remvals = [3, 4]
    a = MatrixDifference(
        modinds, addinds, reminds, modvals, view(addvals, :), view(remvals, :)
    )
    b = MatrixDifference(modinds, addinds, reminds, sparse(modvals), addvals, remvals)
    c = MatrixDifference(modinds, addinds, reminds, modvals, addvals, remvals)

    @testset "constructors" begin
        @test isa(a, MatrixDifference) && isa(b, MatrixDifference)
        @test (
            a.modinds == modinds &&
            a.addinds == addinds &&
            a.reminds == reminds &&
            a.modvals == modvals &&
            a.addvals == addvals &&
            a.remvals == remvals
        )
        @test (
            b.modinds == modinds &&
            b.addinds == addinds &&
            b.reminds == reminds &&
            b.modvals == sparse(modvals) &&
            b.addvals == addvals &&
            b.remvals == remvals
        )
    end

    @testset "equality operator" begin
        @test a == a
        @test a == b && b == a
        @test a == b && b == c && a == c
    end

    @testset "hash" begin
        @test hash(a) == hash(a)
        @test a == b && hash(a) == hash(b)
    end

    @testset "accessors" begin
        @test modified_indices(a) == modinds
        @test added_indices(a) == addinds
        @test removed_indices(a) == reminds
        @test modified(a) == modvals
        @test added(a) == addvals
        @test removed(a) == remvals
    end

    @testset "difference" begin
        A = [1 0 1; 0 1 0; 0 0 1]
        B = [1 1 1 1; 1 1 1 1]
        ia = [1, 2, 5]
        ja = [2, 8, 11]
        ib = [1, 8]
        jb = [2, 3, 4, 11]
        E = Matrix(undef, 0, 0)
        e = []
        ei = Int[]
        se = sparse(e)
        @test diff(E, E) == MatrixDifference((ei, ei), (ei, ei), (ei, ei), se, e, e)
        @test diff(E, A) ==
              MatrixDifference((ei, ei), (ei, ei), (ei, ei), se, vec(A), vec(E))
        @test diff(A, E) ==
              MatrixDifference((ei, ei), (ei, ei), (ei, ei), se, vec(A), vec(E))
        @test diff(A, A) == MatrixDifference(
            ([1, 2, 3], [1, 2, 3]),
            (ei, ei),
            (ei, ei),
            sparse([0, 0, 0, 0, 0, 0, 0, 0, 0]),
            e,
            e,
        )
        @test diff(A, B) == MatrixDifference(
            ([1, 2], [1, 2, 3]),
            (ei, [4]),
            ([3], ei),
            sparse([0, -1, -1, 0, 0, -1]),
            [1, 1],
            [0, 0, 1],
        )
        @test diff(B, A) == MatrixDifference(
            ([1, 2], [1, 2, 3]),
            ([3], ei),
            (ei, [4]),
            sparse([0, 1, 1, 0, 0, 1]),
            [0, 0, 1],
            [1, 1],
        )
        @test diff(A, A, ia, ja, ia, ja) == MatrixDifference(
            ([1, 2, 5], [2, 8, 11]),
            (ei, ei),
            (ei, ei),
            sparse([0, 0, 0, 0, 0, 0, 0, 0, 0]),
            e,
            e,
        )
        @test diff(A, B, ia, ja, ib, jb) == MatrixDifference(
            ([1], [2, 11]),
            ([8], [3, 4]),
            ([2, 5], [8]),
            sparse([0, 0]),
            [1, 1, 1, 1, 1, 1],
            [0, 0, 0, 1, 0, 0, 1],
        )
        @test diff(B, A, ib, jb, ia, ja) == MatrixDifference(
            ([1], [2, 11]),
            ([2, 5], [8]),
            ([8], [3, 4]),
            sparse([0, 0]),
            [0, 0, 0, 1, 0, 0, 1],
            [1, 1, 1, 1, 1, 1],
        )
    end
end

ext: jl


text:
import Gadfly
using Gadfly

function Gadfly.plot(data::Vector{Float64}, changepoints::Vector{Int})
    p = plot(layer(y=data, Geom.line, Theme(default_color=colorant"black")),
             layer(xintercept=changepoints, Geom.vline, Theme(default_color=colorant"red")))
    return p
end

function Gadfly.plot(crops_object::Dict{ASCIIString, Array})
    x = crops_object["number"]
    y = crops_object["constrained"]
    p = plot(x=x, y=y, Guide.title("Elbow Plot for Range of Penalties"),
             Guide.xlabel("Number of Segments"),
             Guide.ylabel("Cost"),
             Geom.line,
             Geom.point,
             Theme(default_color=colorant"black"))
    return p
end


text:
module aoc_21

export IngredientList, read_input
export ingredients, allergies

const ingredients_r = r"([a-z\s]*)\(contains ([a-z,\s]*)\)"

struct IngredientList
    ingredients::Set{String}
    allergies::Set{String}
end

ingredients(list::IngredientList) = list.ingredients
allergies(list::IngredientList) = list.allergies

import Base.intersect
function intersect(s::IngredientList, lists::IngredientList...)
    i = intersect( s.ingredients, ingredients.(lists)... )
    a = intersect( s.allergies, allergies.(lists)... )
    IngredientList( i, a )
end

const IngredientLists = Set{IngredientList}

function read_input(filename="input.txt")
    lines = readlines(filename)
    lists = IngredientLists()
    # count = 0
    for line in lines
        m = match(ingredients_r, line)
        ingredients = split(strip(m.captures[1])," ")
        # count += length(ingredients)
        allergies = split(strip(m.captures[2]),", ")
        push!(lists, IngredientList(Set(ingredients), Set(allergies)) )
    end
    # println(count)
    lists
end

function build_map(lists::IngredientLists)
    allergies_to_lists = Dict{String,IngredientLists}()
    for list in lists
        for allergy in list.allergies
            current_lists = get(allergies_to_lists, allergy, IngredientLists())
            push!(current_lists, list)
            allergies_to_lists[ allergy ] = current_lists
        end
    end
    allergies_to_lists
end

function constrain(lists::IngredientLists,
                   allergies_to_lists::Dict{String,IngredientLists} = build_map(lists))
    ingredients_to_allergies = Dict{String,String}()

    change = true

    #while( length(ingredients_to_allergies) < length(allergies_to_lists) )
    while( change )

    change = false

    for (allergy, lists) in allergies_to_lists
        intersection = intersect(lists...)
        if length(intersection.ingredients) == 1 && first(intersection.allergies) == allergy
            ingredient = first(intersection.ingredients)
            ingredients_to_allergies[ ingredient ] = allergy
            # pop!.( ingredients.(lists), ingredient, nothing )
            # pop!.( allergies.(lists), allergy, nothing )
            pop!( allergies_to_lists, allergy )
            for list in Iterators.flatten( values(allergies_to_lists) )
                pop!( list.ingredients, ingredient, nothing )
                pop!( list.allergies, allergy, nothing )
            end
            change = true
        end
    end

    # println(ingredients_to_allergies)

    end

    for (allergy, lists) in allergies_to_lists
        # remaining_ingredient_list = intersect(lists...)
        # allergies_to_lists[ allergy ] = Set([ remaining_ingredient_list ])
        # [pop!.( ingredients.(lists) , ingredient, nothing) for ingredient in ingredients(remaining_ingredient_list)]
    end

    (lists, ingredients_to_allergies, allergies_to_lists)
end

function solve(input = read_input())
    out = constrain( deepcopy(input) )
    allergens = keys( out[2] )
    # allergens = union( keys(out[2]), ingredients.( first.( values( out[3] ) ) ) ... )
    input_a = [input...]
    [ pop!.( ingredients.(input_a) , allergen, nothing ) for allergen in allergens]
    count = sum(length.(ingredients.(input_a)))
    list_of_allergens = join(first.( sort( collect(out[2]) , by=p->p.second) ),",")
    (count, list_of_allergens)
    # 2262
end

part1() = solve()[1]
part2() = solve()[2]

end


text:
type XGBoost <: DecisionTree
    n_estimators::Int64
    learning_rate::Float64
    min_samples_split::Int64
    min_impurity::Float64
    max_depth::Int64
    init_estimate::Union{String, Array}
    loss::LogisticLoss
    trees::Vector{XGBoostRegressionTree}
end

function XGBoost(;
                 n_estimators::Int64 = 2,
                 learning_rate::Float64 = 0.001,
                 min_samples_split::Int64 = 2,
                 min_impurity::Float64 = 0.00001,
                 max_depth::Int64 = 2,
                 init_estimate = "nothing",
                 loss::LogisticLoss = LogisticLoss())
    trees = []
    single_boost_tree = XGBoostRegressionTree(min_samples_split = min_samples_split,
            min_gain = min_impurity, max_depth = max_depth, loss = loss)
    for i in 1:n_estimators
        push!(trees, single_boost_tree)
    end
    return XGBoost(n_estimators, learning_rate,
                   min_samples_split, min_impurity,
                   max_depth, init_estimate, loss, trees)
end

function train!(model::XGBoost, X::Matrix, y::Vector)
    y = one_hot(y)
    n_sample, n_feature = size(X)
    model.init_estimate = mean(y, 1)
    y_pred = repmat(model.init_estimate, size(y,1),1)
    for i in 1:model.n_estimators
        y_and_pred = hcat(y, y_pred)
        print("$(y_and_pred)")
        train!(model.trees[i], X, y_and_pred)
        update_pred = predict(model.trees[i],X)
        y_pred += model.learning_rate * update_pred
        progress = 1/model.n_estimators * i * 100
        #println("Progress: $(progress) %")
    end
end

function predict(model::XGBoost,
                 x::Matrix)
    n_sample = size(x)
    flag = 1
    y_pred = 0
    for tree in model.trees
        update_pred = predict(tree,x)
        update = model.learning_rate * update_pred
        if flag == 1
            y_pred = update
            flag = 0
        else
            y_pred += update
        end
    end
    res = zeros(size(y_pred,1))
    for i in 1:size(y_pred,1)
        exp_pred = exp(y_pred[i,:])
        y_pred[i,:] = exp_pred / sum(exp_pred)
        print("$(y_pred[i,:])")
        res[i] = indmax(y_pred[i,:])
    end
    return res-1
end

function test_xgboost()
    X_train, X_test, y_train, y_test = make_iris()
    model = XGBoost()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    print("$(predictions)")
    print("classification accuracy", accuracy(y_test, predictions))
end


text:
mutable struct DirichletDist
    #=
    Dirichlet Distribution
    p(mu | alpha)
    = gamma(sum(alpha))
    * Prod_k[ mu_k^(alpha_k-1)]
    / gamma(alpha_1) / ... / gamma(alpha_k)
    =#
    _alpha::AbstractArray{Float64,1}
end

function pdf(dirichlet::DirichletDist, mu::AbstractArray{Float64,1})
    N = size(mu)[1]
    @assert N == size(dirichlet._alpha)[1]

    sum_alpha = sum(dirichlet._alpha)
    prod_k = 1.0
    gamma_k = 1.0
    for k = 1:N
        alpha_k = dirichlet._alpha[k]
        prod_k *= (mu[k]^(alpha_k - 1.0))
        gamma_k *= gamma(alpha_k)
    end
    return gamma(sum_alpha) * prod_k / gamma_k
end

function pdf(dirichlet::DirichletDist, mu::AbstractArray{Float64,2})
    N = size(mu)[1]
    n_samples = size(mu)[2]
    @assert N == size(dirichlet._alpha)[1]

    return [pdf(dirichlet, mu[:, i]) for i = 1:n_samples]
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
using Documenter, GeomDAE

makedocs(
    sitename = "GeomDAE.jl",
    format = :html,
    pages = ["Home" => "index.md",
             "Tutorial" => "tutorial.md",
             "Modules"  => [
                "Equations"         => "modules/equations.md",
                "Integrators"       => "modules/integrators.md",
                "Linear Solvers"    => "modules/solvers_linear.md",
                "Nonlinear Solvers" => "modules/solvers_nonlinear.md",
                "Tableaus"          => "modules/tableaus.md"]
             ]
)


text:
@testset "Spaces" begin
  k = FlintQQ

  Qx, x = PolynomialRing(FlintQQ, "x")
  K1, a1 = NumberField(x^2 - 2, "a1")
  K2, a2 = NumberField(x^3 - 2, "a2")

  K1t, t = PolynomialRing(K1, "t")
  F = GF(3)

  Hecke.change_base_ring(::FlintRationalField, ::Hecke.gfp_mat) = error("asd")
  @test_throws ErrorException quadratic_space(FlintQQ, F[1 2; 2 1])

  Hecke.change_base_ring(::FlintRationalField, x::Hecke.gfp_mat) = x
  @test_throws ErrorException quadratic_space(FlintQQ, F[1 2; 2 1])

  L, b = NumberField(t^2 + a1)

  for K in [k, K1, K2, L]
    V = @inferred quadratic_space(K, 2)
    @test_throws ArgumentError quadratic_space(K, -1)

    V = @inferred quadratic_space(K, identity_matrix(FlintZZ, 2))
    @test (@inferred gram_matrix(V)) == identity_matrix(K, 2)
    @test (@inferred dim(V)) == 2
    @test (@inferred rank(V)) == 2
    @test @inferred isregular(V)
    @test (@inferred involution(V)) === identity

    V = @inferred quadratic_space(K, zero_matrix(K, 2, 2))
    @test (@inferred gram_matrix(V)) == zero_matrix(K, 2, 2)
    @test (@inferred gram_matrix(V)) isa dense_matrix_type(K)
    @test (@inferred dim(V)) == 2
    @test (@inferred rank(V)) == 0
    @test @inferred isquadratic(V)
    @test @inferred !ishermitian(V)
    @test (@inferred fixed_field(V)) === K

    @test_throws ArgumentError quadratic_space(K, zero_matrix(K, 1, 2))

    M = identity_matrix(K, 2)
    M[1, 2] = 2
    M[2, 1] = 1
    @test_throws ArgumentError quadratic_space(K, M)

    # Determinant & discrimimant

    V = quadratic_space(K, 2)
    @test (@inferred discriminant(V)) == -1
    @test (@inferred det(V)) == 1

    V = quadratic_space(K, 4)
    @test (@inferred discriminant(V)) == 1
    @test (@inferred det(V)) == 1

    M = identity_matrix(K, 2)
    M[1, 2] = 2
    M[2, 1] = 2
    V = quadratic_space(K, M)
    @test (@inferred discriminant(V)) == 3
    @test (@inferred det(V)) == -3

    # Gram matrix

    M = identity_matrix(K, 2)
    M[1, 2] = 2
    M[2, 1] = 2
    V = @inferred quadratic_space(K, M)
    N = zero_matrix(K, 4, 2)
    @test (@inferred gram_matrix(V, N)) == zero_matrix(K, 4, 4)

    N = identity_matrix(QQ, 2)
    @test (@inferred gram_matrix(V, N)) == M

    N = zero_matrix(K, 4, 4)
    @test_throws ArgumentError gram_matrix(V, N)

    v = [[1, 0], [0, 1]]
    @test (@inferred gram_matrix(V, v) == M)

    v = [[1, 0, 0], [0, 1]]
    @test_throws ErrorException gram_matrix(V, v)

    B = @inferred orthogonal_basis(V)
    @test isdiagonal(gram_matrix(V, B))

    D = @inferred diagonal(V)
    @test length(D) == 2
    @test issetequal(D, map(K, [1, -3]))

    M = rand(MatrixSpace(K, 4, 4), -10:10)
    M = M + transpose(M)
    while iszero(det(M))
      M = rand(MatrixSpace(K, 4, 4), -10:10)
      M = M + transpose(M)
    end

    V = quadratic_space(K, M)

    F, S = @inferred Hecke._gram_schmidt(M, identity)
    @test gram_matrix(V, S) == F

    M[1, 1] = 0
    M[1, 2] = 0
    M[1, 3] = 0
    M[1, 4] = 0
    M[2, 1] = 0
    M[3, 1] = 0
    M[4, 1] = 0

    V = quadratic_space(K, M)

    @test_throws ErrorException Hecke._gram_schmidt(M, identity)
    F, S = @inferred Hecke._gram_schmidt(M, identity, false)
    @test gram_matrix(V, S) == F
  end

  fl, T =  Hecke.isequivalent_with_isometry(quadratic_space(QQ, matrix(QQ, 2, 2, [4, 0, 0, 1])),
                                            quadratic_space(QQ, matrix(QQ, 2, 2, [1, 0, 0, 1])))
  @test fl

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 2, 2, [1, 0, 0, 3]);
  V = quadratic_space(K, D)
  fl, T = Hecke.isequivalent_with_isometry(V, V)
  @test fl
end


text:
@syms 𝑥
@polyvar 𝑦

function poly(eq, v::Pair)
    α, β = value(first(v)), last(v)
    terms = collect_powers(eq, α)

    if !all(isone, denominator.(collect(keys(terms))))
        error("fractional power not supported")
    end

    sum([val*β^numerator(k) for (k,val) in terms]; init=0)
end

function poly(eq, β)
    xs = get_variables(eq)

    if length(xs) > 1
        error("more than one implicit variable.")
    elseif length(xs) == 0
        return nothing
    else
        return poly(eq, xs[1] => β)
    end
end

poly(eq) = poly(eq, 𝑦)

function sym(p::AbstractPolynomialLike, v::Pair)
    β, α = first(v), value(last(v))
    sum([α^maxdegree(t,β)*c for (t,c) in zip(terms(p), coefficients(p))]; init=0)
end

sym(p::AbstractPolynomialLike, α) = sym(p, var(p) => α)
sym(p::AbstractPolynomialLike) = sym(p, var(p) => 𝑥)

###############################################################################

# pox (power-of-x) is a symbolic function to keep track of the powers of x
# pox(k,n) means k*x^n
@syms pox(k, n)

is_pox(x) = istree(x) && operation(x)==pox
is_not_pox(x) = !is_pox(x)

get_coef(p) = is_pox(p) ? arguments(p)[1] : p
get_power(p) = is_pox(p) ? arguments(p)[2] : 0

replace_x(eq, x) = substitute(eq, Dict(x => pox(1,1)))

iscomplex(x) = x isa Complex

count_rule1 = @rule ^(pox(~k, ~n1), ~n2) => isequal(~k,1) ? pox(1, ~n1 * ~n2) : pox(^(~k,~n2), ~n1 * ~n2)
count_rule2 = @rule pox(~k1, ~n1) * pox(~k2, ~n2) => pox(~k1 * ~k2, ~n1 + ~n2)
count_rule3 = @acrule pox(~k, ~n) * ~u::is_not_pox => pox(~k * ~u, ~n)

"""
    collect_powers separates the powers of x in eq (a polynomial) and returns
    a dictionary of power => term
"""
function collect_powers(eq, x)
    eq = expand(expand_derivatives(eq))
    eq = replace_x(eq, x)
    #eq = Prewalk(PassThrough(count_rule1))(eq)
    eq = Fixpoint(Prewalk(PassThrough(Chain([count_rule1, count_rule2, count_rule3]))))(eq)

    if !istree(eq)
        return Dict{Any, Any}(0 => eq)
    elseif is_pox(eq)
        return Dict{Any, Any}(get_power(eq) => get_coef(eq))
    else
        eqs = Dict{Any, Any}()
        for term in arguments(eq)
            n = get_power(term)
            if haskey(eqs, n)
                eqs[n] = eqs[n] + get_coef(term)
            else
                eqs[n] = get_coef(term)
            end
        end

        return eqs
    end
end


text:
using VoronoiCells
import VoronoiCells: LEFT, RIGHT, LOWER, UPPER
using VoronoiDelaunay
using Base.Test

# Area of square
S = [Point2D(LEFT,LOWER) ; Point2D(RIGHT,LOWER) ; Point2D(RIGHT,UPPER) ; Point2D(LEFT,UPPER)]
@test_approx_eq polyarea(S) 1.0

# Area of triangle
T1 = [Point2D(LEFT,LOWER) ; Point2D(RIGHT,LOWER) ; 0.5*(Point2D(RIGHT,UPPER) + Point2D(LEFT,UPPER))]
@test_approx_eq polyarea(T1) 0.5

T2 = [Point2D(LEFT,LOWER) ; Point2D(LEFT+0.5,LOWER) ;
Point2D(LEFT,LOWER+0.5)]
@test_approx_eq polyarea(T2) 0.5^3


text:
using Test

using LinearAlgebra, Combinatorics, Statistics
using OrdinaryDiffEq, ForwardDiff, DiffResults
using DoubleFloats

include("../src/type.jl")
include("testfunctions.jl")

include("../src/metrics/Minkowski.jl")
include("../src/metrics/KerrSchild.jl")
include("../src/geosol.jl")
include("../src/broyden.jl")
include("../src/outlier.jl")
include("../src/squirrel.jl")

η = X -> ημν(typeof(X[1]),4)

@testset "All tests:" begin

    @time @testset "Minkowski tests:" begin
            include("minkowski_test.jl")
        end

    @time @testset "geosol tests:" begin
            include("geosol_test.jl")
        end

    @time @testset "broyden tests:" begin
            include("broyden_test.jl")
        end

    @time @testset "outlier detection tests:" begin
            include("outlier_test.jl")
        end

    @time @testset "squirrel tests:" begin
            include("squirrel_test.jl")
        end

    @time @testset "seval tests:" begin
            include("seval_test.jl")
    end

end

nothing


text:
mutable struct NBarsEffect{T<:AbstractFloat} <: ConfigurableEffect
    config::EffectConfig
    spec::Array{T}
    maxspec::AbstractFloat
    colors::Array{UInt8, 2}
    backgroundcolor::Array{UInt8}
    normspec::Array{T}
    function NBarsEffect{T}(leddata, config, ana) where T<:AbstractFloat
        spec, maxspec = binFFT(ana, length(leddata.channels), floor(Int, length(ana.spec_bufs[ana.spec_buf_order[1]])/2))
        colors = zeros(UInt8, length(leddata.channels), 3)
        prim_color = config.primary_color
        for i in 1:size(colors, 1)
            colors[i, :] .= round.(UInt8, hsl_to_rgb(prim_color.h/360, prim_color.s*(0.5*i/size(colors,1)+0.5), prim_color.l))
        end
        sec_color = config.secondary_color
        backgroundcolor = round.(UInt8, hsl_to_rgb(sec_color.h/360, sec_color.s, sec_color.l))
        normspec = zeros(length(spec))
        return new(config, spec, maxspec, colors, backgroundcolor, normspec)
    end
end

function update!(leddata::LEDArray, ana::AudioAnalysis, effect::NBarsEffect)
    # Abstract away updating arrays in the object
    update_vars!(leddata, ana, effect)

    for i in eachindex(leddata.channels)

        chan = leddata.channels[i]
        crossover = max(0, min(size(chan,1), floor(Int, effect.normspec[i])))
        for j in 1:crossover
            for k in 1:size(chan,2)
                chan[j,k] = effect.colors[i,k]
            end
        end
        for j in (crossover+1):size(chan,1)
            for k in 1:size(chan, 2)
                chan[j,k] = effect.backgroundcolor[k]
            end
        end
        update!(chan)
    end
end

function process_colors!(leddata::LEDArray, ana::AudioAnalysis, effect::NBarsEffect)
    if !(size(effect.colors) == (length(leddata.channels, 3)))
        effect.colors = zeros(UInt8, length(leddata.channels), 3)
    end

    prim_color = effect.config.primary_color

    for i in 1:size(effect.colors, 1)
        effect.colors[i, :] .= round.(UInt8, hsl_to_rgb(prim_color.h/360, prim_color.s*(0.5*i/size(effect.colors,1)+0.5), prim_color.l))
    end

    sec_color = effect.config.secondary_color
    effect.backgroundcolor = round.(UInt8, hsl_to_rgb(sec_color.h/360, sec_color.s, sec_color.l))
end

function update_vars!(leddata::LEDArray, ana::AudioAnalysis, effect::NBarsEffect)
    effect.spec, effect.maxspec = binFFT(ana, length(leddata.channels), floor(Int, length(ana.spec_bufs[ana.spec_buf_order[1]])/2))
    if length(effect.normspec) != length(effect.spec)
        effect.normspec = zeros(length(effect.spec))
    end

    for i in eachindex(effect.spec)
        ref_len = size(leddata.channels[i],1)
        tmp = effect.spec[i]/fft_scale[]*ref_len
        if !isnan(tmp)
            effect.normspec[i]=abs(tmp)
        end
    end
end


text:
struct ALRandom
    n_samples::Int
    nn_stability::Int
    max_gamma::Float64
    adjusted::Bool
    C_steps::Int
end

function ALRandom(;n_samples, nn_stability, max_gamma, adjusted, C_steps)
    ALRandom(n_samples, nn_stability, max_gamma, adjusted, C_steps)
end

function Base.string(m::ALRandom)
    return "ALRandom_n=$(m.n_samples)_nns=$(m.nn_stability)_adj=$(m.adjusted)"
end

function find_parameters(method::ALRandom, data, labels, oracle, solver, qm, seed)
    @show "Running $method"

    nnpool = NNPool(data, k=method.nn_stability)
    sample = StatsBase.sample(1:size(data,2), method.n_samples, replace=false)
    for i in sample
        add_to_pools!(nnpool, i, ask_oracle(oracle, i), method.nn_stability)
    end

    gamma = find_gamma(nnpool, data, method.nn_stability; max_γ=method.max_gamma, adjusted=method.adjusted)
    @show "Found gamma $gamma"
    max_C = find_max_feasible_C(data, gamma, solver)
    min_C = find_min_feasible_C(data, gamma, solver)

    C_range, model_range = calculate_C_grid(data, labels, oracle, gamma, min_C, max_C, method.C_steps, solver)
    qm = calculate_qm(nnpool, data, labels, model_range, qm)
    qm_best_idx = findmax(qm)[2]
    C = C_range[qm_best_idx]
    return (gamma, C)
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
#Worked with Bill and Aleeze

    using JLD2, Random, LinearAlgebra, Statistics, CSV, DataFrames, FreqTables, Distributions,
    FreqTables, Optim, HTTP, GLM

    #:::::::::::::::::::
    # question 1   :::::
    #:::::::::::::::::::

function ps3()

    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS3-gev/nlsw88w.csv"
    df = CSV.read(HTTP.get(url).body, DataFrame)
    X = [df.age df.white df.collgrad]
    Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
    df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
    y = df.occupation

    K = size(X,2)+1
    J = length(unique(y))

        function mlogit(Beta, X, Z,y)

            K = size(X,2)+1
            J = length(unique(y))
            N = length(y)
            bigY = zeros(N,J)

                for j=1:J
                    bigY[:,j] = y.==j
                end

            bigBeta = [reshape(Beta,K,J-1) zeros(K)]
            bigZ = zeros(N,J)

                for j=1:J
                    bigZ[:,j] = Z[:,j]-Z[:,J]
                end

            num = zeros(N,J)
            dem = zeros(N)

                for j=1:J
                    XZ = cat(X,bigZ[:,j],dims=2)
                    num[:,j] = exp.(XZ*bigBeta[:,j])
                    dem .+= num[:,j]
                end

            P = num./repeat(dem,1,J)

            loglike = -sum( bigY.*log.(P) )

                return loglike
        end

        startHere0 = rand((J-1)*K)
        startNearSol = 1.0000001 * [0.04885858084422522, 0.0009042465857047558, -2.144317825388887,
         0.4425753938680137, 0.04448332173290848, 0.7276857878249758, -3.1545819831597375, -0.028422550899673403,
          0.09091530580736103, -0.13094629047687026, -4.1213637842532025, 0.2578758066405763, 0.025242867291607733,
           0.8513538201292878, -4.026561140941997, -0.5232455580421997, 0.04201112040602136, -0.5535189061525868,
           -4.471599137048158, -0.97027896459939, 0.0856777779538767, -1.1140765374638752, -6.954753655923485,
           -0.3570730023566607, 0.08253331780409119, -0.5940641835229424, -5.325234363329876, -0.8111789371715972]

    mlogit_hat_optim = optimize(b-> mlogit(b,X,Z, y), startNearSol, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=200))
    println(mlogit_hat_optim.minimizer)

    #:::::::::::::::::::
    # question 2   :::::
    #:::::::::::::::::::

 # source: https://stats.idre.ucla.edu/stata/output/multinomial-logistic-regression/
    # the standard interpretation of the multinomial logit is that for a unit change in the predictor
    # variable, the logit of outcome m relative to the referent group is expected to change by its respective
    # parameter estimate (which is in log-odds units) given the variables in the model are held constant.

    # in our example it is the logit of occupation relative to a referenced occupation is expected to change by
    # the estimated parametr, ceteris paribus.

    #:::::::::::::::::::
    # question 3   :::::
    #:::::::::::::::::::
    function nested_logit(alpha, X, Z, y, nesting_structure)

        beta = alpha[1:end-3]
        lambda = alpha[end-2:end-1]
        gamma = alpha[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)

            for j=1:J
                bigY[:,j] = y.==j
            end
        bigBeta = [repeat(beta[1:K],1,length(nesting_structure[1])) repeat(beta[K+1:2K],1,length(nesting_structure[2])) zeros(K)]

        T = promote_type(eltype(X),eltype(alpha))
        num   = zeros(T,N,J)
        lidx  = zeros(T,N,J)
        dem   = zeros(T,N)

        for j=1:J
            if j in nesting_structure[1]
                lidx[:,j] = exp.( (X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[1] )
            elseif j in nesting_structure[2]
                lidx[:,j] = exp.( (X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[2] )
            else
                lidx[:,j] = exp.(zeros(N))
            end
        end

        for j=1:J
            if j in nesting_structure[1]
                num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[1][:]];dims=2).^(lambda[1]-1)
            elseif j in nesting_structure[2]
                num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[2][:]];dims=2).^(lambda[2]-1)
            else
                num[:,j] = lidx[:,j]
            end

            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end
        nesting_structure = [[1 2 3], [4 5 6 7]]
        startvals = [2*rand(2*size(X,2)).-1; 1; 1; .1]

        td2 = TwiceDifferentiable(alpha -> nested_logit(alpha, X, Z, y, nesting_structure), startvals; autodiff = :forward)
        α_hat_nlogit = optimize(td2, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=200))
        α_hat_nlogit_ad = α_hat_nlogit.minimizer

        H  = Optim.hessian!(td2, α_hat_nlogit_ad)
        α_hat_nlogit_ad_se = sqrt.(diag(inv(H)))

        println([α_hat_nlogit_ad α_hat_nlogit_ad_se])

        return nothing
    end

        ps3()


text:
# array.jl

export ndgrid_array

"""
    (xg, yg, ...) = ndgrid_array(v1, v2, ...)

Method to construct a tuple of (dense) `Array`s from a set of vectors.

This tuple can use a lot of memory so should be avoided in general!
It is provided mainly for testing and timing comparisons.

Each input should be an `AbstractVector` of some type.
The corresponding output Array will have the same element type.

This method provides similar functionality as Matlab's `ndarray` function
but is more general because the vectors can be any type.

# Examples

```jldoctest
julia> ndgrid_array(1:3, 1:2)
([1 1; 2 2; 3 3], [1 2; 1 2; 1 2])

julia> ndgrid(1:3, [:a,:b])
([1 1; 2 2; 3 3], [:a :b; :a :b; :a :b])
```
"""
ndgrid_array(vs::AbstractVector...) = Array.(ndgrid(vs...))


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a PTX-based test compiler, and generate reflection methods for it

PTXCompilerJob = CompilerJob{PTXCompilerTarget,TestCompilerParams}

struct PTXKernelState
    data::Int64
end
GPUCompiler.kernel_state_type(@nospecialize(job::PTXCompilerJob)) = PTXKernelState
@inline @generated ptx_kernel_state() = GPUCompiler.kernel_state_value(PTXKernelState)

# a version of the test runtime that has some side effects, loading the kernel state
# (so that we can test if kernel state arguments are appropriately optimized away)
module PTXTestRuntime
    using ..GPUCompiler
    import ..PTXKernelState

    function signal_exception()
        ptx_kernel_state()
        return
    end

    # dummy methods
    # HACK: if malloc returns 0 or traps, all calling functions (like jl_box_*)
    #       get reduced to a trap, which really messes with our test suite.
    malloc(sz) = Ptr{Cvoid}(Int(0xDEADBEEF))
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return
end
GPUCompiler.runtime_module(::PTXCompilerJob) = PTXTestRuntime

function ptx_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false,
                 minthreads=nothing, maxthreads=nothing, blocks_per_sm=nothing,
                 maxregs=nothing, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = PTXCompilerTarget(cap=v"7.0",
                               minthreads=minthreads, maxthreads=maxthreads,
                               blocks_per_sm=blocks_per_sm, maxregs=maxregs)
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function ptx_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function ptx_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function ptx_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function ptx_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    ptx_method = Symbol("ptx_$(method)")
    @eval begin
        $ptx_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $ptx_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function ptx_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = ptx_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end


text:
# juliap -O3 --check-bounds=no --math-mode=fast diffusion_2D_perf_gpu.jl
using Plots, Printf, AMDGPU

# macros to avoid array allocation
macro qx(ix,iy)  esc(:( -D_dx*(C[$ix+1,$iy+1] - C[$ix,$iy+1]) )) end
macro qy(ix,iy)  esc(:( -D_dy*(C[$ix+1,$iy+1] - C[$ix+1,$iy]) )) end

function compute!(C2, C, D_dx, D_dy, dt, _dx, _dy, size_C1_2, size_C2_2)
    ix = (workgroupIdx().x - 1) * workgroupDim().x + workitemIdx().x  # only workgroupIdx starts at 1
    iy = (workgroupIdx().y - 1) * workgroupDim().y + workitemIdx().y  # only workgroupIdx starts at 1
    if (ix<=size_C1_2 && iy<=size_C2_2)
        C2[ix+1,iy+1] = C[ix+1,iy+1] - dt*( (@qx(ix+1,iy) - @qx(ix,iy))*_dx + (@qy(ix,iy+1) - @qy(ix,iy))*_dy )
    end
    return
end

@views function diffusion_2D(; do_visu=false)
    # Physics
    Lx, Ly  = 10.0, 10.0
    D       = 1.0
    ttot    = 1e-4
    # Numerics
    BLOCKX  = 32
    BLOCKY  = 8
    GRIDX   = 16*24
    GRIDY   = 64*24
    nx, ny  = BLOCKX*GRIDX, BLOCKY*GRIDY # number of grid points
    nout    = 50
    # Derived numerics
    dx, dy  = Lx/nx, Ly/ny
    dt      = min(dx, dy)^2/D/4.1
    nt      = cld(ttot, dt)
    xc, yc  = LinRange(dx/2, Lx-dx/2, nx), LinRange(dy/2, Ly-dy/2, ny)
    D_dx    = D/dx
    D_dy    = D/dy
    _dx, _dy= 1.0/dx, 1.0/dy
    # Array initialisation
    C       = ROCArray(exp.(.-(xc .- Lx/2).^2 .-(yc' .- Ly/2).^2))
    C2      = copy(C)
    TBSize  = (BLOCKX, BLOCKY, 1)
    GSize   = (nx, ny, 1)
    size_C1_2, size_C2_2 = size(C,1)-2, size(C,2)-2
    t_tic = 0.0; niter = 0
    # Time loop
    for it = 1:nt
        if (it==11) t_tic = Base.time(); niter = 0 end
        wait( @roc groupsize = TBSize gridsize = GSize compute!(C2, C, D_dx, D_dy, dt, _dx, _dy, size_C1_2, size_C2_2) )
        C, C2 = C2, C # pointer swap
        niter += 1
        if do_visu && (it % nout == 0)
            opts = (aspect_ratio=1, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), clims=(0.0, 1.0), c=:davos, xlabel="Lx", ylabel="Ly", title="time = $(round(it*dt, sigdigits=3))")
            display(heatmap(xc, yc, Array(C)'; opts...))
        end
    end
    t_toc = Base.time() - t_tic
    A_eff = 2/1e9*nx*ny*sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it  = t_toc/niter                  # Execution time per iteration [s]
    T_eff = A_eff/t_it                   # Effective memory throughput [GB/s]
    @printf("Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits=3), niter)
    return
end

diffusion_2D(; do_visu=false)

max_stars_repo_name: luraess/amdgpu-test

max_issues_repo_name: luraess/amdgpu-test

max_forks_repo_name: luraess/amdgpu-test