text:
@testset "feature_transform" begin
    function ind2cart(F)
        s = CartesianIndices(axes(F))
        return map(i -> CartesianIndex(s[i]), F)
    end
    @testset "Square Images" begin
        # (1)
        A = [true false; false true]
        F = feature_transform(A)
        @test F == ind2cart([1 1; 1 4])
        D = distance_transform(F)
        @test D == [0 1; 1 0]

        # (2)
        A = [true true; false true]
        F = feature_transform(A)
        @test F == ind2cart([1 3; 1 4])
        D = distance_transform(F)
        @test D == [0 0; 1 0]

        # (3)
        A = [false false; false true]
        F = feature_transform(A)
        @test F == ind2cart([4 4; 4 4])
        D = distance_transform(F)
        @test D ≈ [sqrt(2) 1.0; 1.0 0.0]

        # (4)
        A = [true false true; false true false; true true false]
        F = feature_transform(A)
        @test F == ind2cart([1 1 7; 1 5 5; 3 6 6])
        D = distance_transform(F)
        @test D == [0 1 0; 1 0 1; 0 0 1]

        # (5)
        A = [false false true; true true false; true true true]
        F = feature_transform(A)
        @test F == ind2cart([2 5 7; 2 5 5; 3 6 9])
        D = distance_transform(F)
        @test D == [1 1 0; 0 0 1; 0 0 0]

        # (6)
        A = [
            true false true true
            false true false false
            false true true false
            true false false false
        ]
        F = feature_transform(A)
        @test F == ind2cart([1 1 9 13; 1 6 6 13; 4 7 11 11; 4 4 11 11])
        D = distance_transform(F)
        @test D ≈ [0.0 1.0 0.0 0.0; 1.0 0.0 1.0 1.0; 1.0 0.0 0.0 1.0; 0.0 1.0 1.0 sqrt(2)]
    end

    @testset "Rectangular Images" begin
        # (1)
        A = [true false true; false true false]
        F = feature_transform(A)
        @test F == ind2cart([1 1 5; 1 4 4])
        D = distance_transform(F)
        @test D == [0 1 0; 1 0 1]

        # (2)
        A = [true false; false false; false true]
        F = feature_transform(A)
        @test F == ind2cart([1 1; 1 6; 6 6])
        D = distance_transform(F)
        @test D == [0 1; 1 1; 1 0]

        # (3)
        A = [
            true false false
            true false false
            false true true
            true true true
            false true false
        ]
        F = feature_transform(A)
        @test F == ind2cart([1 1 1; 2 2 13; 2 8 13; 4 9 14; 4 10 10])
        D = distance_transform(F)
        @test D == [0.0 1.0 2.0; 0.0 1.0 1.0; 1.0 0.0 0.0; 0.0 0.0 0.0; 1.0 0.0 1.0]

        @test feature_transform(Gray.(A)) == F
    end

    @testset "Corner Case Images" begin
        null1 = CartesianIndex((typemin(Int),))
        null2 = CartesianIndex((typemin(Int), typemin(Int)))
        # (1)
        A = [false]
        F = feature_transform(A)
        @test F == [null1]
        D = distance_transform(F)
        @test D == [Inf]

        # (2)
        A = [true]
        F = feature_transform(A)
        @test F == ind2cart([1])
        D = distance_transform(F)
        @test D == [0]

        # (3)
        A = [true false]
        F = feature_transform(A)
        @test F == ind2cart([1 1])
        D = distance_transform(F)
        @test D == [0 1]

        # (4)
        A = [false; false]
        F = feature_transform(A)
        @test F == [null1; null1]
        D = distance_transform(F)
        @test D == [Inf; Inf]

        # (5)
        A = [true; true]
        F = feature_transform(A)
        @test F == ind2cart([1; 2])
        D = distance_transform(F)
        @test D == [0; 0]

        # (6)
        A = [true; true; false]
        F = feature_transform(A)
        @test F == ind2cart([1; 2; 2])
        D = distance_transform(F)
        @test D == [0; 0; 1]

        # (7)
        A = falses(3, 3)
        F = feature_transform(A)
        @test all(x -> x == null2, F)
        D = distance_transform(F)
        @test all(x -> x == Inf, D)

        # (8)
        A = trues(4, 2, 3)
        F = feature_transform(A)
        @test F == ind2cart(reshape(1:length(A), size(A)))
        D = distance_transform(F)
        @test all(x -> x == 0, D)

        # (9)
        A = falses(4, 2, 3)
        A[3, 1, 2] = true
        @test all(==(CartesianIndex(3, 1, 2)), feature_transform(A))
    end

    @testset "Anisotropic images" begin
        A, w = [false false; false true], (3, 1)
        F = feature_transform(A; weights=w)
        @test F == ind2cart([4 4; 4 4])
        D = distance_transform(F, w)
        @test D ≈ [sqrt(10) 3.0; 1.0 0.0]

        A, w = [false false; false true], (1, 3)
        F = feature_transform(A; weights=w)
        @test F == ind2cart([4 4; 4 4])
        D = distance_transform(F, w)
        @test D ≈ [sqrt(10) 1.0; 3.0 0.0]

        A, w = [true false; false true], (3, 1)
        F = feature_transform(A; weights=w)
        @test F == ind2cart([1 1; 4 4])
        D = distance_transform(F, w)
        @test D == [0 1; 1 0]

        A, w = [true false; false true], (1, 3)
        F = feature_transform(A; weights=w)
        @test F == ind2cart([1 4; 1 4])
        D = distance_transform(F, w)
        @test D == [0 1; 1 0]
    end
end

nothing


text:
module Ccp

import ..Clarus

function defaultresources(;params...)
  return Clarus.api_request("CCP","DefaultResources";params...)
end

end #MODULE-END

ext: jl


text:
mutable struct Evaluation
    site :: Vector{Float64}
    result :: Union{Vector{Float64}, Nothing}

    Evaluation( site ) = new( site, nothing );
end

@with_kw struct DataBase <: SomeDataBase
	eval_dict :: Dict{ Int, Evaluation } = Dict{ Int, Evaluation }();
end

# DataBase methods
ids( db :: DataBase ) = collect( keys( db.eval_dict ) );
sites( db :: DataBase, ids :: Vector{Int} ) = [ db.eval_dict[id].sites for id in ids ];
function results( db :: DataBase, ids :: Vector{Int}; include_nothing = true )
    all_results = [ db.eval_dict[id].result for id in ids ];
    if !include_nothing
        deleteat!(all_results, isnothing.(all_results))
    end
    return all_results
end

"Push a single `site` into `db`, set Evaluation result to nothing."
function push_site!(db:: DataBase, site :: Vector{Float64}, id :: Int = -1 )
    if id < 0
        max_id = let db_ids = ids(db);
            isempty(db_ids) ? 0 : maximum( db_ids )
        end;
        id = max_id + 1;
    end

    db.eval_dict[id] = Evaluation(site);
end

"Push multiple `sites` into `db`, set Evaluation results to nothing."
function push_sites!(db:: DataBase, sites :: Vector{Vector{Float64}})
    max_id = let db_ids = ids(db);
        isempty(db_ids) ? 0 : maximum( db_ids )
    end;
    new_ids = Vector{Int}(undef, length(sites) )
    for (site_ind, site) ∈ enumerate(sites)
        id = max_id + site_ind;
        new_ids[ site_ind ] = id;
        push_site!(db, site, id)
    end
    return new_ids
