text:
using ImageCore, Colors, FixedPointNumbers, OffsetArrays
using Base.Test

@testset "reinterpret" begin
    a = rand(Gray{U8}, (4,5))
    for T in (Gray{U8}, Gray{Float32}, Gray{Float64})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpret(eltype(T), b))
        if ImageCore.squeeze1
            @test isa(rb, Array{eltype(T),2})
            @test size(rb) == (4,5)
        else
            @test isa(rb, Array{eltype(T),3})
            @test size(rb) == (1,4,5)
        end
        c = copy(rb)
        rc = @inferred(reinterpret(T, c))
        @test isa(rc, Array{T,2})
        @test size(rc) == (4,5)
    end
    a = rand(AGray{U8}, (4,5))
    for T in (AGray{U8}, GrayA{Float32}, AGray{Float64})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpret(eltype(T), b))
        @test isa(rb, Array{eltype(T),3})
        @test size(rb) == (2,4,5)
        c = copy(rb)
        rc = @inferred(reinterpret(T, c))
        @test isa(rc, Array{T,2})
        @test size(rc) == (4,5)
    end
    a = rand(RGB{U8}, (4,5))
    for T in (RGB{U8}, HSV{Float32}, XYZ{Float64})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpret(eltype(T), b))
        @test isa(rb, Array{eltype(T),3})
        @test size(rb) == (3,4,5)
        c = copy(rb)
        rc = @inferred(reinterpret(T, c))
        @test isa(rc, Array{T,2})
        @test size(rc) == (4,5)
    end
    for a in (rand(RGB{U8}, 4), rand(RGB{U8}, (4,5)))
        b = @inferred(reinterpret(HSV{Float32}, float32.(a)))
        @test isa(b, Array{HSV{Float32}})
        @test ndims(b) == ndims(a)
    end
    a = rand(ARGB{U8}, (4,5))
    for T in (ARGB{U8}, AHSV{Float32}, AXYZ{Float64})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpret(eltype(T), b))
        @test isa(rb, Array{eltype(T),3})
        @test size(rb) == (4,4,5)
        c = copy(rb)
        rc = @inferred(reinterpret(T, c))
        @test isa(rc, Array{T,2})
        @test size(rc) == (4,5)
    end
    # RGB1/RGB4
    a = rand(RGB{U8}, (4,5))
    for T in (RGB1{U8},RGB4{Float32})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpret(eltype(T), b))
        @test isa(rb, Array{eltype(T),3})
        @test size(rb) == (4,4,5)
        c = copy(rb)
        rc = @inferred(reinterpret(T, c))
        @test isa(rc, Array{T,2})
        @test size(rc) == (4,5)
    end
    a = [RGB(1,0,0) RGB(0,0,1);
         RGB(0,1,0) RGB(1,1,1)]
    @test reinterpret(U8, a) == cat(3, [1 0; 0 1; 0 0], [0 1; 0 1; 1 1])
    b = convert(Array{BGR{U8}}, a)
    @test reinterpret(U8, b) == cat(3, [0 0; 0 1; 1 0], [1 1; 0 1; 0 1])
    # indeterminate type tests
    a = Array(RGB{AbstractFloat},3)
    @test_throws ErrorException reinterpret(Float64, a)
    Tu = TypeVar(:T)
    a = Array(RGB{Tu},3)
    @test_throws ErrorException reinterpret(Float64, a)
end

@testset "convert" begin
    a = [RGB(1,0,0) RGB(0,0,1);
         RGB(0,1,0) RGB(1,1,1)]
    c = @inferred(convert(Array{BGR}, a))
    @test eltype(c) == BGR{U8}
    c = @inferred(convert(Array{BGR{Float32}}, a))
    @test eltype(c) == BGR{Float32}
    c = @inferred(convert(Array{Lab}, a))
    @test eltype(c) == Lab{Float32}
end

@testset "eltype conversion" begin
    @test float32(Float64) == Float32
    @test float32(U8)      == Float32
    @test float64(RGB{U8}) == RGB{Float64}

    a = [RGB(1,0,0) RGB(0,0,1);
         RGB(0,1,0) RGB(1,1,1)]
    @test eltype(a) == RGB{U8}
    @test eltype(u8.(a))       == RGB{U8}
    @test eltype(ufixed8.(a))  == RGB{U8}
    @test eltype(ufixed10.(a)) == RGB{UFixed10}
    @test eltype(ufixed12.(a)) == RGB{UFixed12}
    @test eltype(ufixed14.(a)) == RGB{UFixed14}
    @test eltype(ufixed16.(a)) == RGB{U16}
    @test eltype(u16.(a))      == RGB{U16}
#    @test eltype(float16.(a)) == RGB{Float16}
    @test eltype(float32.(a)) == RGB{Float32}
    @test eltype(float64.(a)) == RGB{Float64}

    a = U8[0.1,0.2,0.3]
    @test eltype(a) == U8
    @test eltype(u8.(a))       == U8
    @test eltype(ufixed8.(a))  == U8
    @test eltype(ufixed10.(a)) == UFixed10
    @test eltype(ufixed12.(a)) == UFixed12
    @test eltype(ufixed14.(a)) == UFixed14
    @test eltype(ufixed16.(a)) == U16
    @test eltype(u16.(a))      == U16
#    @test eltype(float16.(a)) == Float16
    @test eltype(float32.(a)) == Float32
    @test eltype(float64.(a)) == Float64

    a = OffsetArray(U8[0.1,0.2,0.3], -1:1)
    @test eltype(a) == U8
    @test eltype(u8.(a))       == U8
    @test eltype(ufixed8.(a))  == U8
    @test eltype(ufixed10.(a)) == UFixed10
    @test eltype(ufixed12.(a)) == UFixed12
    @test eltype(ufixed14.(a)) == UFixed14
    @test eltype(ufixed16.(a)) == U16
    @test eltype(u16.(a))      == U16
#    @test eltype(float16.(a)) == Float16
    @test eltype(float32.(a)) == Float32
    @test eltype(float64.(a)) == Float64
    @test indices(float32.(a)) == (-1:1,)
end

nothing


text:
"""
Google Cloud Platform service-account API credentials.
"""
module credentials

export Credentials, JSONCredentials, MetadataCredentials

import Base: show, print
import JSON
import MbedTLS
import HTTP
using URIParser

using ..error
using ..root

abstract type Credentials end

"""
Get credential/service-account information from GCE metadata server

See [Storing and Retreiving Instance Metadata](https://cloud.google.com/compute/docs/storing-retrieving-metadata)
"""
struct MetadataCredentials <: Credentials
    url::String
    service_account::String
    project_id::String
    client_email::String
    scopes::Vector{String}
    function MetadataCredentials(; url::AbstractString=METADATA_ROOT, service_account::AbstractString="default")
        metadata = new(url, service_account)
        project_id, client_email, scopes = try
            (
                get(metadata, "project-id"; context=:project),
                get(metadata, "email"; context=:service_account),
                split(chomp(get(metadata, "scopes"; context=:service_account)), '\n')
            )
        catch e
            throw(CredentialError("Unable to contact metadata server"))
        end
        new(url, service_account, project_id, client_email, scopes)
    end
end

function Base.get(credentials::MetadataCredentials, path::AbstractString; context::Symbol=:service_account)
    headers = Dict{String, String}("Metadata-Flavor" => "Google")
    if context == :service_account
        url = joinpath(credentials.url, "instance/service-accounts/$(credentials.service_account)")
    elseif context == :instance
        url = joinpath(credentials.url, "instance")
    elseif context == :project
        url = joinpath(credentials.url, "project")
    else
        throw(CredentialError("Unknown metadata context: $context"))
    end
    res = HTTP.get(joinpath(url, path), headers)
    if res.status != 200
        throw(CredentialError("Unable to obtain credentials from metadata server"))
    end
    String(res.body)
end

"""
    JSONCredentials(...)

Parse JSON credentials created for a service-account at
[Google Cloud Platform Console](https://console.cloud.google.com/apis/credentials)
"""
struct JSONCredentials <: Credentials
    account_type::String
    project_id::String
    private_key_id::String
    private_key::MbedTLS.PKContext
    client_email::String
    client_id::String
    auth_uri::URI
    token_uri::URI
    auth_provider_x509_cert_url::URI
    client_x509_cert_url::URI
    function JSONCredentials(account_type::AbstractString, project_id::AbstractString,
                             private_key_id::AbstractString, private_key::AbstractString,
                             client_email::AbstractString, client_id::AbstractString,
                             auth_uri::AbstractString, token_uri::AbstractString,
                             auth_provider_x509_cert_url::AbstractString,
                             client_x509_cert_url::AbstractString)
        ctx = MbedTLS.PKContext()
        MbedTLS.parse_key!(ctx, private_key)
        new(
            account_type, project_id, private_key_id, ctx, client_email, client_id,
            URI(auth_uri), URI(token_uri), URI(auth_provider_x509_cert_url), URI(client_x509_cert_url)
        )
    end
end

"""
    JSONCredentials(data::Dict{Symbol, String})

Initialise credentials from dictionary containing values.
"""
function JSONCredentials(data::AbstractDict{Symbol, <: AbstractString})
    fields = [fieldnames(JSONCredentials)...]
    fields[findfirst(fields, :account_type)] = :type  # type is a keyword!
    missing = setdiff(fields, keys(data))
    if !isempty(missing)
        info(missing)
        throw(CredentialError("Missing fields in key: ", join(missing, ", ")))
    end
    JSONCredentials((data[field] for field in fields)...)
end

"""
    JSONCredentials(filename)

Load credentials from a JSON file.
"""
JSONCredentials(filename::AbstractString) = JSONCredentials(JSON.parsefile(filename; dicttype=Dict{Symbol, String}))
JSONCredentials(io::IO) = JSONCredentials(JSON.parse(io; dicttype=Dict{Symbol, String}))

Base.convert(::Type{JSONCredentials}, x::AbstractString) = JSONCredentials(x)

function print(io::IO, x::Credentials)
    fields = [:project_id, :client_email]
    print(io, join(("$field: $(getfield(x, field))" for field in fields), "\n"))
end
show(io::IO, x::JSONCredentials) = print(io, x)

end

ext: jl


text:
export fgl_deriv

function fgl_deriv( a, y, h )
    #   Computes the fractional derivative of order alpha (a) for the function
    #   y sampled on a regular grid with spacing h, using the Grunwald-Letnikov
    #   formulation.
    #   Note that this implementation is similar to that of Bayat 2007
    #   (FileExchange: 13858-fractional-differentiator), and takes the exact
    #   same inputs, but uses a vectorized formulation to accelerates the
    #   computation in Matlab.
    #   Copyright 2014 Jonathan Hadida
    #   Contact: Jonathan dot hadida [a] dtc.ox.ac.uk
        #using LinearAlgebra SpecialFunctions # added for Julia 1.1
        n     = length(y);
        J     = collect(0:(n-1));
        #G1    = gamma.( J+1 );
        #G2    = gamma.( a+1-J );
        G1    = gamma.( J+convert(Array{Int64,1},ones(size(J))));
        test21 = convert(Array{Float64,1},J);
        test22 = (a+1)*ones(size(J));
        G2    = gamma.(-test21+test22);
        s     = (-1) .^ J;

        #M     = tril( ones(n,n) );
        M     = LowerTriangular( ones(n,n) );
        test  = y;
        R     = zeros(n,n);
        for i = 1:n
            R[i,i:n] = y[1:n - i + 1]
            R[i:n,i] = y[1:n - i + 1]
        end
        #T     = repmat( (gamma.(a+1)/(h^a)) * s ./ (G1.*G2),1,size(M,2));
        T     = repeat( (gamma.(a+1)/(h^a)) * s ./ (G1.*G2),1,size(M,2));
        #Y     = reshape(sum( R .* M .* T, 2 ), size(y));
        Y     = reshape(sum( R .* M .* T, dims=2 ), size(y));
        #Y[isnan.(Y)] = 0;
        inds = findall(isnan,Y);
        for i = 1:length(inds)
        Y[inds[i]] = 0
        end
        return Y
    end


text:
using PyCall

py"""
import json
import os
import pickle
import math

def send_packet(packet: dict) -> bool:
    res = False
    if 'ip' in packet['_source']['layers'].keys():
        n = int(packet['_source']['layers']['ip']['ip.ttl'])
        if (n & (n - 1) == 0) and n != 0:
            res = True
    return res

def prepare_sample(filename='tmp.json'):
    f = open(filename)
    pcap = json.load(f)
    f.close()
    packet_arr = []
    for i, packet in enumerate(pcap):
        # tmp = [direction, time diff to next recieved,
        #        time diff to next send, packet size]
        direction = int(send_packet(packet))
        time_diff_send = None
        time_diff_rcv = None
        time = float(packet['_source']['layers']['frame']['frame.time_relative'])
        for next_packet in pcap[i+1:]:
            tmp_direction = send_packet(next_packet)
            if not time_diff_rcv and not tmp_direction:
                time_diff_rcv = float(next_packet['_source']
                                       ['layers']['frame']
                                       ['frame.time_relative']) - time
            if not time_diff_send and tmp_direction:
                time_diff_send = float(next_packet['_source']
                                      ['layers']['frame']
                                      ['frame.time_relative']) - time

        packet_size = int(packet['_source']['layers']['frame']['frame.len'])
        tmp = [direction, math.log(1+time_diff_rcv),
               math.log(1+time_diff_send), float(packet_size)/2**16]
        tmp += extract_packet_features(packet)
        packet_arr.append(tmp)

        if i == 255:
            break
    return packet_arr

def create_samples(mal_path, clean_path, mal_sample_path, clean_sample_path):
    mal_files = [os.path.join(mal_path, f)
                 for f in os.listdir(mal_path)
                 if os.path.isfile(os.path.join(mal_path, f))]
    clean_files = [os.path.join(clean_path, f)
                   for f in os.listdir(clean_path)
                   if os.path.isfile(os.path.join(clean_path, f))]

    for i, mal in enumerate(mal_files):
        if i > 393:
            print(i, mal)
            try:
                os.system('tshark -r {} -T json > tmp.json'.format(mal))
                with open(mal_sample_path+str(i), 'wb') as out_:
                    pickle.dump(prepare_sample(), out_)
            except:
                print('error')
            os.system('rm tmp.json')

    for i, clean in enumerate(clean_files):
        print(i, clean)
        try:
            os.system('tshark -r {} -T json > tmp.json'.format(clean))
            with open(clean_sample_path+str(i), 'wb') as out_:
                pickle.dump(prepare_sample(), out_)
        except:
            print('error')
        os.system('rm tmp.json')
"""

@doc """
    convert_jsons_2_pickle_samples(mal_path::String, clean_path::Stirng,
        mal_sample_path::String, clean_sample_path::String)

Loads pcap samples stored on paths 'mal_path' and 'clean_path', extract features and
stored them in pickle format on 'mal_sample_path' and 'clean_sample_path'.

From each packet are extract these features: packet direction,
time difference to next recieved, time difference to next send, packet size.
""" ->
function convert_jsons_2_pickle_samples(mal_path::String, clean_path::Stirng,
            mal_sample_path::String, clean_sample_path::String)
    py"create_samples"(mal_path, clean_path, mal_sample_path, clean_sample_path)
end

max_stars_repo_name: kozvojtex/MMWavenet

max_issues_repo_name: kozvojtex/MMWavenet

max_forks_repo_name: kozvojtex/MMWavenet


text:
#-------------------------------------------------------------------------------
# bearing and range available

"""
    $TYPEDEF

Bearing and Range constraint from a Pose2 to Point2 variable.
"""
mutable struct Pose2Point2BearingRange{B <: IIF.SamplableBelief, R <: IIF.SamplableBelief} <: IIF.AbstractRelativeMinimize
    bearing::B
    range::R
end

function getSample(cfo::CalcFactor{<:Pose2Point2BearingRange}, N::Int=1)
  smpls = zeros(2, N)
  smpls[1,:] = rand(cfo.factor.bearing, N)[:]
  smpls[2,:] = rand(cfo.factor.range, N)[:]

  # must return at least first element in `::Tuple` as `::Matrix`
  return (smpls,)
end

function IIF.getParametricMeasurement(s::Pose2Point2BearingRange{<:Normal, <:Normal})

  meas = [mean(s.bearing), mean(s.range)]
  iΣ = [1/var(s.bearing)             0;
                      0  1/var(s.range)]

  return meas, iΣ
end

# TODO consolidate with parametric constraint, follow at #467
function (cfo::CalcFactor{<:Pose2Point2BearingRange})(meas, xi, lm)
  # 1-bearing
  # 2-range
  # world frame
  θ = meas[1] + xi[3]
  mx = meas[2]*cos(θ)
  my = meas[2]*sin(θ)

  ex = lm[1] - (mx + xi[1])
  ey = lm[2] - (my + xi[2])

  # res = [eθ, er]
  eθ = atan((my + xi[2]), (mx + xi[1])) - atan(lm[2], lm[1])  # eθ
  er= sqrt(ex^2 + ey^2) # some wasted computation here       # er

  #Old way
  # rot = meas[1]+xi[3]
  # res[1] = ( lm[1] - (meas[2]*cos( rot ) + xi[1]) )^2
  # res[2] = ( lm[2] - (meas[2]*sin( rot ) + xi[2]) )^2
  # res[1] += res[2]
  # res[2] = 0.0
  # return res[1]

  # IIF v0.21+
  return [eθ, er]
end

# quick check
# pose = (0,0,0),  bear = 0.0,  range = 10.0   ==>  lm = (10,0)
# pose = (0,0,0),  bear = pi/2,  range = 10.0   ==>  lm = (0,10)
# pose = (0,0,pi/2),  bear = 0.0,  range = 10.0   ==>  lm = (0,10)
# pose = (0,0,pi/2),  bear = pi/2,  range = 10.0   ==>  lm = (-10,0)

# Support for database based solving

passTypeThrough(d::FunctionNodeData{Pose2Point2Range}) = d

mutable struct PackedPose2Point2BearingRange <: IncrementalInference.PackedInferenceType
    bearstr::String
    rangstr::String
end

function convert( ::Type{PackedPose2Point2BearingRange}, d::Pose2Point2BearingRange )
  return PackedPose2Point2BearingRange( convert(PackedSamplableBelief, d.bearing), convert(PackedSamplableBelief, d.range) )
end

# TODO -- should not be resorting to string, consider specialized code for parametric distribution types and KDEs
function convert( ::Type{Pose2Point2BearingRange}, d::PackedPose2Point2BearingRange )
  Pose2Point2BearingRange( convert(SamplableBelief, d.bearstr), convert(SamplableBelief, d.rangstr) )
end


text:
# # PTDF with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Sourabh Dalvi

# ## Introduction

# PowerSimulations.jl supports linear PTDF optimal power flow formulation. This example shows a
# single multi-period optimization of economic dispatch with a linearilzed DC-OPF representation of
# using PTDF power flow and how to extract duals values or locational marginal prices for energy.

# ## Dependencies
# We can use the same RTS data and some of the initialization as in
# [OperationsProblem example](../../notebook/3_PowerSimulations_examples/1_operations_problems.ipynb)
# by sourcing it as a dependency.
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(
    pkgpath,
    "test",
    "3_PowerSimulations_examples",
    "01_operations_problems.jl",
));

# Since we'll be retreving duals, we need a solver that returns duals values
# here we use Ipopt.
using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)

# In the [OperationsProblem example](../../notebook/3_PowerSimulations_examples/1_operations_problems.ipynb)
# we defined a unit-commitment problem with a copper plate representation of the network.
# Here, we want do define an economic dispatch (linear generation decisions) with
# linear DC-OPF using PTDF network representation.
# So, starting with the network, we can select from _almost_ any of the endpoints on this
# tree:
#nb TypeTree(PSI.PM.AbstractPowerModel,  init_expand = 10, scopesep="\n")

# For now, let's just choose a standard PTDF formulation.
ed_template = template_economic_dispatch(network = StandardPTDFModel)

# Currently  energy budget data isn't stored in the RTS-GMLC dataset.
ed_template.devices[:Hydro] = DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver)

# Calculate the PTDF matrix.
PTDF_matrix = PTDF(sys)

# Now we can build a 4-hour economic dispatch / PTDF problem with the RTS data.
# Here, we have to pass the keyword argument `constraint_duals` to OperationsProblem
# with the name of the constraint for which duals are required for them to be returned in the results.
problem = OperationsProblem(
    EconomicDispatchProblem,
    ed_template,
    sys,
    horizon = 4,
    optimizer = solver,
    balance_slack_variables = true,
    constraint_duals = [:CopperPlateBalance, :network_flow],
    PTDF = PTDF_matrix,
)

# And solve the problem and collect the results
res = solve!(problem)

# Here we collect the dual values from the results for the `CopperPlateBalance` and `network_flow`
# constraints. In the case of PTDF network formulation we need to compute the final LMP for each bus in the system by
# subtracting the duals (μ) of `network_flow` constraint multipled by the PTDF matrix
# from the  dual (λ) of `CopperPlateBalance` constraint.
# Note:we convert the results from DataFrame to Array for ease of use.
λ = convert(Array, res.dual_values[:CopperPlateBalance])
μ = convert(Array, res.dual_values[:network_flow])

# Here we create Dict to store the calculate congestion component of the LMP which is a product of μ and the PTDF matrix.
buses = get_components(Bus, sys)
congestion_lmp = Dict()
for bus in buses
    congestion_lmp[get_name(bus)] = μ * PTDF_matrix[:, get_number(bus)]