end

function eval!( db :: DataBase, f :: F where F<:Function, id :: Int; reevaluate = true )
    if reevaluate || isnothing( db.eval_dict[id].result )
        db.eval_dict[id].result = f( db.eval_dict[id].site );
    end
    nothing
end


text:
# jinc.jl

using MIRT: jinc

using Plots: plot
using Test: @test, @inferred

@test jinc(0) ≈ pi/4

jinx = x -> jinc.(x)
r = LinRange(-10,10,201)
@inferred jinx(r)
for T in (Int, Float16, Float32, Float64)
	@inferred jinc(T(0))
	@inferred jinc(T(1))
end

r = LinRange(-10,10,201)
plot(r, jinc.(r))


text:
module TopicModelsVB

using DelimitedFiles
using SpecialFunctions
using LinearAlgebra
using Random
using Distributions
using OpenCL
using Crayons

export Document, Corpus
export TopicModel
export check_doc, check_corp, check_model
export LDA, fLDA, CTM, fCTM, CTPF, gpuLDA, gpuCTM, gpuCTPF
export readcorp, writecorp, abridge_corp!, alphabetize_corp!, compact_corp!, condense_corp!, pad_corp!, remove_empty_docs!, remove_redundant!, stop_corp!, trim_corp!, trim_docs!, fixcorp!, showdocs, showtitles, getvocab, getusers
export train!
export @gpu
export gendoc, gencorp, showtopics, showlibs, showdrecs, showurecs, predict, topicdist

include("macros.jl")
include("utils.jl")
include("Corpus.jl")
include("TopicModel.jl")
include("LDA.jl")
include("fLDA.jl")
include("CTM.jl")
include("fCTM.jl")
include("CTPF.jl")
include("gpuLDA.jl")
include("gpuCTM.jl")
include("gpuCTPF.jl")
include("modelutils.jl")

end


text:
"""
    with_changed_bound(args...; kwargs...)

Specifies a model variant that has a new bound set. Forwards arguments to
[`change_bound`](@ref). Intended for usage with [`screen`](@ref).
"""
with_changed_bound(args...; kwargs...) = m -> change_bound(m, args...; kwargs...)

"""
    with_changed_bounds(args...; kwargs...)

Specifies a model variant that has new bounds set. Forwards arguments to
[`change_bounds`](@ref). Intended for usage with [`screen`](@ref).
"""
with_changed_bounds(args...; kwargs...) = m -> change_bounds(m, args...; kwargs...)

"""
    with_removed_metabolite(args...; kwargs...)

Specifies a model variant without a certain metabolite. Forwards arguments to
[`remove_metabolite`](@ref). Intended to be used with [`screen`](@ref).
"""
with_removed_metabolite(args...; kwargs...) = m -> remove_metabolite(m, args...; kwargs...)

"""
    with_removed_metabolites(args...; kwargs...)

Plural version of [`with_removed_metabolite`](@ref), calls
[`remove_metabolites`](@ref) internally.
"""
with_removed_metabolites(args...; kwargs...) =
    m -> remove_metabolites(m, args...; kwargs...)

"""
    with_added_reactions(args...; kwargs...)

Specifies a model variant with reactions added. Forwards the arguments to
[`add_reactions`](@ref). Intended to be used with [`screen`](@ref).
"""
with_added_reactions(args...; kwargs...) = m -> add_reactions(m, args...; kwargs...)

"""
    with_removed_reaction(args...; kwargs...)

Specifies a model variant without a certain reaction. Forwards arguments to
[`remove_reaction`](@ref). Intended to be used with [`screen`](@ref).
"""
with_removed_reaction(args...; kwargs...) = m -> remove_reaction(m, args...; kwargs...)

"""
    with_removed_reactions(args...; kwargs...)

Plural version of [`with_removed_reaction`](@ref), calls
[`remove_reactions`](@ref) internally.
"""
with_removed_reactions(args...; kwargs...) = m -> remove_reactions(m, args...; kwargs...)


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Pkg.TOML
meta = length(ARGS) > 0 ? ARGS :
  keys(TOML.parsefile(joinpath(@__DIR__, "Notebooks.toml")))

convertjl = joinpath(@__DIR__, "convert.jl")

for proj in meta
  run(`$(Base.julia_cmd()) $convertjl $proj`)
end

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
function readMap()
    map = Dict()
    while !eof(stdin)
        line = split(readline(stdin), " -> ")
        map[line[1]] = line[2][1]
    end
    return map
end

function rec(left, right, d)
    memokey = (left, right, d)
    if memokey in keys(memo)
        return memo[memokey]
    end
    k = left*right
    if d == 0 || k ∉ keys(map)
        if left == right
            return Dict{Char, Int64}(left => 2)
        else
            return Dict{Char, Int64}(left => 1, right => 1)
        end
    end
    mid = map[k]
    l = rec(left, mid, d-1)
    r = rec(mid, right, d-1)
    result = merge(+, l, r)
    result[mid] -= 1
    memo[memokey] = result
    return result
end

function solve(p, n)
    cnt = Dict{Char, Int64}(p[1] => 1)
    for i in 2:length(p)
        cnt = merge(+, cnt, rec(p[i-1], p[i], n))
        cnt[p[i-1]] -= 1
    end
    return maximum(values(cnt)) - minimum(values(cnt))
end

poly = readline()
readline()
map = readMap()
memo = Dict{Tuple{Char, Char, Int64}, Dict{Char, Int64}}()
println(solve(poly, 10))
println(solve(poly, 40))

max_stars_repo_name: MKolman/advent-of-code

max_issues_repo_name: MKolman/advent-of-code

max_forks_repo_name: MKolman/advent-of-code


text:
# TODO:
# - 1. Improve argument handling
# - 2. Update code to use UnPack.jl to make it more readable
# - 3. Export new interface
# - 4. Document Interface

using Mustache, Highlights, .WeaveMarkdown, Markdown, Dates, Printf

const FORMATS = Dict{String,WeaveFormat}()

# TODO: do some assertion for necessary fields of `format`
register_format!(format_name::AbstractString, format::WeaveFormat) = push!(FORMATS, format_name => format)
register_format!(_, format) = error("Format needs to be a subtype of WeaveFormat.")

set_format_options!(doc; kwargs...) = set_format_options!(doc.format; kwargs...)

function render_doc(doc::WeaveDoc)
    restore_header!(doc)

    docformat = doc.format
    body = joinlines(render_chunk.(Ref(docformat), copy(doc.chunks)))
    return render_doc(docformat, body, doc)
end

include("common.jl")
include("htmlformats.jl")
include("texformats.jl")
include("pandocformats.jl")
include("miscformats.jl")


text:
#Matrices related to stochastic processes

abstract StochasticProcess
export evolve

immutable BrownianProcess <: StochasticProcess
end

# Generates Brownian motion
evolve(p::BrownianProcess, dx::Real, xend::Real) = cumsum(WhiteNoiseProcess(dx, xend))

# Generates a white noise process
function WhiteNoiseProcess(dx::Real, xend::Real)
  x=[0:dx:xend]
  randn(length(x),1)*sqrt(dx)
end

################
# Airy process #
################

immutable AiryProcess <: StochasticProcess
  beta::Real
end