end
congestion_lmp = DataFrame(congestion_lmp)

# Finally here we get the LMP for each node in a lossless DC-OPF using the PTDF formulation.
LMP = λ .- congestion_lmp


text:
using Test
using Covers

cc=Cover(4,8)
cc(1,1,true)
cc(2,2,true)

@test all(cc .^ Bool[1 1 1 1 1 1 1 1;
                     0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0])

@test all(ncovered(cc) .== [1 2 1 1 1 1 1 1; 0 1 0 0 0 0 0 0; 0 1 0 0 0 0 0 0; 0 1 0 0 0 0 0 0])

@test all(multicovered(cc) .== Bool[0 1 0 0 0 0 0 0; 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0])

@test cc[9] && (!cc[10])

cc[9] =false
@test all(cc .^ Bool[0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0])

cc[10] =true

@test all(cc .^ Bool[0 1 1 0 0 0 0 0;
                     0 1 1 0 0 0 0 0;
                     0 1 1 0 0 0 0 0;
                     0 1 1 0 0 0 0 0])
cc[9] =true

@test all(cc .^ Bool[0 1 1 0 0 0 0 0;
                     0 1 1 0 0 0 0 0;
                     0 1 1 0 0 0 0 0;
                     0 1 1 0 0 0 0 0])

cc[10] =false
@test all(cc .^ Bool[0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0;
                     0 1 0 0 0 0 0 0])

cc=Cover(4,8)
cc(1,1,true)
cc(2,2,true)

@test cc[1,3] && (!cc[2,3])

cc[1,3] =false
@test all(cc .^ Bool[0 1 0 0 0 0 0 0;
                  0 1 0 0 0 0 0 0;
                  0 1 0 0 0 0 0 0;
                  0 1 0 0 0 0 0 0])

cc[2,3] =true

@test all(cc .^ Bool[0 1 1 0 0 0 0 0;
                  0 1 1 0 0 0 0 0;
                  0 1 1 0 0 0 0 0;
                  0 1 1 0 0 0 0 0])

cc.=false

@test all(cc .^ Bool[0 0 0 0 0 0 0 0;
                  0 0 0 0 0 0 0 0;
                  0 0 0 0 0 0 0 0;
                  0 0 0 0 0 0 0 0])


text:
# multiple dispatch:
f(n, m) = "base case"
f(n::Number, m::Number) = "n and m are both numbers"
f(n::Number, m) = "n is a number"
f(n, m::Number) = "m is a number"
f(n::Integer, m::Integer) = "n and m are both integers"
#> f (generic function with 5 methods)
f(1.5, 2) #> "n and m are both numbers"
f(1, "bar") #> "n is a number"
f(1, 2) #> "n and m are both integers"
f("foo", [1,2]) #> "base case"
f(n::Float64, m::Integer) = "n is a float and m is an integer"
#> f (generic function with 6 methods)
f(1.5, 2) #> "n is a float and m is an integer"

# methods:
methods(+)
methods(sort)
methodswith(String)

# measuring execution:
fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)
@time fib(35)
@elapsed fib(35) #> elapsed time: 0.115853481 seconds (45144 bytes allocated)

max_stars_repo_path: Module 1/Chapter03/generic_functions.jl

max_issues_repo_path: Module 1/Chapter03/generic_functions.jl

max_forks_repo_path: Module 1/Chapter03/generic_functions.jl


text:
export knot, knotc, knotu

"""
    knot(npts, ord, centre=false, step=1.0)
Generate a `ord` order B-spline open uniform knot vector over `npts = n+1` control vertices.

An open uniform knot vector has multiplicity of knot values
 at the ends equals to the order `ord` of the B-Spline basis function.
Internal knot values are evenly spaced (uniform).
The resulting open uniform basis functions yield curves
 that behave most nearly like Bézier curves.
`centre` Boolean value indicates if the knots are centered in the
 zero value or if them goes from zero to `m-1` where `m=npts+ord`.
 If `true` is choosen, then `m` must be odd.
Necessary condition is that `npts>=ord` in order to correct build the knot vector.

---

# Arguments
- `npts::Int64`: the number of the `n+1` points of the controll polygon.
- `ord::Int64`: the order of the B-Spline (`degree k-1`).
- `centre::Bool=false`: if the knot is centered in zero.
- `step::Float64=1.0`: the step between the knots.

---

# Examples
```jldoctest
julia> knot(5, 2)
7-element Array{Float64,1}:
 0.0
 0.0
 1.0
 2.0
 3.0
 4.0
 4.0
```

```jldoctest
julia> knot(6, 3, false, 0.25)
9-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.25
 0.5
 0.75
 1.0
 1.0
 1.0
```

```jldoctest
julia> knot(5, 2, true)
7-element Array{Float64,1}:
 -2.0
 -2.0
 -1.0
 0.0
 1.0
 2.0
 2.0
```

---

_By Elia Onofri, Giuseppe Santorelli_
"""

function knot(npts::Int64, ord::Int64, centre::Bool=false, step::Float64=1.0)::Array{Float64}
    # Input Data Check
    @assert npts >= ord ("ERROR: ord > npts")

    local m::Int64 = npts+ord      		   # knot vector dimension
    local backstep::Float64 = 0.0   		   # backstep is the affine transition parameter for 0-central knots
    local x::Array{Float64} = Array{Float64}(m)    # knot vector output

    # Different preparation for centre knot vector
    if centre
        @assert m % 2 == 1 ("ERROR: no centring points available")
        backstep = step * (npts - ord + 1) / 2
    end

    x[1] = 0 - backstep
    for i = 2 : m
        if i > ord && i <= npts + 1
            x[i] = x[i-1] + step
        else
            x[i] = x[i-1]
        end
    end
    return x
end

#--------------------------------------------------------------------------------------------------------------------------------

"""
knotc(npts, ord, b)

Generate a nonuniform open knot vector proportional to the chord lengths between defining polygon vertices.

# Arguments

- `npts::Int64`: the number of the contol polygon vertices.
- `ord::Int64`: order of the basis function.
- `b::Array{Float64}`: 2-dimensional array containing the contol polygon vertices, every column represents a point.

---

# Examples

```jldoctest
julia> b = [0. 2 4 6 8; 0 6 3 6 6; 0 0 0 0 0]
3×5 Array{Float64,2}:
 0.0  2.0  4.0  6.0  8.0
 0.0  6.0  3.0  6.0  6.0
 0.0  0.0  0.0  0.0  0.0

julia> knotc(5, 3, b)
8-element Array{Float64,1}:
 0.0
 0.0
 0.0
 1.45338
 2.38171
 3.0
 3.0
 3.0
```

---
_Giuseppe Santorelli_
"""

function knotc(npts::Int64, ord::Int64, b::Array{Float64,2})::Array{Float64}

    @assert length(b) == 3 * npts ("ERROR: array b not matching with parameters")
    @assert npts >= ord ("ERROR: ord > npts")

    m::Int64 = npts + ord
    x::Array{Float64} = Array{Float64}(m)
    chord::Array{Float64} = Array{Float64}(npts - 1)

    maxchord = 0
    icount = 0

    # determine chord distance between defining polygon vertices and their sum
    for i = 2 : npts
        icount = icount + 1
        xchord = b[1, i] - b[1, i - 1]
        ychord = b[2, i] - b[2, i - 1]
        zchord = b[3, i] - b[3, i - 1]
        chord[icount] = sqrt(xchord * xchord + ychord * ychord + zchord * zchord)
        maxchord = maxchord + chord[icount]
    end

    # multiplicity of c=order zeros at the beginning of the open knot vector
    for i = 1 : ord
        x[i] = 0
    end

    # generate the internal knot values
    csum = 0
    for i = 1 : npts - ord
        csum = csum + chord[i]
        numerator = i / (npts - ord + 1) * chord[i + 1] + csum
        x[ord + i] = (numerator / maxchord) * (npts - ord + 1)
    end

    # multiplicity of c=order zeros at the end of the open knot vector
    for i = npts + 1 : m
        x[i] = npts - ord + 1
    end

    return x

end

#--------------------------------------------------------------------------------------------------------------------------------

"""
    knotu(npts, ord, step=1.0)
Generate a B-spline periodic uniform knot vector over `npts = n+1` control vertices.

In a periodic uniform knot vector the knot values are evenly spaced by `step`.
Moreover each base is a translate of the others as it could be defined like:

```math
    N_{i, k}(t) = N_{i-1, k}(t-1) = N_{i+1, k}(t+1)
```

---

# Arguments
- `npts::Int64`: the number of the `n+1` points of the controll polygon.
- `ord::Int64`: the order of the B-Spline (`deg = ord-1`).
- `step::Float64=1.0`: the step between the knots.

---

# Example
```jldoctest
julia> knotu(5, 2)
7-element Array{Float64,1}:
 0.0
 1.0
 2.0
 3.0
 4.0
 5.0
 6.0
```

```
julia> knotu(5, 4, 0.25)
9-element Array{Float64,1}:
 0.0
 0.25
 0.5
 0.75
 1.0
 1.25
 1.5
 1.75
 2.0
```

---
_By Elia Onofri, Giuseppe Santorelli_
"""

function knotu(npts::Int64, ord::Int64, step::Float64=1.0)::Array{Float64}
    # Input Data Check
    @assert npts >= ord ("ERROR: ord > npts")

    local m::Int64 = npts + ord      		   # knot vector dimension
    local x::Array{Float64} = Array{Float64}(m)    # knot vector output

    for i = 1 : m
        x[i] = (i-1) * step
    end

    return x
end


text:
function polyc_Req(Req::Vector, L0::Real, Kₓ::Real, Rtot::Vector, Cplx::AbstractMatrix, θ::Vector, Kav::Matrix)
    @assert size(Cplx) == (length(θ), size(Kav, 1))
    @assert size(Kav, 2) == length(Rtot)
    θ /= sum(θ)

    ψ = Kav .* Req' .* Kₓ
    ψRS = sum(ψ, dims = 2) .+ 1.0
    L0Kt = L0 / Kₓ * θ
    Lbounds = L0Kt .* expm1.(Cplx * log.(ψRS))
    Rbounds = L0Kt .* Cplx * (ψ ./ ψRS) .* exp.(Cplx * log.(ψRS))
    Lfbnds = L0Kt .* exp.(Cplx * log.(ψRS .- 1.0))
    Lmbnds = L0Kt .* Cplx * (ψRS .- 1.0)

    @assert length(Lbounds) == size(Cplx, 1)
    @assert size(Rbounds) == (length(θ), length(Rtot))
    return Lbounds, Rbounds, Lfbnds, Lmbnds
end

function polyc(args...)
    ansType = promote_type(typeof(args[1]), typeof(args[2]), eltype(args[3]), eltype(args[4]), eltype(args[5]))
    f! = (F, x) -> F .= args[3] .- vec(sum(polyc_Req(x, args...)[2], dims = 1)) .- x
    Req = rootSolve(f!, convert(Vector{ansType}, args[3]))
    return polyc_Req(Req, args...)
end

polycm = (Kₓ, Rtot, Cplx, Ltheta, Kav) -> polyc(sum(Ltheta), Kₓ, Rtot, Cplx, Ltheta ./ sum(Ltheta), Kav)


text:
using ChoroplethMaps
using DataFrames
using Base.Test

###############################################################################
# User-facing functionality:

pv = Provider.STATESUMMARY()
@test Provider.names(pv) == [:STUSPS, :AWATER, :LSAD, :AFFGEOID, :STATENS,
                             :GEOID, :ALAND, :STATEFP, :NAME]
@test Provider.projection(pv) == "epsg:4269"
@test size(Provider.polygons(pv), 1) == 52
@test Provider.polygons(pv)[1][1][1] == -118.593969

df = DataFrame(NAME=["Alabama", "California"],
       feature=[26.6673430422, 6.549492945] )

mp = mapify(df, pv, key=:NAME)
@test names(mp) == [:CM_X, :CM_Y, :NAME, :CM_ORIG_KEY, :feature]

# Does data look basically correct?
@test mp[1, :NAME] == "Alabama_1"
@test mp[1, :CM_X] == -9.848286458886575e6
@test mp[1, :CM_Y] == 3.7493855901445053e6
@test mp[1, :feature] == 26.6673430422

# High-level test of polygon splitting:
@test unique(mp[:NAME]) == ["Alabama_1", "California_1", "California_2",
                            "California_3", "California_4", "California_5",
                            "California_6"]

###############################################################################
# Non-user-facing functionality:

indices = ChoroplethMaps.findindices(DataFrame(CM_POLYGONS=Provider.polygons(pv)))
@test size(indices, 1) == 52
@test indices[1] == Array[[1,10], [11,19], [20,26], [27,43], [44,56], [57,467]]
@test indices[2] == Array[[1,10]]


text:
# script to run inference for the OU SDEMEM model

using Pkg
using LinearAlgebra
using DataFrames

println("start run script for kalman")

# load functions
include(pwd()*"/src/SDEMEM OU neuron data/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU neuron data/mcmc.jl")

M_subjects = 1

#seed = parse(Int,ARGS[1])

seed = 100

y,prior_parameters_η,prior_parameters_σ_ϵ = set_up(M_subjects,seed)

job = string(M_subjects)

# run MH-Gibbs
R = 100000 #500000
burn_in = 20000 #250000

startval_ϕ = ones(M_subjects,3)

for j = 1:3

    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

    startval_ϕ[:,j,1] = μ_0_j*ones(M_subjects)

end

startval_η = zeros(6)

for j = 1:3

    # get paramaters
    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

    startval_η[j] = μ_0_j # set start value to mean
    startval_η[j+3] = (α_j - 1)/β_j

end

# we should transform this parameter, to have the "walk" in a better space
startval_log_σ_ϵ = log(0.2)

# set time step
dt = diff(y[1].mSec)[1]

# estimate parameters using exact Gibbs sampling

ll = @time kalman(y, exp(startval_log_σ_ϵ), startval_ϕ, dt)
sum(ll)

Σ_i_σ_ϵ = 1^2

Σ_i_ϕ = Matrix{Float64}[]
for i in 1:M_subjects; push!(Σ_i_ϕ, [0.04 0 0;0 0.4 0;0 0 0.04]); end

γ_ϕ_0 = 1.
γ_σ_ϵ_0 = 1.
μ_i_ϕ = startval_ϕ
μ_i_σ_ϵ = startval_log_σ_ϵ
α_star_ϕ = 0.25
α_star_σ_ϵ = 0.25
log_λ_i_ϕ = log.(2.4/sqrt(3)*ones(M_subjects))
log_λ_i_σ_ϵ = log(2.4)
update_interval = 1
α_power = 0.7
start_update = 100

# set random numbers
Random.seed!(seed)

run_time_kalman = @elapsed chain_ϕ_kalman, chain_σ_ϵ_kalman, chain_η_kalman, accept_vec_kalman, loglik_vec = gibbs_exact(R,
                                                                                y,
                                                                                dt,
                                                                                Σ_i_σ_ϵ,
                                                                                Σ_i_ϕ,
                                                                                γ_ϕ_0,
                                                                                γ_σ_ϵ_0,
                                                                                μ_i_ϕ,
                                                                                μ_i_σ_ϵ,
                                                                                α_star_ϕ,
                                                                                α_star_σ_ϵ,
                                                                                log_λ_i_ϕ,
                                                                                log_λ_i_σ_ϵ,
                                                                                update_interval,
                                                                                start_update,
                                                                                α_power,
                                                                                startval_ϕ,
                                                                                startval_log_σ_ϵ,
                                                                                startval_η,
                                                                                prior_parameters_η,
                                                                                prior_parameters_σ_ϵ);

println(run_time_kalman)
println(sum(accept_vec_kalman[1,:])/(M_subjects*R))
println(sum(accept_vec_kalman[2,:])/R)
println(sum(accept_vec_kalman[3,:])/(3*R))

sim_data = zeros(4,1)

sim_data[1] = run_time_kalman
sim_data[2] = sum(accept_vec_kalman[1,:])/(M_subjects*R)
sim_data[3] = sum(accept_vec_kalman[2,:])/R
sim_data[4] = sum(accept_vec_kalman[3,:])/(3*R)

chain_ϕ_export = zeros(M_subjects*3, R)

idx = 0
for m = 1:M_subjects
    for j = 1:3
        global idx = idx + 1
        chain_ϕ_export[idx,:] = chain_ϕ_kalman[m,j,:]
    end
end

loglik_data = zeros(R,1)
loglik_data[:] = loglik_vec
DataFrame(loglik_data)

CSV.write("data/SDEMEM OU neuron data/kalman/sim_data_"*string(seed)*"_"*job*".csv", DataFrame(sim_data))
CSV.write("data/SDEMEM OU neuron data/kalman/chain_sigma_epsilon_"*string(seed)*"_"*job*".csv", DataFrame(chain_σ_ϵ_kalman'))
CSV.write("data/SDEMEM OU neuron data/kalman/chain_eta_"*string(seed)*"_"*job*".csv", DataFrame(chain_η_kalman'))
CSV.write("data/SDEMEM OU neuron data/kalman/chain_phi_"*string(seed)*"_"*job*".csv", DataFrame(chain_ϕ_export'))
CSV.write("data/SDEMEM OU neuron data/kalman/loglik_data_"*string(seed)*"_"*job*".csv", DataFrame(loglik_data))

println(Threads.nthreads())

println("end run script")

max_stars_repo_path: src/SDEMEM OU neuron data/run_script_kalman_one_subject.jl


text:
module DiffImages
using ImageCore,
      Zygote,
      ImageTransformations,
      StaticArrays,
      CoordinateTransformations,
      Interpolations,
      ChainRulesCore,
      LinearAlgebra,
      Rotations

using Zygote: @adjoint
using ChainRulesCore: NoTangent

export colorify, channelify
include("colors/conversions.jl")
include("geometry/warp.jl")
include("geometry/adjoints.jl")

end


text:
const respOk = Dict{String, Any}("status" => "ok")

export
    initDfg,
    toggleMockServer,
    registerRobot,
    registerSession

function initDfg(config, dfg, request)::Dict{String, Any}
end

"""
Configure to just return "{"status": "OK"}" for all requests.
"""
function toggleMockServer(config, dfg, request)::Dict{String, Any}
    @show request
    @show config
    if haskey(request, "payload")
        if haskey(request["payload"], "isMockServer")
            config["isMockServer"] = request["payload"]["isMockServer"]
            return Dict{String, Any}("isMockServer" => config["isMockServer"])
        else
            error("Please provide a 'payload' body containing a field 'isMockServer': [\"true\"/\"false\"] for this request.")
        end
    else
        error("Please provide a 'payload' body containing a field 'isMockServer': [\"true\"/\"false\"] for this request.")
    end
end

function registerRobot(config, dfg, request)::Dict{String, Any}
  @info "This is Caesar.registerRobot"
  config["robotId"] = request["payload"]["robot"]
  return respOk
end

function registerSession(config, dfg, request)::Dict{String, Any}
  @info "This is registerSession"
  if haskey(request, "robot")
    if config["robotId"] == request["payload"]["robot"]
      config["sessionId"] = request["payload"]["session"]
      return respOk
    else
      error("Currently only supports one robot, config has $(config["robotId"]), but you asked for robotId=$(request["robot"])")
    end
  end
end


text:
Pkg.add("Cairo")
Pkg.add("Luxor")
Pkg.add("ProgressMeter")
Pkg.add("ForwardDiff")
Pkg.add("FFTW")
Pkg.add("Zygote")
Pkg.add("DataStructures")
Pkg.add("StatsBase")
Pkg.add(Pkg.PackageSpec(name="CuArrays", version="1.7.3"))
Pkg.add("CuArrays")
Pkg.add("Flux")
Pkg.add(["CUDAdrv", "CUDAnative", "CuArrays"])
Pkg.test(["CUDAdrv", "CUDAnative", "CuArrays"])
Pkg.add("HTTP")
Pkg.add("JSON")
Pkg.add("ProfileView")
Pkg.add("PackageCompiler")
Pkg.add("GeometricalPredicates")

# https://github.com/JuliaGeometry/Clipper.jl
Pkg.add("Clipper")

Pkg.rm("CuArrays")

Pkg.update()

import Profile
function test_profile()
    Profile.clear()
    Profile.init(n = 10^7, delay = 0.01)

    open("/tmp/prof.txt", "w") do s
        # I want to ignore all low-cost operations (e.g. <10)
        Profile.print(IOContext(s, :displaysize => (24, 200)), mincount=10)
    end

    # other visualization methods
    #
    # Profile.print()
    # import ProfileView
    # ProfileView.view()
end

import Flux
import CuArrays
Flux.gpu(x) = x
Flux.gpu(x) = fmap(CuArrays.cu, x)

max_stars_repo_name: lihebi/bhdl

max_issues_repo_name: lihebi/bhdl

max_forks_repo_name: lihebi/bhdl


text:
@testset "435.non-overlapping-intervals.jl" begin
    @test erase_overlap_intervals([(1, 2), (2, 3), (3, 4), (1, 3)]) == 1
    @test erase_overlap_intervals([(1, 2), (1, 2), (1, 2)]) == 2
    @test erase_overlap_intervals([(1, 2), (2, 3)]) == 0
    @test erase_overlap_intervals(Tuple{Int, Int}[]) == 0
end


text:
using LinearAlgebra, DiffEqOperators, Random, Test

# Generate random parameters
al = rand(5)
bl = rand(5)
cl = rand(5)
dx = rand(5)
ar = rand(5)
br = rand(5)
cr = rand(5)

# Construct 5 arbitrary RobinBC operators
for i in 1:5
    Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]), dx[i])

    Q_L, Q_b = Array(Q,5i)

    #Check that Q_L is is correctly computed
    @test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
    @test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
    @test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]

    #Check that Q_b is computed correctly
    @test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]

    # Construct the extended operator and check that it correctly extends u to a (5i+2)
    # vector, along with encoding boundary condition information.
    u = rand(5i)

    Qextended = Q*u
    CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
    @test length(Qextended) ≈ 5i+2

    # Check concretization
    @test Array(Qextended) ≈ CorrectQextended

    # Check that Q_L and Q_b correctly compute BoundaryPaddedVector
    @test Q_L*u + Q_b ≈ CorrectQextended

    @test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended

end

# Construct 5 arbitrary RobinBC operators w/non-uniform grid
al = rand(5)
bl = rand(5)
cl = rand(5)
dx = rand(5)
ar = rand(5)
br = rand(5)
cr = rand(5)
for j in 1:2
    for i in 1:5
        if j == 1
            Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]),
                        dx[i] .* ones(5 * i))
        else
            Q = RobinBC([al[i], bl[i], cl[i]], [ar[i], br[i], cr[i]],
                        dx[i] .* ones(5 * i))
        end
        Q_L, Q_b = Array(Q,5i)

        #Check that Q_L is is correctly computed
        @test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
        @test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
        @test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]

        #Check that Q_b is computed correctly
        @test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]

        # Construct the extended operator and check that it correctly extends u to a (5i+2)
        # vector, along with encoding boundary condition information.
        u = rand(5i)

        Qextended = Q*u
        CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
        @test length(Qextended) ≈ 5i+2

        # Check concretization
        @test Array(Qextended) ≈ CorrectQextended

        # Check that Q_L and Q_b correctly compute BoundaryPaddedVector
        @test Q_L*u + Q_b ≈ CorrectQextended

        @test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended

    end