# Calculates the largest eigenvalue of a stochastic Airy process
# with Brownian noise
function evolve(p::AiryProcess, dx::Real, xend::Real)
    x=[0:dx:xend]
    N=length(x)

    #Discretized Airy operator
    a=-(2/dx^2)*ones(N) - x
    b=+(1/dx^2)*ones(N-1)
    #Plus noise
    dW=WhiteNoiseProcess(dx, xend)
    a+=(2/sqrt(p.beta))*dW/dx

    maxeig(SymTridiagonal(a,b))
end

#Sample program
#t=10000 #number of trials
#v=[AiryProcess(0.001, 10, 2) for i=1:t]
#binsize=.2
#grid=[-5:binsize:2]
#x=hist(v, grid)
#for i=1:length(grid)
#    @printf("%10.5f  %8d  %s\n",grid[i],x[i],repeat("*",int(x[i]/(t*binsize)*200)))
#end


text:
"""
    nn_interpolation(...)

* array_in is the array/map, which will be interpolated in a higher resolution.
* nlongitude is the new width of the mesh, nlatitude will be calculated as nlongitude/2 + 1
"""
function nn_interpolation(array_in, nlongitude=256)

     (width_in, height_in) = size(array_in)
     width_out = nlongitude
     height_out = convert(Int64, width_out/2 + 1)
     array_out = fill(0, (width_out, height_out))

     # calculate the new distance between each dot/pixle (ratio)
     # "-1" because indicies starts by 1 in Julia
     ratio_w = (width_in-1)/(width_out-1)
     ratio_h = (height_in-1)/(height_out-1)

     row = ones(Float64, width_out)
     for i = 1 : length(row)
          row[i] = row[i] + (i-1) * ratio_w
     end
     row[end] = round(width_in)

     col = ones(Float64, height_out)
     for i = 1 : length(col)
          col[i] = col[i] + (i-1) * ratio_h
     end
     col[end] = round(height_in)

     width_positions = round.(Int64, row)
     # col .+= 0.0000001  -> correction for function round, because round(2.5) equals 2 ???!!!
     height_positions = round.(Int64, col)

     wide = 1
     high = 1

     for b in width_positions
          high = 1
          for h in height_positions
               array_out[wide, high] = array_in[b, h]
               high += 1

          end
          wide += 1
     end
     return array_out
end

"""
    upscale_world(...)

* scales up the discrete world map in width nlongitude and height nlongitude/2 + 1
* saves upscaled world in the input array's directory (dirpath)
"""
function upscale_world(dirpath = joinpath(@__DIR__, "..","input","world"), filename = "The_World128x65.dat",nlongitude=256)

     world = Klimakoffer.read_geography(joinpath(dirpath, filename), 128, 65)
     upscaled_map = nn_interpolation(world, nlongitude)
     (nlongitude, nlatitude) = size(upscaled_map)

     new_fname = string("The_World", nlongitude, "x", nlatitude, ".dat")
     new_fp = joinpath(@__DIR__, dirpath, new_fname)

     if !isfile(new_fp)
          touch(new_fp)
     end

     open(new_fp,"w") do file
          for lat = 1:nlatitude
               for long = 1:nlongitude
                    write(file, string(upscaled_map[long,lat]))
               end
               write(file, "\n")
           end
     end
end

"""
    bilinear_interpolation(...)

* array_in is the (albedo-)array, which will be interpolated in a higher resolution.
* nlongitude is the new width of the mesh, nlatitude will be calculated as nlongitude/2 + 1
"""
function bilinear_interpolation(array_in, nlongitude =256)

     (width_in, height_in) = size(array_in)

     width_out = nlongitude
     height_out = convert(Int64, (width_out/2 + 1))
     array_out = Array{Float64}(ones((width_out, height_out)))

     # calculate the new distance between each dot/pixle (ratio)
     # "-1" because the indicies in Julia
     ratio_w = (width_in-1)/(width_out-1)
     ratio_h = (height_in-1)/(height_out-1)

     high = 1
     wide = 1

     for w = 1 : width_out
          high = 1
          for h = 1 : height_out

               # the coordinates from the current pixel/knot
               y = 1 + (h-1) * ratio_h
               x = 1 + (w-1) * ratio_w

               # calculate coordinates from surrounding pixels/knots
               x_floor = floor(Int64, x)
               y_floor = floor(Int64, y)
               x_ceil = min(ceil(Int64, x), width_in)
               y_ceil = min(ceil(Int64, y), height_in)

               # get values of the surrounding pixels
               # case 1: x and y are intergers => x_ceil = x_floor, y_ceil = y_floor
               if x_ceil == x_floor && y_ceil == y_floor
                    val = array_in[x_ceil , y_ceil]

               # case 2: x is an integer
               elseif x_ceil == x_floor
                    val1 = array_in[x_ceil, y_floor]
                    val2 = array_in[x_ceil, y_ceil]
                    val = val1 * (y_ceil - y) + val2 * (y- y_floor)

               #case 3: y is an integer
               elseif y_ceil == y_floor
                    val1 = array_in[x_floor, y_ceil]
                    val2 = array_in[x_ceil, y_ceil]
                    val = val1 * (x_ceil - x) + val2 * (x - x_floor)
               # case 4: x and y are floating point numbers
               else
                    v1 = array_in[x_floor, y_floor]
                    v2 = array_in[x_ceil, y_floor]
                    v3 = array_in[x_floor, y_ceil]
                    v4 = array_in[x_ceil, y_ceil]

                    # calculate the pixels (x, y) value
                    val1 = (v1 * (x_ceil - x) + v2 * (x - x_floor))  # first linear interpolation in x direction between v1 and v2
                    val2 = (v3 * (x_ceil - x) + v4 * (x - x_floor))  # second linear interpolation in x direction between v3 and v4
                    val = val1 * (y_ceil - y) + val2 * (y- y_floor)  # final "bilinear interpolation"
               end

               array_out[wide, high] = round(val, digits=2)
               high += 1
          end
          wide += 1
     end
     return array_out
end

"""
    upscale_albedo(...)

* scales up the albedo map in width nlongitude and height nlongitude/2 + 1
* saves upscaled albedo in the input array's directory (dirpath)
"""
function upscale_albedo(dirpath = joinpath(@__DIR__,"..","input","albedo"), filename = "albedo128x65.dat", nlongitude=256)

     albedo = Klimakoffer.read_albedo(joinpath(@__DIR__,dirpath, filename), 128, 65)
     upscaled_albedo = bilinear_interpolation(albedo, nlongitude)
     (nlongitude, nlatitude) = size(upscaled_albedo)

     new_fname = string("albedo", nlongitude, "x", nlatitude, ".dat")
     new_fp = joinpath(@__DIR__,dirpath, new_fname)
     if !isfile(new_fp)
          touch(new_fp)
     end

     open(new_fp,"w") do file
          for lat = 1:nlatitude
               for long = 1:nlongitude
                    write(file, string("      ",upscaled_albedo[long,lat]))
               end
               write(file,"\n")
           end
     end
end

"""
    save_outline_from_path(...)

* saves the outline from discrete world, which is load from the given filepath (=dirpath+filename)
"""
function save_outline_from_path(dirpath = joinpath(@__DIR__, "..", "input", "world"), filename = "The_World128x65.dat",nlongitude=128)

     nlatitude = convert(Int64, nlongitude/2+1)
     world = Klimakoffer.read_geography(joinpath(dirpath, filename), nlongitude, nlatitude)
     outline = get_outline_from_world(world)

     filename = string("The_World_Outline", nlongitude, "x", nlatitude, ".dat")
     new_fp = joinpath(@__DIR__,dirpath, filename)
     if !isfile(new_fp)
          touch(new_fp)
     end

     open(new_fp,"w") do file
          for lat = 1:nlatitude
               for long = 1:nlongitude
                    write(file, string(outline[long,lat]))
               end
               write(file, "\n")
          end
     end