end

#3rd order RobinBC, calculated with left stencil [-11/6 3 -3/2 1/3], right stencil [-1/3 3/2 -3 11/6] and [α,β,γ] = [1 6 10]
u0 = -4/10
uend = 125/12
u = Vector(1.0:10.0)
Q = RobinBC((1.0, 6.0, 10.0), (1.0, 6.0, 10.0), 1.0, 3)
urobinextended = Q*u
@test urobinextended.l ≈ u0
@test urobinextended.r ≈ uend
# General BC should be equivalent
G = GeneralBC([-10.0, 1.0, 6.0], [-10.0, 1.0, 6.0], 1.0, 3)
ugeneralextended = G*u
@test ugeneralextended.l ≈ u0
@test ugeneralextended.r ≈ uend

#TODO: Implement tests for BC's that are contingent on the sign of the coefficient on the operator near the boundary


text:
#=
   Euler #2 in Julia.

  Problem 2
  """
  Each new term in the Fibonacci sequence is generated by adding the
  previous two terms. By starting with 1 and 2, the first 10 terms will be:

  1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

  Find the sum of all the even-valued terms in the sequence which do not
  exceed four million.
  """

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

using Memoization
include("Euler.jl")

# 7.5558e-5
function euler2a()
    n = 100 # some arbitrary large value
    f = zeros(BigInt,(n,1))
    f[1:2] = [0,1]
    [f[i] = f[i-1]+f[i-2] for i in 3:n]
    sum([f[i]*(f[i]%2 == 0 && f[i] < 4000000) for i in 1:n])
end

const fibmem2b = Dict{Int,Int}()
function fib2b(n::Int)
    get!(fibmem2b, n) do
        n < 3 ? 1 : fib2b(n-1) + fib2b(n-2)
    end
end

# 8.108e-6s
function euler2b()
    n = 1
    f = 0
    s = 0
    while f < 4_000_000
        f = fib2b(n)
        if f % 2 == 0
            s += f
        end
        n += 1
    end
    return s
end

#
# This is slower than using a Dict!
#
@memoize function fib2c(n::Int)
    n < 3 ? 1 : fib2c(n-1) + fib2c(n-2)
end

# 0.008s
function euler2c()
    n = 1
    s = 0
    while (f = fib2c(n)) < 4_000_000
        if f % 2 == 0
            s += f
        end
        n += 1
    end
    return s
end

# Using no cache
function fib_no_cache(n::Int)
  if n <= 2
     return 1
  end
  return fib_no_cache(n-1)+fib_no_cache(n-2)
end

# 0.0379s
function euler2d()
    n = 1
    s = 0
    while (f = fib_no_cache(n)) < 4_000_000
        if f % 2 == 0
            s += f
        end
        n += 1
    end
    return s
end

# run_euler(euler2a)
run_euler(euler2b)
# run_euler(euler2c)
# run_euler(euler2d)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
const RAIL_COLORS = (
    # preserve yellow for future performance linting
    45, # light blue
    123, # light cyan
    150, # ???
    215, # orange
    231, # white
)
const N_RAILS = length(RAIL_COLORS)
const TYPE_ANNOTATION_COLOR = :light_cyan

printlnstyled(args...; kwarg...) = printstyled(args..., '\n'; kwarg...)

function print_rails(io, depth)
    for i = 1:depth
        color = RAIL_COLORS[i%N_RAILS+1]
        printstyled(io, '│'; color)
    end
end

max_stars_repo_name: aviatesk/AbstractInterpreterPlayground

max_issues_repo_name: aviatesk/AbstractInterpreterPlayground

max_forks_repo_name: aviatesk/AbstractInterpreterPlayground


text:
# Optimization utilities:
# - FISTA: Beck, A., and Teboulle, M., 2009, A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems
# - Anderson acceleration: Wei et al., 2021, Stochastic Anderson Mixing for Nonconvex Stochastic Optimization

export FISTA, FISTA!, Anderson
export reset!, linesearch_backtracking, spectral_radius

## FISTA

abstract type AbstractFISTA{T}<:Flux.Optimise.AbstractOptimiser end

mutable struct FISTA{T}<:AbstractFISTA{T}
    L::T
    prox::Function
    t::T
    Nesterov::Bool
end

mutable struct FISTA!{T}<:AbstractFISTA{T}
    L::T
    prox!::Function
    t::T
    Nesterov::Bool
end

FISTA(L::T, prox::Function; t::T=T(1), Nesterov::Bool=true) where {T<:Real} = FISTA{T}(L, prox, t, Nesterov)

FISTA!(L::T, prox!::Function; t::T=T(1), Nesterov::Bool=true) where {T<:Real} = FISTA!{T}(L, prox!, t, Nesterov)

reset!(opt::AbstractFISTA{T}) where {T<:Real} = (opt.t = T(1))

function Flux.Optimise.apply!(opt::AbstractFISTA{T}, x::AbstractArray{CT,N}, g::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}

    # Gradient + proxy update
    steplength = T(1)/opt.L
    if opt isa FISTA
        g .= x-opt.prox(x-steplength*g, steplength)
    elseif opt isa FISTA!
        opt.prox!(x-steplength*g, steplength, g); g .= x.-g
    end

    # Nesterov momentum
    if opt.Nesterov
        t = T(0.5)*(T(1)+sqrt(T(1)+T(4)*opt.t^2))
        g .*= (t+opt.t-T(1))/t
        opt.t = t
    end

    return g

end

function spectral_radius(A::AT, x::AbstractArray{T,N}; niter::Int64=10) where {T,N,AT<:Union{AbstractMatrix{T},AbstractLinearOperator{T,N,N}}}
    x = x/norm(x)
    y = similar(x)
    ρ = real(T)(0)
    for _ = 1:niter
        y .= A*x
        ρ = norm(y)/norm(x)
        x .= y/norm(y)
    end
    return ρ
end

## Anderson acceleration

mutable struct Anderson{T}<:Flux.Optimise.AbstractOptimiser
    lr::Union{Nothing,T}
    hist_size::Integer
    β::T
    ΔX::Union{Nothing,Vector{<:AbstractArray}}
    ΔG::Union{Nothing,Vector{<:AbstractArray}}
    x0::Union{Nothing,AbstractArray}
    g0::Union{Nothing,AbstractArray}
    init::Bool
end

Anderson(; lr::Union{Nothing,T}=nothing, hist_size::Integer=5, β::T=1.0) where {T<:Real} = Anderson{T}(lr, hist_size, β, nothing, nothing, nothing, nothing, false)

reset!(opt::Anderson) = (opt.ΔX = nothing; opt.ΔG = nothing; opt.x0 = nothing; opt.g0 = nothing; opt.init = false)

function Flux.Optimise.apply!(opt::Anderson{T}, x::AbstractArray{CT,N}, g::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}

    # Rescaling gradient w/ learning rate
    ~isnothing(opt.lr) && (g .*= opt.lr)

    # Initialize if necessary
    if ~opt.init
        opt.ΔX = Vector{typeof(x)}(undef,0)
        opt.ΔG = Vector{typeof(g)}(undef,0)
        opt.x0 = deepcopy(x)
        opt.g0 = deepcopy(g)
        opt.init = true
        return g
    end

    # Update model/residual history
    if length(opt.ΔX) < opt.hist_size
        push!(opt.ΔX, x-opt.x0)
        push!(opt.ΔG, g-opt.g0)
    else
        for i = 1:opt.hist_size-1
            opt.ΔX[i] .= opt.ΔX[i+1]
            opt.ΔG[i] .= opt.ΔG[i+1]
        end
        opt.ΔX[end] .= x-opt.x0
        opt.ΔG[end] .= g-opt.g0
    end
    opt.x0 .= x
    opt.g0 .= g

    # Computing steplength
    Γ = solve_Anderson_steplength(opt.ΔG, g)

    # Gradient update
    return g .= sum(opt.ΔX.*Γ)+opt.β*(g-sum(opt.ΔG.*Γ))

end

function solve_Anderson_steplength(ΔR::AbstractArray{AT,1}, r::AbstractArray{CT}) where {T<:Real,CT<:RealOrComplex{T},AT<:AbstractArray{CT}}
    n = length(ΔR)
    nx = prod(size(ΔR[1]))
    A = Array{CT,2}(undef,nx,n)
    @inbounds for j = 1:n
        A[:,j] = vec(ΔR[j])
    end
    return A\vec(r)
end

function linesearch_backtracking(obj::Function, x0::AbstractArray{CT,N}, p::AbstractArray{CT,N}, lr::T; fx0::Union{Nothing,T}=nothing, niter::Integer=3, mult_fact::T=T(0.5), verbose::Bool=false) where {T<:Real,N,CT<:RealOrComplex{T}}

    isnothing(fx0) && (fx0 = obj(x0))
    for n = 1:niter
        fx = obj(x0+lr*p)
        fx < fx0 ? break : (verbose && print("."); lr *= mult_fact)
    end
    return lr

end

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
export MtlComputeCommandEncoder
export set_function!, set_buffer!, dispatchThreads!, endEncoding!
export set_buffers!, append_current_function!