end

"""
    get_outline_from_world(world)

* get the outline from our discrete world, so coast is presented by ones
"""
function get_outline_from_world(world)
     (nlongitude, nlatitude) = size(world)
     world = clear_map(world)
     outline = Array{Int64}(zeros((nlongitude, nlatitude)))

     for b = 1:nlongitude
          for h = 1:nlatitude
               if world[b,h] == 1 || world[b,h] == 3 # 1 = land, 3 = permanent snow cover

                    if b == 1 && h == 1
                         if !(world[b+1,h] == 1 && world[b,h+1] == 1)
                              outline[b,h] = 1
                         end

                    elseif b == 1 && h != 1 && h != nlatitude
                         if world[b,h+1] == 0 || world[b,h-1] == 0 || world[b+1,h] == 0
                              outline[b,h] = 1
                         end

                    elseif b == 1 && h == nlatitude
                         if !(world[b+1,h] == 1 && world[b,h-1] == 1)
                              outline[b,h] = 1
                         end

                    elseif h == nlatitude && b != 1 && b != nlongitude
                         if world[b-1,h] == 0 || world[b, h-1] == 0 || world[b+1,h] == 0
                              outline[b,h] = 1
                         end

                    elseif h == nlatitude && b == nlongitude
                         if !(world[b-1,h] == 1 && world[b,h-1] == 1)
                              outline[b,h] = 1
                         end

                    elseif b == nlongitude && h != 1 && h != nlatitude
                         if  world[b-1,h] == 0 || world[b,h-1] == 0 || world[b,h+1] == 0
                              outline[b,h] = 1
                         end

                    elseif b == nlongitude && h == 1
                         if !(world[b-1,h] == 1 && world[b,h+1] == 1)
                              outline[b,h] = 1
                         end

                    elseif h == 1 && b != nlongitude && b != 1
                         if world[b-1,h] == 0 || world[b,h+1] == 0 || world[b+1,h] == 0
                              outline[b,h] = 1
                         end

                    else
                         if world[b-1,h] == 0 || world[b+1,h] == 0 || world[b,h-1] == 0 || world[b,h+1] == 0
                              outline[b,h] = 1
                         end

                    end
               end
          end
     end
     return outline # new
end

"""
    clear_map(map)

* a help function, which presents land (1) and permanent snow cover (3) by ones and
*    perennial sea ice (2), lakes, inland seas (4) and ocean (5) by zeros
"""
function clear_map(map)
     (nlongitude, nlatitude) = size(map)
     for b =1:nlongitude
          for h=1:nlatitude
               if map[b,h] == 3
                    map[b,h] = 1

               elseif map[b,h] == 2 || map[b,h] >= 4
                    map[b,h] = 0
               end
          end
     end
     return map
end


text:
julia> data = {["a", "b", "c"], ["", "q", "z"], ["zap", "zip", "Zot"]}
3-element Array{Any,1}:
 ["a","b","c"]
 ["","q","z"]
 ["zap","zip","Zot"]

julia> sorttable(data, column=2, reverse=true) # named arguments
3-element Array{Any,1}:
 ["zap","zip","Zot"]
 ["","q","z"]
 ["a","b","c"]

julia> sorttable(data, >, 2) # the same thing, with positional arguments
3-element Array{Any,1}:
 ["zap","zip","Zot"]
 ["","q","z"]
 ["a","b","c"]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
	HighPass(type::Symbol, f, dBperOct, fs=96e3)

Return Biquad of sampling frequency fs as DSP.SecondOrderSections that resembles the Highpass
functions of the Lake EQ.

`type` must be one of [:BT, :BS, :LR] resembling Butterworth, Bessel and
Linkwitz-Riley respectively. `dBperOct` must be one of [6, 12, 18, 24, 30, 36, 42, 48].

# Implementation

## Butterworth

The design for the Butterworth filters uses the implementation by DSP.jl.

## Linkwitz-Riley

The Linkwitz-Riley filters are constructed from two Butterworth stages
at the same frequency.

## Bessel

The Bessel filters are constructed from the Polynomial Representation given in:

https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119011866.app1 (12.6.2019)

(Appendix 1 of *Richard W. Middlestead*:
Digital Communications with Emphasis on Data Modems: Theory, Analysis, Design, Simulation, Testing, and Applications.
1st edition, Published 2017 by John Wiley & Sons, Inc.)

The resulting filter is normalized to be -3 dB at 1 rad/s. The -3 dB frequencies
for all orders of Bessel filters in the Lake processor at a frequency of 1 kHz
were measured and are used to shift the Bessel prototype to fit the Lake
filter after adjustment to the given frequency and filter type.

The transformation from lowpass prototype to highpass is then achieved by simply adding
n zeros to the filter at 0+0im where n is the filter order.

This is contrary to the standard method of converting a lowpass to
highpass where the substitution s -> 1/s is performed. However, while
this transformation does accurately mirror the Level response, the group delay
of the newly created highpass is not flat anymore but peaked above orders of 2.
This peak gets stronger the higher the order of the filter.