const MTLComputeCommandEncoder = Ptr{MtComputeCommandEncoder}

mutable struct MtlComputeCommandEncoder <: MtlCommandEncoder
    handle::MTLComputeCommandEncoder
    cmdbuf::MtlCommandBuffer
end

Base.unsafe_convert(::Type{MTLComputeCommandEncoder}, q::MtlComputeCommandEncoder) = q.handle

function MtlComputeCommandEncoder(cmdbuf::MtlCommandBuffer; dispatch_type::Union{Nothing,MtDispatchType} = nothing)
    if isnothing(dispatch_type)
        handle = mtNewComputeCommandEncoder(cmdbuf)
    else
        handle = mtNewComputeCommandEncoderWithDispatchtype(cmdbuf, dispatchtype)
    end
    obj = MtlComputeCommandEncoder(handle, cmdbuf)
    #finalizer(unsafe_destroy!, obj)
    return obj
end

device(cce::MtlComputeCommandEncoder) = cce.cmdbuf.device

set_function!(cce::MtlComputeCommandEncoder, pip::MtlComputePipelineState) =
    mtComputeCommandEncoderSetComputePipelineState(cce, pip)

set_buffer!(cce::MtlComputeCommandEncoder, buf::MtlBuffer, offset::Integer, index::Integer) =
    mtComputeCommandEncoderSetBufferOffsetAtIndex(cce, buf, offset, index - 1)
#set_bufferoffset!(cce::MtlComputeCommandEncoder, offset::Integer, index::Integer) =
#    mtComputeCommandEncoderBufferSetOffsetAtIndex(cce, offset, index)
set_buffers!(cce::MtlComputeCommandEncoder, bufs::Vector{T},
             offsets::Vector{Int}, indices::UnitRange{Int}) where {T<:MtlBuffer} =
    mtComputeCommandEncoderSetBuffersOffsetsWithRange(cce, handle_array(bufs), offsets, indices .- 1)
#=set_buffers!(cce::MtlComputeCommandEncoder, bufs::Vector{MtlPtr{T}},
             offsets::Vector{Int}, indices::UnitRange{Int}) where {T} =
    mtComputeCommandEncoderSetBuffersOffsetsWithRange(cce, bufs, offsets, indices .- 1)=#

set_bytes!(cce::MtlComputeCommandEncoder, ptr, length::Integer, index::Integer) =
    mtComputeCommandEncoderSetBytesLengthAtIndex(cce, ptr, length, index - 1)

dispatchThreads!(cce::MtlComputeCommandEncoder, gridSize::MtSize, threadGroupSize::MtSize) =
    mtComputeCommandEncoderDispatchThreadgroups_threadsPerThreadgroup(cce, gridSize, threadGroupSize)

# higher level
set_argument!(cce::MtlComputeCommandEncoder, index::Integer, buf::MtlBuffer, offset::Integer = 0) =
    set_buffer!(cce, buf, offset, index)

function set_argument!(cce::MtlComputeCommandEncoder, index::Integer, val::T) where T
    @assert T.isbitstype

    ref = Base.RefValue(val)
    ptr = Base.unsafe_convert(Ptr{T}, ref)
    set_bytes!(cce, ptr, sizeof(T), index)
end

function set_arguments!(cce::MtlComputeCommandEncoder, args...)
    for (i,arg) = enumerate(args)
        set_argument!(cce, i, arg)
    end
end
#####
# encode in the Command Encoder
function MtlComputeCommandEncoder(f::Base.Callable, cmdbuf::MtlCommandBuffer; kwargs...)
    encoder = MtlComputeCommandEncoder(cmdbuf; kwargs...)
    f(encoder)
    close(encoder)
    return encoder
end

append_current_function!(cce::MtlComputeCommandEncoder, gridSize::MtSize, threadGroupSize::MtSize) =
    dispatchThreads!(cce, gridSize, threadGroupSize)

#### use
use!(cce::MtlComputeCommandEncoder, buf::MtlBuffer, mode::MtResourceUsage=ReadWriteUsage) =
    mtComputeCommandEncoderUseResourceUsage(cce, buf, mode)

use!(cce::MtlComputeCommandEncoder, buf::Vector{MtlBuffer}, mode::MtResourceUsage=ReadWriteUsage) =
    mtComputeCommandEncoderUseResourceCountUsage(cce, handle_array(buf), length(buf), mode)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the transformations between Geodetic and Geocentric frames.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/geodetic_geocentric.jl
# ==================================================

# Function: ecef_to_geodetic
# --------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-3: Converting ECEF to Lat Lon [1, p. 173].
#
# According to this example, using:
#
#   r_ecef = 6524.834 i + 6862.875 j + 6448.296 k [km]
#
# one gets:
#
#   Geodetic Latitude  = 34.352496°
#            Longitude = 46.4464°
#            Altitude  = 5085.22 km
#
# Scenario 02
# ===========
#
# At the poles, we have a singularity. We know that, if:
#
#   r_ecef = 0 i + 0 j + Z k
#
# then
#
#   Geodetic Latitude  = 90° for Z > 0 and -90° for Z < 0
#            Longitude =  0°
#            Altitude  = Z - b_wgs84
#
################################################################################

@testset "Function ecef_to_geodetic" begin
    # Scenario 01
    # ===========

    r = [6524.834e3, 6862.875e3, 6448.296e3]

    ϕ_gd, λ_gd, h = ecef_to_geodetic(r)

    @test rad2deg(ϕ_gd) ≈ 34.352496 atol = 1e-6
    @test rad2deg(λ_gd) ≈ 46.4464   atol = 1e-4
    @test h/1000        ≈ 5085.22   atol = 1e-2

    # Scenario 02
    # ===========

    aux = rand(0:1000)

    Z = R0 + aux
    ϕ_gd, λ_gd, h = ecef_to_geodetic([0;0;Z])

    @test rad2deg(ϕ_gd) ≈ 90
    @test rad2deg(λ_gd) ≈ 0
    @test h             ≈ Z - SatelliteToolbox.b_wgs84

    Z = -R0 + aux
    ϕ_gd, λ_gd, h = ecef_to_geodetic([0;0;Z])

    @test rad2deg(ϕ_gd) ≈ -90
    @test rad2deg(λ_gd) ≈ 0
    @test h             ≈ -Z - SatelliteToolbox.b_wgs84
end

################################################################################
#                                 Test Results
################################################################################
#
# Using the FORTRAN code available in:
#
#   https://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm
#
# we obtained the following matrix of data:
#
# | Geocentric lat [rad] | Distance [m] | Geodetic lat [rad]   | Altitude [m]        |
# |----------------------|--------------|----------------------|---------------------|
# |  19.86               | R0 + 1987    |  1.0134554245512695  | 17352.756962650223  |
# |  π/2                 | R0 + 190686  |  1.5707963267948966  | 212070.68575482070  |
# |  31.35 * π/180       | R0 + 752000  |  0.54808101129276687 | 757773.37237201189  |
# | -19.86               | R0 + 1987    | -1.0134554245512695  | 17352.756962650223  |
# | -π/2                 | R0 + 190686  | -1.5707963267948966  | 212070.68575482070  |
# | -31.35 * π/180       | R0 + 752000  | -0.54808101129276687 | 757773.37237201189  |
# |  0                   | R0           |  0.0000000000000000  | 0.0000000000000000  |
# |  15                  | 10000        |  1.3567978765139961  | -6353137.8454352869 |
# | -15                  | 10000        | -1.3567978765139961  | -6353137.8454352869 |
#
################################################################################
#
# FORTRAN code used to obtain the results
# ==============================================================================
#
#      subroutine GEOD(r,z,fi,h)
#          implicit real*8(a-h,o-z)
#          data a,frec /6378137.d0,298.257223563d0/
#          b=dsign(a-a/frec,z)
#          if(r.eq.0d0) return
#          E=((z+b)*b/a-a)/r
#          F=((z-b)*b/a+a)/r
#          P=(E*F+1.)*4d0/3.
#          Q=(E*E-F*F)*2.
#          D=P*P*P+Q*Q
#          if(D.ge.0d0) then
#              s=dsqrt(D)+Q
#              s=dsign(dexp(dlog(dabs(s))/3d0),s)
#              v=P/s-s
#              v=-(Q+Q+v*v*v)/(3*P)
#          else
#              v=2.*dsqrt(-P)*dcos(dacos(Q/P/dsqrt(-P))/3.)
#          endif
#          G=.5*(E+dsqrt(E*E+v))
#          t=dsqrt(G*G+(F-v*G)/(G+G-E))-G
#          fi=datan((1.-t*t)*a/(2*b*t))
#          h=(r-a*t)*dcos(fi)+(z-b)*dsin(fi)
#          end
#
#      program TEST
#          implicit real*8(a-h,o-z)
#          R0 = 6378137.d0
#          PI = 4 * atan(1.d0)
#
#          dlat = 19.86d0
#          r    = R0 + 1987.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = 90.00d0 * PI / 180.d0
#          r    = R0 + 190686.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = 31.25d0 * PI / 180.d0
#          r    = R0 + 752000.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = -19.86d0
#          r    = R0 + 1987.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = -90.00d0 * PI / 180.d0
#          r    = R0 + 190686.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = -31.25d0 * PI / 180.d0
#          r    = R0 + 752000.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = 0.d0
#          r    = R0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = 15.d0 * PI / 180.d0
#          r    = 10000.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#
#          dlat = -15.d0 * PI / 180.d0
#          r    = 10000.d0
#          call GEOD(r * dcos(dlat), r * dsin(dlat), fi, h)
#          PRINT *, fi, h
#      end program
#
################################################################################

# Function: geocentric_to_geodetic
# --------------------------------

@testset "Function geocentric_to_geodetic" begin
    # North hemisphere
    # ==========================================================================

    ϕ_gd, h = geocentric_to_geodetic(19.86, R0 + 1987)
    @test ϕ_gd ≈ 1.0134554245512695
    @test h    ≈ 17352.756962650223

    ϕ_gd, h = geocentric_to_geodetic(π/2, R0 + 190686)
    @test ϕ_gd ≈ 1.5707963267948966
    @test h    ≈ 212070.68575482070

    ϕ_gd, h = geocentric_to_geodetic(deg2rad(31.25), R0 + 752000)
    @test ϕ_gd ≈ 0.54808101129276687
    @test h    ≈ 757773.37237201189

    # South hemisphere
    # ==========================================================================

    ϕ_gd, h = geocentric_to_geodetic(-19.86, R0 + 1987)
    @test ϕ_gd ≈ -1.0134554245512695
    @test h    ≈ 17352.756962650223

    ϕ_gd, h = geocentric_to_geodetic(-π/2, R0 + 190686)
    @test ϕ_gd ≈ -1.5707963267948966
    @test h    ≈ 212070.68575482070

    ϕ_gd, h = geocentric_to_geodetic(-deg2rad(31.25), R0 + 752000)
    @test ϕ_gd ≈ -0.54808101129276687
    @test h    ≈ 757773.37237201189

    # Special cases
    # ==========================================================================

    ϕ_gd, h = geocentric_to_geodetic(0, R0)
    @test ϕ_gd ≈ 0.0
    @test h    ≈ 0.0

    # D < 0
    ϕ_gd, h = geocentric_to_geodetic(deg2rad(15), 10e3)
    @test ϕ_gd ≈ 1.3567978765139961
    @test h    ≈ -6353137.8454352869

    ϕ_gd, h = geocentric_to_geodetic(-deg2rad(15), 10e3)
    @test ϕ_gd ≈ -1.3567978765139961
    @test h    ≈ -6353137.8454352869
end


text:
using DancingLinks.Sudoku

using Test

DUMMY = zeros(Int, 9, 9)

# From https://norvig.com/sudoku.html
A_TYPICAL_PUZZLE = [
    4 0 0 0 0 0 8 0 5
    0 3 0 0 0 0 0 0 0
    0 0 0 7 0 0 0 0 0
    0 2 0 0 0 0 0 6 0
    0 0 0 0 8 0 4 0 0
    0 0 0 0 1 0 0 0 0
    0 0 0 6 0 3 0 7 0
    5 0 0 2 0 0 0 0 0
    1 0 4 0 0 0 0 0 0
]

# Arto Inkala's "the most difficult sudoku-puzzle known so far"
HARDEST_PUZZLE_2006 = [
    8 5 0 0 0 2 4 0 0
    7 2 0 0 0 0 0 0 9
    0 0 4 0 0 0 0 0 0
    0 0 0 1 0 7 0 0 2
    3 0 5 0 0 0 9 0 0
    0 4 0 0 0 0 0 0 0
    0 0 0 0 8 0 0 7 0
    0 1 7 0 0 0 0 0 0
    0 0 0 0 3 6 0 4 0
]

# Arto Inkala's "the most difficult puzzle I've ever created"
HARDEST_PUZZLE_2010 = [
    0 0 5 3 0 0 0 0 0
    8 0 0 0 0 0 0 2 0
    0 7 0 0 1 0 5 0 0
    4 0 0 0 0 5 3 0 0
    0 1 0 0 7 0 0 0 6
    0 0 3 2 0 0 0 8 0
    0 6 0 5 0 0 0 0 9
    0 0 4 0 0 0 0 3 0
    0 0 0 0 0 9 7 0 0
]

NO_SOLUTION = ones(Int, 9, 9)

dummy = create_sudoku(DUMMY)
a_typical_puzzle = create_sudoku(A_TYPICAL_PUZZLE)
hardest_puzzle_2006 = create_sudoku(HARDEST_PUZZLE_2006)
hardest_puzzle_2010 = create_sudoku(HARDEST_PUZZLE_2010)
no_solution = create_sudoku(NO_SOLUTION)

@test check_solution(solve_sudoku!(dummy), DUMMY)
@test check_solution(solve_sudoku!(a_typical_puzzle), A_TYPICAL_PUZZLE)
@test check_solution(solve_sudoku!(hardest_puzzle_2006), HARDEST_PUZZLE_2006)
@test check_solution(solve_sudoku!(hardest_puzzle_2010), HARDEST_PUZZLE_2010)

@test solve_sudoku!(no_solution) == nothing

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
"""
Package: Forecast

    seaborne(query = "imports")

Return estimates of world seaborne trade from AIS data collected by Marine Traffic.

The data contains deadweight imports and exprts for France, Germany and the United Kingdom from 2015-04-01 to 2021-05-02. The fields are:

num_pc:     number of port calls
mtc:        metric tons of cargo
dwt:        deadweight tonnage
suffix_ma:  30-day moving averages

Data available at UN COMTRADE Monitor.Cerdeiro, Komaromi, Liu and Saeed (2020).

# Arguments
- `query`: Defaults to "paper" and accepts the following values:
    - "full": Full original dataset from from 2015-04-01 to 2021-05-02
    - "imports": Subset with imports.

# Returns
Dataframe containing information based on the parameter `query`

# Examples
```julia-repl
julia> seaborne()
[ Info: Seaborne deadweight trade imports from AIS
2199×4 DataFrame
  Row │ Date        France   Germany  UK
      │ Date        Int64    Int64    Int64
──────┼──────────────────────────────────────
    1 │ 2015-04-01   507946   878377  599573
    2 │ 2015-04-02   332043  1501614  772714
    3 │ 2015-04-03   810077   941663  262994
   [...]
```
"""
function seaborne(query::String = "imports")

    @assert query in ["full","imports"]

    path = art_path("seaborne")

    sb_df = open(joinpath(path, "seaborne.csv")) do file
        CSV.read(file,DataFrame,
                 dateformat = "mm/dd/yyyy HH:MM:SS pp",
                 types = Dict(:date => Date))
    end

    if query == "full"
        @info "Full dataset with estimates of world seaborne trade from AIS"
        return sb_df
    end

    if query == "imports"
        @info "Seaborne deadweight trade imports from AIS"

        # Group by Country and Flow to select Imports
        gbt = groupby(sb_df, [:country_name, :flow], sort=true)
        return DataFrame(Dict(:Date => gbt[1].date ,
                              :UK => gbt[1].dwt, :Germany => gbt[3].dwt, :France => gbt[5].dwt))

    end
end


text:
using DashBootstrapComponents

alert = html_div([
    dbc_button("Toggle", id = "alert-toggle-auto", className = "me-1", n_clicks = 0),
    html_hr(),
    dbc_alert(
        "Hello! I am an auto-dismissing alert!",
        id = "alert-auto",
        is_open = true,
        duration = 4000,
    ),
]);