Lake provides flat group delay response in Bessel highpass filters by adding zeros
in the nominator, but not reversing the coefficient order in the denominator, both
of which would happen in the substitution s -> 1/s. The Lake version is still
approximately equal to the subtitution method because Bessel polynomials are approximately
reversible with the right normalization for s. (see https://www.rane.com/note147.html (22.6.2019))

Note that in contrast to the Lake Bessel lowpass filters, odd orders are not phase-inverted.
"""
function HighPass(type::Symbol, f, dBperOct, fs=96e3)
	responsetype = Highpass(f, fs=fs)
	return PassFilter(responsetype, type, dBperOct)
end

"""
	LowPass(type::Symbol, f, dBperOct, fs=96e3)

Return Biquad of sampling frequency fs as DSP.SecondOrderSections that resembles the Lowpass
functions of the Lake EQ.

`type` must be one of [:BT, :BS, :LR] resembling Butterworth, Bessel and
Linkwitz-Riley respectively. `dBperOct` must be one of [6, 12, 18, 24, 30, 36, 42, 48].

# Implementation

## Butterworth

The design for the Butterworth filters uses the implementation by DSP.jl.

## Linkwitz-Riley

The Linkwitz-Riley filters are constructed from two Butterworth stages
at the same frequency.

## Bessel

The Bessel filters are constructed from the Polynomial Representation given in:

https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119011866.app1 (12.6.2019)

(Appendix 1 of *Richard W. Middlestead*:
Digital Communications with Emphasis on Data Modems: Theory, Analysis, Design, Simulation, Testing, and Applications.
1st edition, Published 2017 by John Wiley & Sons, Inc.)

The resulting filter is normalized to be -3 dB at 1 rad/s. The -3 dB frequencies
for all orders of Bessel filters in the Lake processor at a frequency of 1 kHz
were measured and are used to shift the Bessel prototype to fit the Lake
filter after adjustment to the given frequency and filter type.

Note that for odd orders of Bessel filters, the polarity is reversed.
"""
function LowPass(type::Symbol, f, dBperOct, fs=96e3)
	responsetype = Lowpass(f, fs=fs)
	return PassFilter(responsetype, type, dBperOct)
end

"""
	Bessel(N::Int)

Return analog Bessel lowpass prototype as DSP.ZeroPoleGain that fits Lake EQ
after transformation.
"""
function Bessel(N::Int)
	b(n::Int) = factorial(2*N-n)//(2^(N-n)*factorial(N-n)*factorial(n))
	denominator = [b(n)|>Float64 for n=N:-1:0]
	nominator = zeros(N+1)
	nominator[end] = denominator[end]

	if isodd(N)
		nominator = -nominator
	end

	myfilter = PolynomialRatio(nominator, denominator)
	myfilter2 = convert(ZeroPoleGain, myfilter)

	#construct function for NLsolve to find -3dB point
	myfilterdB(w) = 20*log10.(abs.(freqs(myfilter2, w)))
	function func!(dB, w)
	   dB[1] = myfilterdB(w[1]).+3
	end

	res = nlsolve(func!, [3.], autodiff=:forward)
	w3dB = res.zero[1]

	BesselCorrectionFactor = 1/w3dB

	LakeCorrectionFactor = ([997.2, 784.8, 710.4, 659.2, 615.2, 577.5, 544.6, 516.7]./1000)[N]

	correctionFactor = LakeCorrectionFactor*BesselCorrectionFactor

	myfilter3 = ZeroPoleGain(myfilter2.z .* correctionFactor,
							 myfilter2.p .* correctionFactor,
				myfilter2.k * correctionFactor^(length(myfilter2.p)-length(myfilter2.z)))
	return myfilter3
end

function BesselInvert(prototype)
	N = length(prototype.p)
	if isodd(N)
		newGain = -prototype.k
	else
		newGain = prototype.k
	end
	return ZeroPoleGain(fill(0+0im, N), prototype.p, newGain)
end

function PassFilter(responsetype, type, dBperOct)
	if dBperOct in [6,12,18,24,30,36,42,48]
		Poles = dBperOct/6|>Int
	else
		error("dBperOct is not one of [6,12,18,24,30,36,42,48]")
	end

	if type == :BT
		designmethod = Butterworth(Poles)
		Filter = convert(SecondOrderSections, digitalfilter(responsetype, designmethod))
	elseif type == :BS
		designmethod = Bessel(Poles)
		if typeof(responsetype)<:Lowpass
			Filter = convert(SecondOrderSections, digitalfilter(responsetype, designmethod))
		else #Highpass
			hp_responsetype = Lowpass(responsetype.w)
			hpbs_designmethod = BesselInvert(designmethod)
			Filter = convert(SecondOrderSections, digitalfilter(hp_responsetype, hpbs_designmethod))
		end
	elseif type == :LR
		if Poles in [2,4,6,8]
			designmethod = Butterworth(Poles/2|>Int)
			OneStage = convert(ZeroPoleGain, digitalfilter(responsetype, designmethod))
			TwoStages = OneStage*OneStage
			Filter = SecondOrderSections(TwoStages)
		else
			error("Linkwitz-Riley must have dBperOct of 12, 24, 36 or 48")
		end
	else
		error("unknown Low- or Highpass type")
	end
	return Filter
end


text:
include("FractalDimensions.jl")

#S_0 = [[0.0,0.0], [1.0,0.0], [0.5,1.0]]
S_0 = [[0.0,0.0], [1.0,0.0], [0.5,2.0]]
#S_0 = [[0.0,0.0], [1.0,0.0], [0.0,1.0]]

f1(x) = (S_0[1]+x)/2
f2(x) = (S_0[2]+x)/2
f3(x) = (S_0[3]+x)/2

# WARNING This will use up your memory and slow down your computer for k>10
function sierpinsky_endpoints(k)
        s::Vector{Vector{Float64}} = []
        s1 = S_0
        for i in 1:k
            s2 = setdiff([f1.(s1); f2.(s1); f3.(s1)], s1)
            append!(s, s1)
            s1 = s2
        end
    return s
end

s = sierpinsky_endpoints(10)

using Plots
x = [u[1] for u in s]
y = [u[2] for u in s]

plotlyjs()
scatter(x,y, markersize=1, markerstrokealpha=0.0, legend=false)
savefig("sierpinsky.eps")

d, N, ε = box_counting_dimension(s, 9, 0.0)

loglog_regression(N,ε)

log(3)/log(2)

plot(-log.(ε), log.(N), xlabel="-ln ε", ylabel="ln N", marker=:d, legend=false)
savefig("boxdim_sierpinsky.eps")

C, r = correlation_dimension(s)

loglog_regression(C,r)

plot(log.(r), log.(C), xlabel="ln r", ylabel="ln C", marker=:d, legend=false)
savefig("cordim_sierpinsky.eps")

max_stars_repo_name: csimal/chaos-and-fractals

max_issues_repo_name: csimal/chaos-and-fractals

max_forks_repo_name: csimal/chaos-and-fractals


text:
module LDAforCMB

using PyCall

export	trans_to_pico,
		custom_real_eigs,
		samp_cov

const pathtosrc = dirname(@__FILE__())
@pyimport pypico
const pico = pypico.load_pico(joinpath(pathtosrc,"chain_data/pico3_tailmonty_v33.dat"))
picoget(x::Dict{Symbol,Float64}) = pico[:get](;x...)["cl_TT"]
function trans_to_pico(lcdm_pars)
	# the lcdm_pars should be ordered as follows
	lcdm_names = [:tau,  :omch2,  :logA, :theta  , :ns,  :ombh2]
	lcdmD = Dict(zip(lcdm_names, lcdm_pars))
	picoinputs = Dict(
		:re_optical_depth => lcdmD[:tau],
		:massive_neutrinos => 3.046,
		symbol("scalar_amp(1)") => 1e-10*exp(lcdmD[:logA]),
		:helium_fraction => 0.248,
		symbol("scalar_nrun(1)") => 0.0,
		:theta => lcdmD[:theta],
		:omnuh2 => 0.0,
		:ombh2 => lcdmD[:ombh2],
		:omch2 => lcdmD[:omch2],
		symbol("scalar_spectral_index(1)") => lcdmD[:ns]
	)
	picoget(picoinputs)
end

function custom_real_eigs(M::Matrix)
	d, v = eig(M)
	perm = sortperm(real(d), rev=true)
	d = real(d[perm])
	v = real(v[:,perm])
	d,v
end
function custom_real_eigs(M::Matrix, nbasis)
	d, v, = eigs(M, nev = nbasis)
	perm = sortperm(real(d), rev=true)
	d = real(d[perm])
	v = real(v[:,perm])
	d, v
end

function samp_cov(S)
	d, v =  custom_real_eigs(S)
	sample = v * (sqrt(pos(d)) .* randn(size(S,1)))
end
function samp_cov(S, nb)
	d, v =  custom_real_eigs(S, nb)
	sample = v * (sqrt(pos(d)) .* randn(length(d)))
end

function pos(x)
	y=real(x)
	y[y.<0.0] = 0.0
	y
end

end # module

max_stars_repo_name: EthanAnderes/LDAforCMB

max_issues_repo_name: EthanAnderes/LDAforCMB

max_forks_repo_name: EthanAnderes/LDAforCMB


text:
using DamagedShearBand ; const DSB = DamagedShearBand
using MAT

# imposed values
ϵ̇₁₁ = -1e-4
σ₂₂ = -1e-6 # negative in compression

# initial higher compressive stress multiplier (S = σ₁₁/σ₂₂)
S = 1.0

# Create a Rheology type instance containing elastic moduli and damage parameters. Change default values by  supplying keywords arguments.
r = DSB.Rheology(D₀=0.3)

# initial damage
D_i = r.D₀

# time integration parameters :
tspan = (0.0,15) # initial and final simulation time
Δt_i = 0.1 # initial timestep

# Build stress and strain tensors
σᵢⱼ_i = DSB.build_principal_stress_tensor(r,S,σ₂₂,D_i ; abstol=1e-15) # takes care of the plane strain constraint by solving non linear out of plane strain wrt σ₃₃ using Newton algorithm
ϵᵢⱼ_i = DSB.compute_ϵij(r,D_i,σᵢⱼ_i)

# Integrate over time, by solving for the unknowns [σ₁₁next, σ₃₃next, ϵ₂₂next] at each timestep
# -- e₀ is the absolute tolerance used to adapt time stepping, may be a scalar or a NamedTuple with keys (D,σ,ϵ)
# -- abstol is the absolute tolerance on the unknowns in the Newton solver
t_vec, σᵢⱼ_vec, ϵᵢⱼ_vec, D_vec = DSB.adaptative_time_integration(r,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇₁₁,Δt_i,tspan ;
                                                                 abstol=1e-12,
                                                                 time_maxiter=nothing,
                                                                 newton_maxiter=100,
                                                                 e₀=(D=1e-4, σ=1.0, ϵ=1e-5),
                                                                 print_frequency=10000)

# make a dictionary out of time series :
vars_dict = Dict("sigma1"=>[σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec],
                 "sigma1"=>[σᵢⱼ[2,2] for σᵢⱼ in σᵢⱼ_vec],
                 "sigma1"=>[σᵢⱼ[3,3] for σᵢⱼ in σᵢⱼ_vec],
                 "eps1"=>[ϵᵢⱼ[1,1] for ϵᵢⱼ in ϵᵢⱼ_vec],
                 "eps2"=>[ϵᵢⱼ[2,2] for ϵᵢⱼ in ϵᵢⱼ_vec],
                 "eps2"=>[ϵᵢⱼ[3,3] for ϵᵢⱼ in ϵᵢⱼ_vec],
                 "D"=>D_vec,
                 "t"=>t_vec)

# write that dict into a .mat file :
matwrite("./examples/singlepoint_D0=$(r.D₀)_epsdot=$(ϵ̇₁₁).mat", vars_dict)


text:
#!/usr/bin/env julia

using Keldysh, Test, LinearAlgebra

@testset "branch" begin
  let tmax = 2.0, β = 3.0
    fwd = Branch(forward_branch, tmax)
    back = Branch(backward_branch, tmax)
    imag = Branch(imaginary_branch, β)

    @test fwd.domain == forward_branch
    @test back.domain == backward_branch
    @test imag.domain == imaginary_branch

    @test length(fwd) == tmax
    @test length(back) == tmax
    @test length(imag) == β

    @test fwd(0.0).val == 0.0
    @test fwd(1.0).val == tmax

    @test back(0.0).val == tmax
    @test back(1.0).val == 0.0

    @test imag(0.0).val == 0.0
    @test imag(1.0).val == -1.0im * β
  end
end

@testset "contour" begin
  let c = Contour(full_contour, tmax=2.0, β=5.0)
    for i in 1:3
      @test c.branches[i].domain == BranchEnum(i)
    end

    c = twist(c)
    for i in 1:3
      @test c.branches[i].domain == BranchEnum(mod1(i+1, 3))
    end

    @test nbranches(full_contour) == 3
    @test nbranches(keldysh_contour) == 2
    @test nbranches(imaginary_contour) == 1

    for b in instances(BranchEnum)
      @test get_branch(c, b).domain == b
    end
  end

  let tmax = 2.0, β = 5.0, npts_real = 21
    c = Contour(keldysh_contour, tmax=tmax)
    @test forward_branch ∈ c
    @test backward_branch ∈ c
    @test imaginary_branch ∉ c
  end
end

@testset "time_grid" begin
  let tmax = 2.0, β = 5.0, npts_real=21, npts_imag=51
    c = Contour(full_contour, tmax=tmax, β=β)
    grid = TimeGrid(c, npts_real=npts_real, npts_imag=npts_imag)
    @test grid.step[1] ≈ 0.1
    @test grid.step[2] ≈ -0.1
    @test grid.step[3] ≈ -0.1im

    @test length(grid, forward_branch) == npts_real
    @test length(grid, backward_branch) == npts_real
    @test length(grid, imaginary_branch) == npts_imag

    @test map(p -> p.idx, grid) == 1:(2npts_real + npts_imag)

    for i in 1:3
      @test grid.branch_bounds[i][1].val == grid.contour.branches[i](0.0)
      @test grid.branch_bounds[i][2].val == grid.contour.branches[i](1.0)
    end

    @test β == Keldysh.get_beta(grid, nothing)
    @test_throws AssertionError Keldysh.get_beta(grid, β)

    @test integrate(t -> 1, grid) ≈ -1.0im * β
  end

  let tmax = 2.0, β = 5.0, npts_real = 21
    c = Contour(keldysh_contour, tmax=tmax)
    grid = TimeGrid(c, npts_real=npts_real)
    @test grid.step[1] ≈ 0.1
    @test grid.step[2] ≈ -0.1
    @test length(grid, forward_branch) == npts_real
    @test length(grid, backward_branch) == npts_real

    @test map(p -> p.idx, grid) == 1:(2npts_real)

    for i in 1:2
      @test grid.branch_bounds[i][1].val == grid.contour.branches[i](0.0)
      @test grid.branch_bounds[i][2].val == grid.contour.branches[i](1.0)
    end

    @test_throws AssertionError Keldysh.get_beta(grid, nothing)
    @test β == Keldysh.get_beta(grid, β)
  end

  let tmax = 2.0, c = Contour(keldysh_contour, tmax=tmax)
    grid = TimeGrid(c, npts_real=51)

    Δt1 = TimeGF(grid) do t1, t2
      t1.val.val - t2.val.val
    end

    Δt2 = TimeGF(grid) do t1, t2
      integrate(t -> 1.0, grid, t1, t2)
    end
    @test Δt1 ≈ Δt2
  end

  let tmax = 2.0, c = twist(Contour(keldysh_contour, tmax=tmax))
    grid = TimeGrid(c, npts_real=51)

    Δt1 = TimeGF(grid) do t1, t2
      t1.val.val - t2.val.val
    end

    Δt2 = TimeGF(grid) do t1, t2
      integrate(t -> 1.0, grid, t1, t2)
    end
    @test Δt1 ≈ Δt2
  end

end

@testset "generate_gf" begin
  let tmax = 1.0, β = 1.0, ν = 1/1000, ϵ = 2.0
    c = twist(Contour(full_contour, tmax=tmax, β=β))
    grid = TimeGrid(c, npts_real=51, npts_imag=51)
    dos = Keldysh.gaussian_dos(ν=ν, ϵ=ϵ)

    hyb1 = dos2gf(dos, grid)
    hyb2 = gf_1level(grid, ϵ=ϵ)

    # gf_1level is gf for a delta function spectrum
    @test isapprox(hyb1, hyb2, atol=ν, norm=x -> norm(x, Inf))

    # can't also supply β if it is part of the contour
    @test_throws AssertionError dos2gf(dos, grid, β = β)
  end
end

@testset "dos" begin
  # Flat DOS
  let ν=5.0, D=2.0, dos = Keldysh.flat_dos(ν=ν, D=D)
    α = ν*D
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:3]
    moments_ref = (1/π)*[D*(1 + coth(α)),
                         0.,
                         D^3*(π^2 + α^2)*(1 + coth(α))/(3α^2),
                         0.]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Gaussian DOS
  let ν = 2.0, ϵ = 1.0, dos = Keldysh.gaussian_dos(ν=ν, ϵ=ϵ)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0,
                   ϵ,
                   ϵ^2 + 2ν,
                   ϵ^3 + 6ϵ*ν,
                   ϵ^4 + 12ϵ^2*ν + 12*ν^2,
                   ϵ^5 + 20ϵ^3*ν + 60ϵ*ν^2,
                   ϵ^6 + 30ϵ^4*ν + 180ϵ^2*ν^2 + 120ν^3
                  ]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Bethe DOS
  let t=2.0, dos = Keldysh.bethe_dos(t=t)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0, 0., t^2, 0., 2t^4, 0., 5t^6]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Linear chain DOS
  let t=2.0, dos = Keldysh.chain_dos(t=t)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0, 0., 2t^2, 0., 6t^4, 0., 20t^6]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Square lattice DOS
  let t=2.0, dos = Keldysh.square_dos(t=t)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0, 0., 4t^2, 0., 36t^4, 0., 400t^6]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end
end


text:
using OscillatoryIntegralsODE
using SpecialFunctions
using StaticArrays
using BenchmarkTools
using PyPlot
using OrdinaryDiffEq

ω = 100.
f(x) = exp(-x^2/16)

clf()
xs = 0:0.001:5.0
plot(xs, real.([f(x) * exp(1im * ω * x) for x in xs]), "-", label=raw"$e^{-x^2/16} e^{-i 100 x}$")
legend()
gcf()

##
hi = HarmonicIntegral{Complex{Float64}}(ω)
result = levintegrate(hi, f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)

ref = 0.003798145404262353+0.009953606082753418im
print("relative error = ", (result - ref) / abs(ref), "\n")

## timing
@btime levintegrate($si, $f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)

max_stars_repo_name: xzackli/OscillatoryIntegralsODE

max_issues_repo_name: xzackli/OscillatoryIntegralsODE

max_forks_repo_name: xzackli/OscillatoryIntegralsODE


text:
using Ipopt
using Base.Test

# First of all, test that hs071 example works
include("hs071_test.jl")

# Test setting some options
# String option
println("\nString option")
addOption(prob, "hessian_approximation", "exact")
@test_throws addOption(prob, "hessian_approximation", "badoption")
println("\nInt option")
# Int option
addOption(prob, "file_print_level", 3) == nothing
@test_throws addOption(prob, "file_print_level", -1)
# Double option
println("\nFloat option")
addOption(prob, "derivative_test_tol", 0.5)
@test_throws addOption(prob, "derivative_test_tol", -1.0)

# Test opening an output file
openOutputFile(prob, "blah.txt", 5)

# Test MathProgBase stuff
#include(joinpath(Pkg.dir("MathProgBase"),"test","linprog.jl"))
#linprogtest(IpoptSolver())
    using MathProgBase
    solver = IpoptSolver()
    sol = linprog([-1,0],[2 1],'<',1.5,solver)
    @test sol.status == :Optimal
    @test abs(sol.objval - -0.75) <= 1e-6
    @test norm(sol.sol - [0.75,0.0]) <= 1e-6

    sol = linprog([-1,0],sparse([2 1]),'<',1.5,solver)
    @test sol.status == :Optimal
    @test abs(sol.objval - -0.75) <= 1e-6
    @test norm(sol.sol - [0.75,0.0]) <= 1e-6

    # test infeasible problem:
    # min x
    # s.t. 2x+y <= -1
    # x,y >= 0
    sol = linprog([1,0],[2 1],'<',-1,solver)
#    @test sol.status == :Infeasible

    # test unbounded problem:
    # min -x-y
    # s.t. -x+2y <= 0
    # x,y >= 0
    sol = linprog([-1,-1],[-1 2],'<',[0],solver)
#    @test sol.status == :Unbounded

include(joinpath(Pkg.dir("MathProgBase"),"test","nlp.jl"))
nlptest(IpoptSolver())
convexnlptest(IpoptSolver())
rosenbrocktest(IpoptSolver())


text:
"Hi " * "there!"       # 字符串连接
"Ho " ^ 3              # 重复字符串
string("a = ", 123.3)  # 使用打印方式构造
repr(123.3)            # 返回show函数的输出
occursin("CD", "ABCD") # 检查第一个字符串是否包含于第二个中
"\"\n\t\$"             # C 语言风格的转义（escaping）字符串（注意 \$ 的转义）
x = 123
"$x + 3 = $(x+3)"      # 没有转义的 $ 用来原位求值（interpolation）
"\$199"                # 要表示 $ 符号，你必须转义
raw"D:\path"           # 原始（raw）字符串字面值，在表示 Windows 下的路径很有用
s = "abc"              # 类型为 String 的字符串
chop(s)                # 删除 s 的最后一个字符，返回一 SubString 类型对象
r = r"A|B"             # 创建新的正则表达式（语法同 Perl）
occursin(r, "CD")      # false，不匹配
m = match(r, "ACBD")   # 找到第一处正则表达式匹配

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
export prox_l1!

"""
simple scirpt to compute the l1 proximal map.
our rho is 1./rho compared to other definitions of the proximal map (just a convention)
"""

function prox_l1!(x::Vector{TF},rho::TF) where {TF<:Real}
  x .= sign.(x) .* max.(TF(0.0), abs.(x) .- (1 ./ rho))
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module CUSTATEVEC

using CUDA
using CUDA: CUstream, cudaDataType, @checked, HandleCache, with_workspace, libraryPropertyType
using CUDA: libcustatevec, unsafe_free!, @retry_reclaim, initialize_context

using CEnum: @cenum

const cudaDataType_t = cudaDataType

# core library
include("libcustatevec_common.jl")
include("error.jl")
include("libcustatevec.jl")
include("types.jl")
include("statevec.jl")

# cache for created, but unused handles
const idle_handles  = HandleCache{CuContext,custatevecHandle_t}()

function handle()
    cuda = CUDA.active_state()

    # every task maintains library state per device
    LibraryState = @NamedTuple{handle::custatevecHandle_t, stream::CuStream}
    states = get!(task_local_storage(), :CUQUANTUM) do
        Dict{CuContext,LibraryState}()
    end::Dict{CuContext,LibraryState}

    # get library state
    @noinline function new_state(cuda)
        new_handle = pop!(idle_handles, cuda.context) do
            handle = Ref{custatevecHandle_t}()
            custatevecCreate(handle)
            handle[]
        end

        finalizer(current_task()) do task
            push!(idle_handles, cuda.context, new_handle) do
                context!(cuda.context; skip_destroyed=true) do
                    custatevecDestroy(new_handle)
                end
            end
        end

        custatevecSetStream(new_handle, cuda.stream)

        (; handle=new_handle, cuda.stream)
    end
    state = get!(states, cuda.context) do
        new_state(cuda)
    end

    # update stream
    @noinline function update_stream(cuda, state)
        custatevecSetStream(state.handle, cuda.stream)
        (; state.handle, cuda.stream)
    end
    if state.stream != cuda.stream
        states[cuda.context] = state = update_stream(cuda, state)
    end

    return state.handle
end

end


text:
#*******************************************************/
#* Copyright(c) 2021 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use Knitro to solve the following
# simple nonlinear optimization problem.  This model is test problem
# HS40 from the Hock & Schittkowski collection.
#
# max   x0*x1*x2*x3        (obj)
# s.t.  x0^3 + x1^2 = 1    (c0)
#       x0^2*x3 - x2 = 0   (c1)
#       x3^2 - x1 = 0      (c2)
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO, Test

function example_nlp2noderivs(; verbose=true)
    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalFC                                      *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "obj" and "c" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalFC(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate nonlinear term in objective
        evalResult.obj[1] = x[1] * x[2] * x[3] * x[4]

        # Evaluate nonlinear terms in constraints
        evalResult.c[1] = x[1] * x[1] * x[1]
        evalResult.c[2] = x[1] * x[1] * x[4]

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Initialize Knitro with the problem definition.

    # Add the variables and specify initial values for them.
    # Note: any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    xIndices = KNITRO.KN_add_vars(kc, 4)
    for x in xIndices
        KNITRO.KN_set_var_primal_init_value(kc, x, 0.8)
    end

    # Add the constraints and set the rhs and coefficients
    KNITRO.KN_add_cons(kc, 3)
    KNITRO.KN_set_con_eqbnds_all(kc, [1.0, 0.0, 0.0])

    # Coefficients for 2 linear terms
    lconIndexCons = Int32[1, 2]
    lconIndexVars = Int32[2, 1]
    lconCoefs = [-1.0, -1.0]
    KNITRO.KN_add_con_linear_struct(kc, lconIndexCons, lconIndexVars, lconCoefs)

    # Coefficients for 2 quadratic terms

    # 1st term:  x1^2 term in c0
    # 2nd term:  x3^2 term in c2
    qconIndexCons = Int32[0, 2]
    qconIndexVars1 = Int32[1, 3]
    qconIndexVars2 = Int32[1, 3]
    qconCoefs = [1.0, 1.0]

    KNITRO.KN_add_con_quadratic_struct(
        kc,
        qconIndexCons,
        qconIndexVars1,
        qconIndexVars2,
        qconCoefs,
    )

    # Add callback to evaluate nonlinear(non-quadratic) terms in the model:
    #    x0*x1*x2*x3  in the objective
    #    x0^3         in first constraint c0
    #    x0^2*x3      in second constraint c1
    cb = KNITRO.KN_add_eval_callback(kc, true, Int32[0, 1], callbackEvalFC)

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MAXIMIZE)

    # Set option to println output after every iteration.
    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ITER : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    varbndInfeas, varintInfeas, varviols = KNITRO.KN_get_var_viols(kc, Cint[0, 1, 2, 3])
    coninfeas, conviols = KNITRO.KN_get_con_viols(kc, Cint[0, 1, 2])
    err = KNITRO.KN_get_presolve_error(kc)

    if verbose
        println()
        println("Knitro converged with final status = ", nStatus)
        # An example of obtaining solution information.
        println("  optimal objective value  = ", objSol)
        println("  optimal primal values x  = ", x)
        println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
        println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))
        println("Variables bound violations = ", varbndInfeas)
        println("Variables integrality violations = ", varintInfeas)
        println("Variables violation values = ", varviols)
        println("Constraints bound violations = ", coninfeas)
        println("Constraints violation values = ", conviols)
    end

    @testset "Example HS40 nlp1noderivs" begin
        @test varbndInfeas == [0, 0, 0, 0]
        @test varintInfeas == [0, 0, 0, 0]
        @test varviols ≈ [0.0, 0.0, 0.0, 0.0] atol = 1e-6
        @test coninfeas == [0, 0, 0]
        @test conviols ≈ [0.0, 0.0, 0.0] atol = 1e-6
        @test KNITRO.KN_get_abs_feas_error(kc) == max(conviols...)
    end

    # Delete the Knitro solver instance.
    return KNITRO.KN_free(kc)
end

if KNITRO.KNITRO_VERSION >= v"12.4"
    example_nlp2noderivs(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)
else
    println("Example `nlp2noderivs.jl` is only available with Knitro >= 12.4")
end


text:
using FFTW

@testset "functions" begin
    ag = rand(Gray{Float32}, 4, 5)
    ac = rand(RGB{Float32}, 4, 5)
    for (f, args) in ((fft, (ag,)), (fft, (ag, 1:2)), (plan_fft, (ag,)),
                      (rfft, (ag,)), (rfft, (ag, 1:2)), (plan_rfft, (ag,)),
                      (fft, (ac,)), (fft, (ac, 1:2)), (plan_fft, (ac,)),
                      (rfft, (ac,)), (rfft, (ac, 1:2)), (plan_rfft, (ac,)))
        ret = @test_throws ErrorException f(args...)
        @test occursin("channelview", ret.value.msg)
        @test occursin(eltype(args[1])<:Gray ? "1:2" : "2:3", ret.value.msg)
    end
    for (a, dims) in ((ag, 1:2), (ac, 2:3))
        @test ifft(fft(channelview(a), dims), dims) ≈ channelview(a)
        ret = @test_throws ErrorException rfft(a)
        @test occursin("channelview", ret.value.msg)
        @test occursin("$dims", ret.value.msg)
        @test irfft(rfft(channelview(a), dims), 4, dims) ≈ channelview(a)
    end

    a = [RGB(1,0,0), RGB(0,1,0), RGB(0,0,1)]
    @test a' == [RGB(1,0,0) RGB(0,1,0) RGB(0,0,1)]

    a = [RGB(1,0,0) RGB(0,1,0); RGB(0,0,1) RGB(0, 0, 0)]
    @test a' == [RGB(1,0,0) RGB(0,0,1); RGB(0,1,0) RGB(0, 0, 0)]
end

nothing


text:
######################################################################################################
# Copy from https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/ by Wouter Koolen  #
######################################################################################################
# Tracking helper types

struct CTracking
    sumw;
    CTracking(N) = new(
        Float64.(N)
    ); # makes a copy of the starting situation
end
abbrev(_::Type{CTracking}) = "C";
# add a weight vector and track it
function track(t::CTracking, N, w)
    @assert all(N .≤ t.sumw.+1) "N $N  sumw $(t.sumw)";
    @assert all(N .≥ t.sumw.-log(length(N))) "N $N  sumw $(t.sumw)"; # This line is removed in https://github.com/xuedong/LinBAI.jl
    @assert sum(N) ≈ sum(t.sumw);
    t.sumw .+= w;
    argmin(N .- t.sumw);
end

struct DTracking
    DTracking(N) = new();
end
abbrev(_::Type{DTracking}) = "D";
function track(t::DTracking, N, w)
    argmin(N .- sum(N).*w);
end

# Wrapper to add forced exploration to a tracking rule
struct ForcedExploration
    t;
end
function track(fe::ForcedExploration, N, w)
    t = sum(N);
    K = length(N);
    undersampled = N .≤ sqrt(t) .- K/2;
    if any(undersampled)
        track(fe.t, N, undersampled/sum(undersampled));
    else
        track(fe.t, N, w);
    end
end