callback!(
    app,
    Output("alert-auto", "is_open"),
    Input("alert-toggle-auto", "n_clicks"),
    State("alert-auto", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
# ------------------------------------------------------------------------------------------
# # Metaprogramming with Julia
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Scare Quotes
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia supports *metaprogramming*. This is similar to symbolic programming, where we deal
# with expressions (like $2+2$) as opposed to values (like $4$).
#
# Normally, Julia takes all code we give it as a set of instructions, and carries them out.
# If we type `2+2` it will faithfully add those numbers and give us the result.
# ------------------------------------------------------------------------------------------

2+2

# ------------------------------------------------------------------------------------------
# We can prevent this from happening with quotation marks. Surrounding our code with `"`
# treats it as a literal string of characters, without seeing it as code at all.
# ------------------------------------------------------------------------------------------

x = "2+2"

# ------------------------------------------------------------------------------------------
# We can then explicitly tell Julia to evaluate it later on.
# ------------------------------------------------------------------------------------------

eval(Meta.parse(ans))

# ------------------------------------------------------------------------------------------
# Why go through this complicated way to add numbers? The trick is that when we have the
# expression `2+2`, we can modify it in various interesting ways. As a simple example,
# imagine replacing `+` with `-`.
# ------------------------------------------------------------------------------------------

x = replace(x, "+"=>"-")

eval(Meta.parse(x))

# ------------------------------------------------------------------------------------------
# We don't actually want to work with strings here; Julia has a much more powerful way to
# quote code, the frowny face operator `:()`.
# ------------------------------------------------------------------------------------------

x = :(2+2)

eval(x)

# ------------------------------------------------------------------------------------------
# We can quote larger expression, including blocks and entire function definitions. The
# keyword `quote` is an alternative to `begin` that returns the quoted block.
#
# In larger blocks, Julia will preserve line number information, which appears as a comment.
# ------------------------------------------------------------------------------------------

quote
  x = 2 + 2
  hypot(x, 5)
end

:(function mysum(xs)
    sum = 0
    for x in xs
      sum += x
    end
  end)

# ------------------------------------------------------------------------------------------
# ## Fruit of the Expression Tree
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Strings support "interpolation", which allows us to easily build larger strings from
# smaller components.
# ------------------------------------------------------------------------------------------

x = "yields falsehood when preceded by its quotation"
println(x)

y = "'$x' $x"
println(y)

x = :(2+2)

:($x * $x)

eval(ans)

# ------------------------------------------------------------------------------------------
# ## The Root of all Eval
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `eval` can do more than just returning a result. What happens if we quote something like a
# function definition?
# ------------------------------------------------------------------------------------------

ex = :(foo() = println("I'm foo!"))

# ------------------------------------------------------------------------------------------
# It doesn't actually do anything; yet.
# ------------------------------------------------------------------------------------------

foo()  # throws UndefVarError

# ------------------------------------------------------------------------------------------
# But evaluating `ex` brings `foo` to life!
# ------------------------------------------------------------------------------------------

eval(ex)

foo()

# ------------------------------------------------------------------------------------------
# Using interpolation, we can construct a function definition on-the-fly; in fact, we can
# make a whole series of functions at once.
# ------------------------------------------------------------------------------------------

for name in [:foo, :bar, :baz]
  println(:($name() = println($("I'm $(name)!"))))
end

# ------------------------------------------------------------------------------------------
# And then bring them to life with `eval`, too.
# ------------------------------------------------------------------------------------------

for name in [:foo, :bar, :baz]
  eval(:($name() = println($("I'm $(name)!"))))
end

bar()

baz()

# ------------------------------------------------------------------------------------------
# This can be an *extremely* useful trick when wrapping APIs (say, from a C library or over
# HTTP). APIs often define a list of available functions, so you can grab that and generate
# the whole wrapper at once! See Clang.jl, TensorFlow.jl, or the Base linear algebra
# wrappers for examples.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Original sin
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Here's a more practical example. Consider the following definition of the `sin` function,
# based on the Taylor series.
#
# $$sin(x) = \sum_{k=1}^{\infty} \frac{(-1)^k}{(1+2k)!} x^{1+2k}$$
# ------------------------------------------------------------------------------------------

mysin(x) = sum((-1)^k/factorial(1+2k) * x^(1+2k) for k = 0:5)

mysin(0.5), sin(0.5)

# ------------------------------------------------------------------------------------------
# To see where we are right now, we'll benchmark it.
# ------------------------------------------------------------------------------------------

using BenchmarkTools
@benchmark mysin(0.5)

# ------------------------------------------------------------------------------------------
# Right now, this is much slower than it could be. The reason is that we're looping over
# `k`, which is relatively expensive. It'd be much faster to write out:
# ------------------------------------------------------------------------------------------

mysin(x) = x - x^3/6 + x^5/120 # + ...

# ------------------------------------------------------------------------------------------
# But this is tedious to write, and no longer looks like the original Taylor series. It's
# harder to tell if we've made a mistake, and we easily modify it. Is there a way to get the
# best of both worlds?
#
# How about getting Julia to write out that code for us?
#
# To start with, let's consider a symbolic version of the `+` function.
# ------------------------------------------------------------------------------------------

plus(a, b) = :($a + $b)

plus(1, 2)

# ------------------------------------------------------------------------------------------
# With `plus` we can do more interesting things, like symbolic `sum`:
# ------------------------------------------------------------------------------------------

reduce(+, 1:10)

reduce(plus, 1:10)

eval(ans)

# ------------------------------------------------------------------------------------------
# Given that, we can also sum over symbolic variables.
# ------------------------------------------------------------------------------------------

reduce(plus, [:(x^2), :x, 1])

# ------------------------------------------------------------------------------------------
# This gives us an important piece of the puzzle, but we also need to figure out _what_
# we're summing. Let's crate a symbolic version of the Taylor series above, which
# interpolates the value of `k`.
# ------------------------------------------------------------------------------------------

k = 2
:($((-1)^k) * x^$(1+2k) / $(factorial(1+2k)))

# ------------------------------------------------------------------------------------------
# Now we have one term, we can generate as many as we like.
# ------------------------------------------------------------------------------------------

terms = [:($((-1)^k) * x^$(1+2k) / $(factorial(1+2k))) for k = 0:5]

# ------------------------------------------------------------------------------------------
# And sum them –
# ------------------------------------------------------------------------------------------

reduce(plus, ans)

# ------------------------------------------------------------------------------------------
# And create a function definition out of it:
# ------------------------------------------------------------------------------------------

:(mysin(x) = $ans)

eval(ans)

mysin(0.5), sin(0.5)

@benchmark mysin(0.5)

# ------------------------------------------------------------------------------------------
# On my machine `sin2` takes about 50 *nano*seconds to run – not bad for a naive
# implementation. If we challenged a photon to a twenty metre sprint, we'd win!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Beneath the Expression
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# It's really just a normal tree data structure, and we can take a look inside it.
# ------------------------------------------------------------------------------------------

ex = :(f(x, y))
ex = Expr(:call, :f, :x, :y)
@show ex.head;
@show ex.args;

# ------------------------------------------------------------------------------------------
# In our example above, we replace `+` with `-` in a string. In an expression, we can do
# that by altering the `.args` of the expression.
# ------------------------------------------------------------------------------------------

ex = :(2+2)

ex.args[1] = :-
ex

eval(ex)

# ------------------------------------------------------------------------------------------
# Note that larger, more complex expressions are a bit trickier than this. They are
# *nested*, which means that the expression `2+3` is contained inside the larger expression
# `1 + (2 + 3)`.
# ------------------------------------------------------------------------------------------

ex = :(1 + (2 + 3))

ex.args

# ------------------------------------------------------------------------------------------
# A package called MacroTools provides a way to deal with this. It shows us *all* sub-
# expressions in turn, allowing us to decide how we want to change them. It can be thought
# of as a bit like a find-and-replace operation. Here's an example that finds all integers
# in an expression, and increments them.
# ------------------------------------------------------------------------------------------

# import Pkg; Pkg.add("MacroTools")
using MacroTools
using MacroTools: postwalk

postwalk(ex) do x
  x isa Integer ? x+1 : x
end

# ------------------------------------------------------------------------------------------
# To get a feel for what's happening, you can use `@show` to see what `postwalk` sees.
#
# (`@show` is Julia's single most useful feature; if you don't understand what code is
# doing, throw it in and see what's happening.)
# ------------------------------------------------------------------------------------------

map(x -> @show(x), [1,2,3])

postwalk(ex) do x
  @show x
end

# ------------------------------------------------------------------------------------------
# MacroTools also provides tools for *pattern matching* over expressions. `a_ + b_` acts as
# a template; if the expression provided look like the template, `a` and `b` will match the
# two things be added. If not, they'll just be `nothing`.
# ------------------------------------------------------------------------------------------

ex = :(2 + 3)
@capture(ex, a_ + b_)

a, b

ex = :(f(2,3))
@capture(ex, a_ + b_)

a, b

# ------------------------------------------------------------------------------------------
# We can finally use this to replace _all_ `+`s with `-`s in an expression, rather than just
# one.
# ------------------------------------------------------------------------------------------

ex = :(3x^2 + (2x + 1))

postwalk(ex) do x
  @capture(x, a_ + b_) || return x
  :($a - $b)
end

# ------------------------------------------------------------------------------------------
# ## Macro Agressions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# You have probably seen macros already – basic features in Julia like `@show`, `@time` and
# `@inline` are really macros. For basic usage, it's enough to think of them as simple
# annotations that alter how code is run. But we now know enough to dig into how they work
# under the hood.
#
# Normal functions never see _code_, only _values_. If we pass `2+2` to a function, it sees
# `4`.
# ------------------------------------------------------------------------------------------

function foo(x)
  @show x
  return x
end

foo(2+2)

# ------------------------------------------------------------------------------------------
# Macros are much like function, but they see _code_ that was passed to them, like the
# expressions that we saw above. Macros get a chance to manipulate this code and alter how
# it behaves.
# ------------------------------------------------------------------------------------------

macro foo(x)
  @show x
  return x
end

@foo(2+2)

# ------------------------------------------------------------------------------------------
# As a simple example, we can replace one of the arguments to `+` and get a different result
# back.
# ------------------------------------------------------------------------------------------

macro foo(x)
  x.args[2] = 5
  return x
end

@foo(2+2)

# ------------------------------------------------------------------------------------------
# MacroTools provides a useful tool, `@expand`, to see what's going on inside a macro; it
# reveals the code that the macro returns without running it.
# ------------------------------------------------------------------------------------------

@expand @foo(2+2)

# ------------------------------------------------------------------------------------------
# You can of course use this on the macros that come with Julia, and it's a good way to
# learn how they work.
# ------------------------------------------------------------------------------------------

@time 2+2

@expand @time 2+2

# ------------------------------------------------------------------------------------------
# Let's do something more advanced. Using the techniques we developed above, we can find-
# and-replace `+` expressions in code. We'll still do the addition, but we'll also log what
# we added for debugging purposes.
#
# First, let's make an example expression using quotation.
# ------------------------------------------------------------------------------------------

ex = quote
  s = 0
  for x in xs
    s = s + x
  end
  return s
end

# ------------------------------------------------------------------------------------------
# We can use `ex` to develop the code transformation we want, and check that the code does
# the right thing.}
# ------------------------------------------------------------------------------------------

postwalk(ex) do x
  @capture(x, a_ + b_) || return x
  quote
    println("Adding " * string($a) * " to " * string($b))
    $x
  end
end

# ------------------------------------------------------------------------------------------
# Now we just wrap this in a macro, and we can add it to a normal function definition!
# ------------------------------------------------------------------------------------------

macro log_adds(ex)
  postwalk(ex) do x
    @capture(x, a_ + b_) || return x
    quote
      println("Adding " * string($a) * " to " * string($b))
      $x
    end
  end
end

@log_adds function mysum(xs)
  sum = 0
  for x in xs
    sum = sum + x
  end
  return sum
end

mysum(1:10)

# ------------------------------------------------------------------------------------------
# Adding and removing `@log_adds` is much nicer than inserting the debug calls by hand,
# especially if you have lots of `+`s in the code. Perhaps you can extend the macro to
# support logging other operators, like `-`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Generated Functions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# _Generated functions_ are a new metaprogramming tool unique to Julia. This section will
# briefly describe and motivate them, and [this blog
# post](http://mikeinnes.github.io/2017/08/24/cudanative.html) gives more examples for the
# interested reader.
#
# Essentially, a generated function is like a macro that operates on _types_ rather than
# expression trees. If we try to inspect arguments, we'll see their types rather than their
# values.
# ------------------------------------------------------------------------------------------

@generated function gadd(a, b)
  Core.println("a = $a, b = $b")
  :(a+b)
end

gadd(5, 2.5)

# ------------------------------------------------------------------------------------------
# Generated functions become more powerful when working with more complex types. For
# example, notice how arrays store their rank inside the type, so this is something we can
# generate code for.
# ------------------------------------------------------------------------------------------

rand(2,2)

typeof(ans)

# ------------------------------------------------------------------------------------------
# Why might this be useful? If you write code to deal with arrays, you'll notice that you
# often need a `for` loop for each dimension of the array. If you want to work with a 7D
# array, you need to write seven nested loops!
# ------------------------------------------------------------------------------------------

function mysum(xs::Array{<:Any,1})
  sum = 0
  for i = 1:length(xs)
    sum += xs[i]
  end
  return sum
end

function mysum(xs::Array{<:Any,2})
  sum = 0
  for i = 1:size(xs,1)
    for j = 1:size(xs, 2)
      sum += xs[i]
    end
  end
  return sum
end

# ------------------------------------------------------------------------------------------
# Where other languages simply hard-code a version of each function for vectors, matrices
# and perhaps 3D arrays, Julia allows us to easily write N-dimensional algorithms by
# generating nested loops as needed.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
alphafoldfile(uniprotXname) = "AF-$uniprotXname-F1-model_v1.pdb"

"""
    try_download_alphafold(uniprotXname, path=alphafoldfile(uniprotXname))

Attempt to download an [AlphaFold](https://alphafold.com/) structure.
Returns `nothing` if no entry corresponding to `uniprotXname` exists; otherwise
it return `path`, the pathname of the saved file.
"""
function try_download_alphafold(uniprotXname::AbstractString, path::AbstractString=alphafoldfile(uniprotXname))
    fn = alphafoldfile(uniprotXname)
    try
        Downloads.download("https://alphafold.ebi.ac.uk/files/$fn", path)
    catch err
        if isa(err, RequestError) && err.response.status == 404
            return nothing
        else
            throw(err)
        end
    end
    return path
end

"""
    download_alphafolds(msa::AbstractMultipleSequenceAlignment; dirname=pwd())

Download all available [AlphaFold](https://alphafold.com/) structures for the sequences in `msa`.
Missing entries are silently skipped.
"""
function download_alphafolds(msa::AbstractMultipleSequenceAlignment; dirname=pwd())
    for name in sequencenames(msa)
        uname = uniprotX(name)
        path = joinpath(dirname, alphafoldfile(uname))
        if !isfile(path)
            try_download_alphafold(uname, path)
        end
    end
end

"""
    getchain(filename::AbstractString; model="1", chain="A")

Read a PDB file `filename` and extract the specified chain.
"""
getchain(filename::AbstractString; model="1", chain="A") = read(filename, PDBFile; model, chain)

function countitems(list)
    count = Dict{eltype(list),Int}()
    for item in list
        count[item] = get(count, item, 0) + 1
    end
    return count
end

function findall_subseq(subseq, seq)
    starts = Int[]
    n = length(subseq)
    for i in eachindex(seq)
        if seq[i] == subseq[begin] && i + n - 1 <= lastindex(seq)
            issame = true
            for (ss, s) in zip(view(seq, i:i+n-1), subseq)
                issame &= ss == s
            end
            if issame
                push!(starts, i)
            end
        end
    end
    return starts
end