text:
#= You are a thief who has a bag with some capacity, you broke into a hous
and you found some valuables, you want to take the maximum value back
with you. You are confused what to pick, no worries, coders are here
to help you in your burglary. =#

## Function

function knapsack(weights, values, items, capacity)
    if (items == 0 || capacity == 0)
        return 0
    end
    if (weights[items] > capacity)
        return knapsack(weights, values, items - 1, capacity)
    end
    return max(
        knapsack(weights, values, items - 1, capacity - weights[items]) + values[items],
        knapsack(weights, values, items - 1, capacity),
    )
end

## Input

println("Enter the number of items")
n = readline()
n = parse(Int64, n)
println("Enter the capacity of the knapsack")
w = readline()
w = parse(Int64, w)
values = Int64[]
println("Enter the values of each item")
for i = 1:n
    temp = readline()
    temp = parse(Int64, temp)
    push!(values, temp)
end
weights = Int64[]
println("Enter the weights of each item")
for i = 1:n
    temp = readline()
    temp = parse(Int64, temp)
    push!(weights, temp)
end

## Calling the function

print("The maximum value you could take is ")
print(knapsack(weights, values, n, w))

#=
Sample Test Case:
Input:
    Enter the number of items
    3
    Enter the capacity of the knapsack
    50
    Enter the values of each item
    60 100 120
    Enter the weights of each item
    10 20 30
Ouput:
    The maximum value you could take is 220
=#

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
push!(LOAD_PATH, pwd()*"/../../src/")
using Distributions
using PyPlot
using BVB

x = linspace(-4,4,500)
max_iter = 90
lower = -10.
upper = 10.

init_approx_d = Normal(0.0, 10)

# example 1
@show true_d = TDist(1.0)
ll(x) = logpdf(true_d, x)

include("run_laplacian_gradient_boosting_exclusive.jl")
show()

println("Done.")
readline(STDIN)

ext: jl

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
"""
    stack(t, by = pkeynames(t); select = Not(by), variable = :variable, value = :value)`

Reshape a table from the wide to the long format. Columns in `by` are kept as indexing columns.
Columns in `select` are stacked. In addition to the id columns, two additional columns labeled
`variable` and `value` are added, containing the column identifier and the stacked columns.
See also [`unstack`](@ref).

# Examples

    t = table(1:4, names = [:x], pkey=:x)
    t = pushcol(t, :xsquare, :x => x -> x^2)
    t = pushcol(t, :xcube  , :x => x -> x^3)

    stack(t)
"""
function stack(t::D, by = pkeynames(t); select = isa(t, NDSparse) ? valuenames(t) : excludecols(t, by),
    variable = :variable, value = :value) where {D<:Dataset}

    (by != pkeynames(t)) && return stack(reindex(t, by, select); variable = :variable, value = :value)

    valuecols = columns(t, select)
    valuecol = [valuecol[i] for i in 1:length(t) for valuecol in valuecols]

    labels = fieldnames(typeof(valuecols))
    labelcol = [label for i in 1:length(t) for label in labels]

    bycols = map(arg -> repeat(arg, inner = length(valuecols)), columns(t, by))
    convert(collectiontype(D), Columns(bycols), Columns((labelcol, valuecol), names = [variable, value]))
end

"""
    unstack(t, by = pkeynames(t); variable = :variable, value = :value)

Reshape a table from the long to the wide format. Columns in `by` are kept as indexing columns.
Keyword arguments `variable` and `value` denote which column contains the column identifier and
which the corresponding values.  See also [`stack`](@ref).

# Examples

    t = table(1:4, [1, 4, 9, 16], [1, 8, 27, 64], names = [:x, :xsquare, :xcube], pkey = :x);

    long = stack(t)

    unstack(long)
"""
function unstack(t::D, by = pkeynames(t); variable = :variable, value = :value) where {D<:Dataset}
    tgrp = groupby((value => identity,), t, by, select = (variable, value))
    S = eltype(colnames(t))
    cols = S.(union(columns(t, variable)))
    T = eltype(columns(t, value))
    unstack(D, Base.nonmissingtype(T), pkeys(tgrp), columns(tgrp, value), cols)
end

function unstack(::Type{D}, ::Type{T}, key, val, cols::AbstractVector{S}) where {D <:Dataset, T, S}
    dest_val = Columns(Tuple(Array{Union{T, Missing}}(undef, length(val)) for i in cols); names = cols)
    for (i, el) in enumerate(val)
        for (k, v) in el
            ismissing(columns(dest_val, S(k))[i]) || error("Repeated values with same label are not allowed")
            columns(dest_val, S(k))[i] = v
        end
    end
    convert(collectiontype(D), key, dest_val)
end


text:
importall DataReader

fasta_record  = FastaRecord("test string1", "ACGT")
fasta_record2 = FastaRecord("test string1", "ACGT")

println("testing FastaRecord...")
@test fasta_record.description == "test string1"
@test fasta_record.sequence == "ACGT"

@test fasta_record == fasta_record2

println("testing FastaRecords loading from file...")

sequences = readSequences(dirname(@__FILE__()) * "/../data/input_test_sequences.faa")
@test length(sequences) > 0
@test length(sequences) == 10
@test sequences[1].description == "1"
@test sequences[1].sequence == "LWYVRMHMYR"
@test sequences[10].description == "10"
@test sequences[10].sequence == "QNHHVGNGPA"

println("testing ScoreMatrix loading from sample BLOSUM62 matrix file from ncbi ftp...")

matrix = readMatrix(dirname(@__FILE__()) * "/../data/blosum62.txt")
@test length(matrix.keys) == length(matrix.hsh) == 24

for s = 1 : length(matrix.hsh)
  @test length(matrix.hsh[s]) == s
end


text:
module TestNestedness
    using Test
    using EcologicalNetworks

    # Generate some data

    A = BipartiteProbabilisticNetwork([1.0 0.0 0.0; 0.0 0.1 0.0; 0.0 0.0 0.1])
    B = BipartiteProbabilisticNetwork([1.0 1.0 1.0; 1.0 0.1 0.0; 1.0 0.0 0.0])
    C = BipartiteProbabilisticNetwork([1.0 1.0 1.0; 1.0 0.1 0.3; 0.4 0.2 0.0])
    D = BipartiteProbabilisticNetwork([1.0 1.0 1.0; 0.0 0.1 1.0; 0.0 0.0 1.0])

    @test η(A) ≈ 0.0
    @test η(B) ≈ 1.0
    @test η(C) ≈ 0.9153846153846155
    @test η(D) ≈ 1.0
    @test_throws ArgumentError η(D; dims=3)

    # Almeida-Neto original example
    AN = BipartiteNetwork([1 0 1 1 1; 1 1 1 0 0; 0 1 1 1 0; 1 1 0 0 0; 1 1 0 0 0].==1)
    @test nodf(AN) ≈ 0.58 atol=0.01
    @test nodf(AN; dims=2) ≈ 0.63 atol=0.01
    @test nodf(AN; dims=1) ≈ 0.53 atol=0.01

    # Some WNODF examples
    WN1A = BipartiteQuantitativeNetwork([5 4 3 2 1; 4 3 2 1 0; 3 2 1 0 0; 2 1 0 0 0; 1 0 0 0 0])
    WN1B = BipartiteQuantitativeNetwork([5 4 3 2 1; 5 4 3 2 0; 5 4 3 0 0; 5 4 0 0 0; 5 0 0 0 0])
    WN1D = BipartiteQuantitativeNetwork([5 1 1 1 1; 1 0 0 0 0; 1 0 0 0 0; 1 0 0 0 0; 1 0 0 0 0])

    @test nodf(WN1A) ≈ 1.0 atol=0.01
    @test nodf(WN1A; dims=2) ≈ 1.0 atol=0.01
    @test nodf(WN1A; dims=1) ≈ 1.0 atol=0.01

    @test nodf(WN1B) ≈ 0.5 atol=0.01
    @test nodf(WN1B; dims=2) ≈ 1.0 atol=0.01
    @test nodf(WN1B; dims=1) ≈ 0.0 atol=0.01

    @test nodf(WN1D) ≈ 0.4 atol=0.01
    @test nodf(WN1D; dims=2) ≈ 0.4 atol=0.01
    @test nodf(WN1D; dims=1) ≈ 0.4 atol=0.01

end


text:
function jarquebera(varargin)
    #  JARQUEBERA MATLAB code for jarquebera.fig
    #       JARQUEBERA, by itself, creates a new JARQUEBERA or raises the existing
    #       singleton*.
    #
    #       H = JARQUEBERA returns the handle to a new JARQUEBERA or the handle to
    #       the existing singleton*.
    #
    #       JARQUEBERA('CALLBACK',hObject,eventData,handles,...) calls the local
    #       function named CALLBACK in JARQUEBERA.M with the given input arguments.
    #
    #       JARQUEBERA('Property','Value',...) creates a new JARQUEBERA or raises the
    #       existing singleton*.  Starting from the left, property value pairs are
    #       applied to the GUI before jarquebera_OpeningFcn gets called.  An
    #       unrecognized property name or invalid value makes property application
    #       stop.  All inputs are passed to jarquebera_OpeningFcn via varargin.
    #
    #       *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
    #       instance to run (singleton)".
    #
    #  See also: GUIDE, GUIDATA, GUIHANDLES

    #  Edit the above text to modify the response to help jarquebera

    #  Last Modified by GUIDE v2.5 12-Nov-2013 14:15:13

    #  Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @jarquebera_OpeningFcn, ...
                       'gui_OutputFcn',  @jarquebera_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end

    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    #  End initialization code - DO NOT EDIT
    return varargout

end

#  --- Executes just before jarquebera is made visible.
function jarquebera_OpeningFcn(hObject, ~, handles, varargin)
    #  This function has no output args, see OutputFcn.
    #  hObject    handle to figure
    #  eventdata  reserved - to be defined in a future version of MATLAB
    #  handles    structure with handles and user data (see GUIDATA)
    #  varargin   command line arguments to jarquebera (see VARARGIN)

    #  Choose default command line output for jarquebera
    handles.output = hObject;

    #  Update handles structure
    guidata(hObject, handles);

    #  UIWAIT makes jarquebera wait for user response (see UIRESUME)
    #  uiwait(handles.figure1);
    # JARQUE-BERA test for normality
    # Author:
    # Date: 29-June-2020

    # import XL data
    returns = getXLData('hfma_matlab_data.xlsx','CTA Index');

    sig = 0.95; # significance level (95# )
    dof = 2; # degrees of freedom
    [jb,critval,pval] = fJBTest(returns,sig,dof); # get Jarque-Bera test statistics

    # OUTPUT
    # format data
    data = num2cell([jb;sig*100;dof;critval;pval]);
    for i = 1:numel(data)
            data{i} = sprintf('# 12.4f',data{i});
    end
    # set column headings
    headings = {'Jarque-Bera Value','Sig. Level','Degrees of Freedom','Critical Value','p-Value'};
    # set table output
    set(handles.table_jarque_bera_test,'ColumnName','','RowName',headings,'Data',data);

end

#  --- Outputs from this function are returned to the command line.
function jarquebera_OutputFcn(~, ~, handles)
    #  varargout  cell array for returning output args (see VARARGOUT);
    #  hObject    handle to figure
    #  eventdata  reserved - to be defined in a future version of MATLAB
    #  handles    structure with handles and user data (see GUIDATA)

    #  Get default command line output from handles structure
    varargout{1} = handles.output;
    return varargout
end


text:
# MIRT.jl

"""
`MIRT` is the Michigan Image Reconstruction Toolbox
"""
module MIRT
	using Reexport

	@reexport using MIRTio # make I/O routines available

    const RealU = Number # Union{Real, Unitful.Length}
	include("z-all.jl")

end # module


text:
"""
    steady_analysis(surfaces, reference, freestream; kwargs...)

Perform a steady vortex lattice method analysis.  Return an object of type
[`System`](@ref) containing the system state.

# Arguments
 - `surfaces`:
   - Vector of grids of shape (3, nc+1, ns+1) which represent lifting surfaces
   or
   - Vector of matrices of shape (nc, ns) containing surface panels (see
   [`SurfacePanel`](@ref))
   where `nc` is the number of chordwise panels and `ns` is the number of
   spanwise panels
 - `reference`: Reference parameters (see [`Reference`](@ref))
 - `freestream`: Freestream parameters (see [`Freestream`](@ref))

# Keyword Arguments
 - `symmetric`: Flag for each surface indicating whether a mirror image across
    the X-Z plane should be used when calculating induced velocities. Defaults to
    `false` for each surface
 - `wakes`: Matrix of wake panels (see [`WakePanel`](@ref)) for each surface.  Each
    matrix has shape (nw, ns) where `nw` is the number of chordwise wake panels
    and `ns` is the number of spanwise panels for each surface, defaults to no
    wake panels for each surface
 - `nwake`: Number of chordwise wake panels to use from each wake in `wakes`,
    defaults to all wake panels for each surface
 - `surface_id`: Surface ID for each surface.  The finite core model is disabled
    when calculating the influence of surfaces/wakes that share the same ID.
    By default, all surfaces are assigned their own IDs
 - `wake_finite_core`: Flag for each wake indicating whether the finite core
    model should be enabled when calculating a wake's influence on itself and
    surfaces/wakes with the same surface ID.  Defaults to `true` for each surface.
 - `trailing_vortices`: Flags to enable/disable trailing vortices for each surface,
    defaults to `true` for each surface
 - `xhat`: Direction in which to shed trailing vortices, defaults to `[1, 0, 0]`
 - `additional_velocity`: Function which defines additional velocity as a
    function of location.
 - `fcore`: function which sets the finite core size for each surface based on
    the chord length and/or the panel width. Defaults to `(c, Δs) -> 1e-3`.
    Only used for grid inputs.
 - `calculate_influence_matrix`: Flag indicating whether the aerodynamic influence
    coefficient matrix has already been calculated.  Re-using the same AIC matrix
    will reduce calculation times when the underlying geometry has not changed.
    Defaults to `true`.  Note that this argument is only valid for the pre-allocated
    version of this function.
 - `near_field_analysis`: Flag indicating whether a near field analysis should be
    performed to obtain panel velocities, circulation, and forces. Defaults to `true`.
 - `derivatives`: Flag indicating whether the derivatives with respect
    to the freestream variables should be calculated. Defaults to `true`.
"""
function steady_analysis(surfaces, reference, freestream; kwargs...)

    # pre-allocate system storage
    system = System(surfaces)

    return steady_analysis!(system, surfaces, reference, freestream; kwargs...,
        calculate_influence_matrix = true)
end

"""
    steady_analysis!(system, surfaces, reference, freestream; kwargs...)

Pre-allocated version of `steady_analysis`.
"""
function steady_analysis!(system, surfaces, ref, fs;
    symmetric = fill(false, length(surfaces)),
    wakes = [Matrix{WakePanel{Float64}}(undef, 0, size(surfaces[i],2)) for i = 1:length(surfaces)],
    nwake = size.(wakes, 1),
    surface_id = 1:length(surfaces),
    wake_finite_core = fill(true, length(surfaces)),
    trailing_vortices = fill(true, length(surfaces)),
    xhat = SVector(1, 0, 0),
    additional_velocity = nothing,
    fcore = (c, Δs) -> 1e-3,
    calculate_influence_matrix = true,
    near_field_analysis = true,
    derivatives = true)

    # number of surfaces
    nsurf = length(surfaces)

    # check if input is a grid
    grid_input = typeof(surfaces) <: AbstractVector{<:AbstractArray{<:Any, 3}}

    # if only one value is provided, use it for all surfaces
    symmetric = isa(symmetric, Number) ? fill(symmetric, nsurf) : symmetric
    nwake = isa(nwake, Number) ? fill(nwake, nsurf) : nwake
    surface_id = isa(surface_id, Number) ? fill(surface_id, nsurf) : surface_id
    wake_finite_core = isa(wake_finite_core, Number) ? fill(wake_finite_core, nsurf) : wake_finite_core
    trailing_vortices = isa(trailing_vortices, Number) ? fill(trailing_vortices, nsurf) : trailing_vortices

    # update surface panels stored in `system`
    if grid_input
        # generate and store surface panels
        for isurf = 1:nsurf
            update_surface_panels!(system.surfaces[isurf], surfaces[isurf]; fcore)
        end
    else
        # store provided surface panels
        for isurf = 1:nsurf
            system.surfaces[isurf] .= surfaces[isurf]
        end
    end

    # update other parameters stored in `system`
    system.reference[] = ref
    system.freestream[] = fs
    system.symmetric .= symmetric
    system.wakes .= wakes
    system.nwake .= nwake
    system.surface_id .= surface_id
    system.wake_finite_core .= wake_finite_core
    system.trailing_vortices .= trailing_vortices
    system.xhat[] = xhat

    # unpack variables stored in `system`
    surfaces = system.surfaces
    AIC = system.AIC
    w = system.w
    Γ = system.Γ
    dw = system.dw
    dΓ = system.dΓ
    properties = system.properties
    dproperties = system.dproperties

    # see if wake panels are being used
    wake_panels = nwake .> 0

    # calculate/re-calculate AIC matrix (if necessary)
    if calculate_influence_matrix
        influence_coefficients!(AIC, surfaces;
            symmetric = symmetric,
            surface_id = surface_id,
            trailing_vortices = trailing_vortices .& .!wake_panels,
            xhat = xhat)
    end

    # calculate RHS
    if derivatives
        normal_velocity_derivatives!(w, dw, surfaces, wakes, ref, fs;
            additional_velocity = additional_velocity,
            Vcp = nothing, # no velocity at control points due to surface motion
            symmetric = symmetric,
            surface_id = surface_id,
            nwake = nwake,
            wake_finite_core = wake_finite_core,
            trailing_vortices = trailing_vortices,
            xhat = xhat)
    else
        normal_velocity!(w, surfaces, wakes, ref, fs;
            additional_velocity = additional_velocity,
            Vcp = nothing, # no velocity at control points due to surface motion
            symmetric = symmetric,
            surface_id = surface_id,
            nwake = nwake,
            wake_finite_core = wake_finite_core,
            trailing_vortices = trailing_vortices,
            xhat = xhat)
    end

    # solve for the circulation distribution
    if derivatives
        circulation_derivatives!(Γ, dΓ, AIC, w, dw)
    else
        circulation!(Γ, AIC, w)
    end

    if near_field_analysis
        # perform a near field analysis to obtain panel properties
        if derivatives
            near_field_forces_derivatives!(properties, dproperties, surfaces, wakes,
                ref, fs, Γ, dΓ;
                dΓdt = nothing, # no unsteady forces
                additional_velocity = additional_velocity,
                Vh = nothing, # no velocity due to surface motion
                Vv = nothing, # no velocity due to surface motion
                symmetric = symmetric,
                nwake = nwake,
                surface_id = surface_id,
                wake_finite_core = wake_finite_core,
                wake_shedding_locations = nothing, # shedding location at trailing edge
                trailing_vortices = trailing_vortices,
                xhat = xhat)
        else
            near_field_forces!(properties, surfaces, wakes, ref, fs, Γ;
                dΓdt = nothing, # no unsteady forces
                additional_velocity = additional_velocity,
                Vh = nothing, # no velocity due to surface motion
                Vv = nothing, # no velocity due to surface motion
                symmetric = symmetric,
                nwake = nwake,
                surface_id = surface_id,
                wake_finite_core = wake_finite_core,
                wake_shedding_locations = nothing, # shedding location at trailing edge
                trailing_vortices = trailing_vortices,
                xhat = xhat)
        end
    end

    # save flags indicating whether certain analyses have been performed
    system.near_field_analysis[] = near_field_analysis
    system.derivatives[] = derivatives

    # return the modified system
    return system
end

"""
    unsteady_analysis(surfaces, reference, freestream, dt; kwargs...)

Perform a unsteady vortex lattice method analysis.  Return an object of type
[`System`](@ref) containing the final system state, a matrix of surface panels
(see [`SurfacePanel`](@ref) for each surface at each time step, a matrix of surface
panel properties (see [`PanelProperties`](@ref)) for each surface at each time step,
and a matrix of wake panels (see [`WakePanel`](@ref)) for each surface at each time
step.

# Arguments
 - `surfaces`:
   - Grids of shape (3, nc+1, ns+1) which represent lifting surfaces
     or
   - Matrices of surface panels (see [`SurfacePanel`](@ref)) of shape
    (nc, ns)
    where `nc` is the number of chordwise panels and `ns` is the number of
    spanwise panels. Alternatively, a vector containing surface shapes/positions
    at each time step (including at `t=0`) may be provided to model
    moving/deforming lifting surfaces.
 - `reference`: Reference parameters (see [`Reference`](@ref))
 - `freestream`: Freestream parameters for each time step (see [`Freestream`](@ref))
 - `dt`: Time step vector

# Keyword Arguments
 - `symmetric`: Flags indicating whether a mirror image (across the X-Z plane) should
    be used when calculating induced velocities, defaults to `false` for each surface
 - `initial_wakes`: Vector of initial wakes corresponding to each surface, represented
    by matrices of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns) where
    `nw` is the number of chordwise wake panels and `ns` is the number of
    spanwise panels. Defaults to no wake panels for each surface
 - `initial_circulation`: Vector containing the initial circulation of all surface
    panels in the system.  Defaults to `zeros(N)` where `N` is the total number
    of surface panels in `surfaces`.
 - `nwake`: Maximum number of wake panels in the chordwise direction for each
    surface.  Defaults to `length(dx)` for all surfaces.
 - `surface_id`: Surface ID for each surface.  The finite core model is disabled
    when calculating the influence of surfaces/wakes that share the same ID.
 - `wake_finite_core`: Flag for each wake indicating whether the finite core
    model should be enabled when calculating the wake's influence on itself and
    surfaces/wakes with the same surface ID.  Defaults to `true` for each surface.
 - `save`: Time indices at which to save the time history, defaults to `1:length(dx)`
 - `calculate_influence_matrix`: Flag indicating whether the aerodynamic influence
    coefficient matrix needs to be calculated.  Re-using the same AIC matrix
    will (slightly) reduce calculation times when the underlying geometry has not changed.
    Defaults to `true`.  Note that this argument only affects the pre-allocated
    version of this function.
 - `near_field_analysis`: Flag indicating whether a near field analysis should be
    performed to obtain panel velocities, circulation, and forces for the final
    time step. Defaults to `true`.
 - `derivatives`: Flag indicating whether the derivatives with respect to the
    freestream variables should be calculated for the final time step. Defaults
    to `true`.
"""
unsteady_analysis

# same grids/surfaces at each time step
function unsteady_analysis(surfaces::AbstractVector{<:AbstractArray}, ref, fs, dt;
    nwake = fill(length(dt), length(surfaces)), kwargs...)

    # pre-allocate system storage
    system = System(surfaces; nw = nwake)

    return unsteady_analysis!(system, surfaces, ref, fs, dt;
        kwargs..., nwake, calculate_influence_matrix = true)
end

# different grids/surfaces at each time step
function unsteady_analysis(surfaces::AbstractVector{<:AbstractVector{<:AbstractArray}},
    ref, fs, dt; nwake = fill(length(dt), length(surfaces[1])), kwargs...)

    # pre-allocate system storage
    system = System(surfaces[1]; nw = nwake)

    return unsteady_analysis!(system, surfaces, ref, fs, dt;
        kwargs..., nwake, calculate_influence_matrix = true)
end

"""
    unsteady_analysis!(system, surfaces, reference, freestream, dt; kwargs...)

Pre-allocated version of `unsteady_analysis`.
"""
function unsteady_analysis!(system, surfaces, ref, fs, dt;
    symmetric = fill(false, length(system.surfaces)),
    initial_circulation = zero(system.Γ),
    initial_wakes = [Matrix{WakePanel{Float64}}(undef, 0, size(surfaces[i], 2)) for i = 1:length(system.surfaces)],
    nwake = fill(length(dt), length(system.surfaces)),
    surface_id = 1:length(system.surfaces),
    wake_finite_core = fill(true, length(system.surfaces)),
    additional_velocity = nothing,
    fcore = (c, Δs) -> 1e-3,
    save = 1:length(dt),
    calculate_influence_matrix = true,
    near_field_analysis = true,
    derivatives = true)

    # float number type
    TF = eltype(system)

    # number of surfaces
    nsurf = length(system.surfaces)

    # surface motion?
    surface_motion = eltype(surfaces) <: AbstractVector

    # --- Input Pre-Processing --- #

    # convert scalar inputs to appropriately sized vectors
    symmetric = isa(symmetric, Number) ? fill(symmetric, nsurf) : symmetric
    nwake = isa(nwake, Number) ? fill(nwake, nsurf) : nwake
    surface_id = isa(surface_id, Number) ? fill(surface_id, nsurf) : surface_id
    wake_finite_core = isa(wake_finite_core, Number) ? fill(wake_finite_core, nsurf) : wake_finite_core
    fs = isa(fs, Freestream) ? fill(fs, length(dt)) : fs

    # extract initial surface panels from input
    if surface_motion
        # surface moves, store initial surface panel locations
        initial_surfaces = surfaces[1]
    else
        # surface doesn't move
        initial_surfaces = surfaces
    end

    # --- Update System Parameters --- #

    system.reference[] = ref
    system.symmetric .= symmetric
    system.surface_id .= surface_id
    system.wake_finite_core .= wake_finite_core
    system.trailing_vortices .= false

    # check if existing wake panel storage is sufficient, replace if necessary
    for isurf = 1:nsurf
        if size(system.wakes[isurf], 1) < nwake[isurf]
            # get surface/wake dimensions
            nc, ns = size(surfaces[isurf])
            nw = nwake[isurf]

            # update wake panel storage
            system.wakes[isurf] = Matrix{WakePanel{TF}}(undef, nwake[isurf], size(surfaces[isurf], 2))
        end
    end

    # find intersecting surfaces
    repeated_points = repeated_trailing_edge_points(initial_surfaces)

    # --- Set Initial Simulation Variables --- #

    # store initial surface panels in `system`
    if eltype(initial_surfaces) <: AbstractArray{<:Any, 3}
        # initial surfaces are input as a grid, convert to surface panels
        for isurf = 1:nsurf
            update_surface_panels!(system.surfaces[isurf], initial_surfaces[isurf]; fcore)
        end
    else
        # initial surfaces are input as matrices of surface panels
        for isurf = 1:nsurf
            system.surfaces[isurf] .= initial_surfaces[isurf]
        end
    end

    # store initial wake panels in `system`
    for isurf = 1:nsurf
        for I in CartesianIndices(initial_wakes[isurf])
            system.wakes[isurf][I] = initial_wakes[isurf][I]
        end
    end

    # store initial freestream parameters in `system`
    system.freestream[] = fs[1]

    # store initial circulation parameters in `system`
    system.Γ .= initial_circulation

    # set the initial number of wake panels for each surface
    iwake = [min(size(initial_wakes[isurf], 1), nwake[isurf]) for isurf = 1:nsurf]

    # --- Begin Simulation --- #

    # initialize solution history for each time step
    surface_history = Vector{Vector{Matrix{SurfacePanel{TF}}}}(undef, length(save))
    property_history = Vector{Vector{Matrix{PanelProperties{TF}}}}(undef, length(save))
    wake_history = Vector{Vector{Matrix{WakePanel{TF}}}}(undef, length(save))
    isave = 1

    # # loop through all time steps
    for it = 1 : length(dt)

        first_step = it == 1
        last_step = it == length(dt)

        if surface_motion
            propagate_system!(system, surfaces[1+it], fs[it], dt[it];
                additional_velocity,
                repeated_points,
                nwake = iwake,
                eta = 0.25,
                calculate_influence_matrix = first_step && calculate_influence_matrix,
                near_field_analysis = it in save || (last_step && near_field_analysis),
                derivatives = last_step && derivatives)
        else
            propagate_system!(system, fs[it], dt[it];
                additional_velocity,
                repeated_points,
                nwake = iwake,
                eta = 0.25,
                calculate_influence_matrix = first_step && calculate_influence_matrix,
                near_field_analysis = it in save || (last_step && near_field_analysis),
                derivatives = last_step && derivatives)
        end

        # increment wake panel counter for each surface
        for isurf = 1:nsurf
            if iwake[isurf] < nwake[isurf]
                iwake[isurf] += 1
            end
        end

        # save the surface shape, properties, and resulting shed wake
        if it in save
            surface_history[isave] = [copy(system.surfaces[isurf]) for isurf = 1:nsurf]
            property_history[isave] = [copy(system.properties[isurf]) for isurf = 1:nsurf]
            wake_history[isave] = [system.wakes[isurf][1:iwake[isurf], :] for isurf = 1:nsurf]
            isave += 1
        end

    end

    # return the modified system and time history
    return system, surface_history, property_history, wake_history
end

"""
    propagate_system!(system, [surfaces, ] freestream, dt; kwargs...)

Propagate the state variables in `system` forward one time step using the
unsteady vortex lattice method system of equations.

# Arguments
 - `system`: Object of type `system` which contains the current system state
 - `surfaces`: Surface locations at the end of this time step. If omitted,
   surfaces are assumed to be stationary.
 - `freestream`: Freestream parameters corresponding to this time step.
 - `dt`: Time increment

# Keyword Arguments
 - `additional_velocity`: Function which defines additional velocity as a
    function of location.
 - `repeated_points`: Dictionary of the form `Dict((isurf, i) => [(jsurf1, j1),
    (jsurf2, j2)...]` which defines repeated trailing edge points.  Trailing edge
    point `i` on surface `isurf` is repeated on surface `jsurf1` at point `j1`,
    `jsurf2` at point `j2`, and so forth. See [`repeated_trailing_edge_points`](@ref)
 - `nwake`: Number of wake panels in the chordwise direction for each surface.
 - `eta`: Time step fraction used to define separation between trailing
    edge and wake shedding location.  Typical values range from 0.2-0.3.
 - `calculate_influence_matrix`: Flag indicating whether the aerodynamic influence
    coefficient matrix needs to be calculated. If argument `surfaces` is provided
    the influence matrix will always be recalculated.
 - `near_field_analysis`: Flag indicating whether a near field analysis should be
    performed to obtain panel velocities, circulation, and forces.
 - `derivatives`: Flag indicating whether the derivatives with respect to the
    freestream variables should be calculated.
"""
propagate_system!

# stationary surfaces
propagate_system!(system, fs, dt; kwargs...) = propagate_system!(system,
    nothing, fs, dt; kwargs...)

# moving/deforming surfaces
function propagate_system!(system, surfaces, fs, dt;
    additional_velocity,
    repeated_points,
    nwake,
    eta,
    calculate_influence_matrix,
    near_field_analysis,
    derivatives)

    # NOTE: Each step models the transition from `t = t[it]` to `t = [it+1]`
    # (e.g. the first step models from `t = 0` to `t = dt`).  Properties are
    # assumed to be constant during each time step and resulting transient
    # forces/moments are provided at `t[it+1]` (e.g. the properties returned
    # for the first step correspond to `t = dt`).

    nsurf = length(system.surfaces)

    # unpack constant system parameters
    ref = system.reference[]
    symmetric = system.symmetric
    surface_id = system.surface_id
    wake_finite_core = system.wake_finite_core
    trailing_vortices = system.trailing_vortices
    xhat = system.xhat[]

    # unpack system storage (including state variables)
    previous_surfaces = system.previous_surfaces
    current_surfaces = system.surfaces
    properties = system.properties
    dproperties = system.dproperties
    wakes = system.wakes
    wake_velocities = system.V
    wake_shedding_locations = system.wake_shedding_locations
    AIC = system.AIC
    w = system.w
    dw = system.dw
    Γ = system.Γ
    dΓ = system.dΓ
    dΓdt = system.dΓdt
    Vcp = system.Vcp
    Vh = system.Vh
    Vv = system.Vv
    Vte = system.Vte

    # check if the surfaces are moving
    surface_motion = !isnothing(surfaces)

    # move geometry and calculate velocities for this time step
    if surface_motion

        # check if grids are used to represent the new surfaces
        grid_input = eltype(surfaces) <: AbstractArray{<:Any, 3}

        for isurf = 1:nsurf

            # save current surfaces as previous surfaces
            previous_surfaces[isurf] .= current_surfaces[isurf]

            # set new surface shape...
            if grid_input
                # ...based on grid inputs
                update_surface_panels!(current_surfaces[isurf], surfaces[isurf]; fcore)
            else
                # ...based on surface panels
                current_surfaces[isurf] .= surfaces[isurf]
            end

            # calculate surface velocities
            get_surface_velocities!(Vcp[isurf], Vh[isurf], Vv[isurf], Vte[isurf],
                current_surfaces[isurf], previous_surfaces[isurf], dt)
        end
    else
        # zero out surface motion
        for isurf = 1:nsurf
            system.Vcp[isurf] .= Ref(SVector(0.0, 0.0, 0.0))
            system.Vh[isurf] .= Ref(SVector(0.0, 0.0, 0.0))
            system.Vv[isurf] .= Ref(SVector(0.0, 0.0, 0.0))
            system.Vte[isurf] .= Ref(SVector(0.0, 0.0, 0.0))
        end
    end

    # update stored freestream parameters for this time step
    system.freestream[] = fs

    # update number of wake panels for each surface for this time step
    system.nwake .= nwake

    # update the wake shedding location for this time step
    update_wake_shedding_locations!(wakes, wake_shedding_locations,
        current_surfaces, ref, fs, dt, additional_velocity, Vte,
        nwake, eta)

    # calculate/re-calculate AIC matrix (if necessary)
    if surface_motion || calculate_influence_matrix
        influence_coefficients!(AIC, current_surfaces;
            symmetric = symmetric,
            wake_shedding_locations = wake_shedding_locations,
            surface_id = surface_id,
            trailing_vortices = trailing_vortices,
            xhat = xhat)
    end

    # update the AIC matrix to use the new wake shedding locations
    update_trailing_edge_coefficients!(AIC, current_surfaces;
        symmetric = symmetric,
        wake_shedding_locations = wake_shedding_locations,
        trailing_vortices = trailing_vortices)

    # calculate RHS
    if derivatives
        normal_velocity_derivatives!(w, dw, current_surfaces, wakes,
            ref, fs; additional_velocity, Vcp, symmetric, nwake,
            surface_id, wake_finite_core, trailing_vortices, xhat)
    else
        normal_velocity!(w, current_surfaces, wakes, ref, fs;
            additional_velocity, Vcp, symmetric, nwake, surface_id,
            wake_finite_core, trailing_vortices, xhat)
    end

    # save (negative) previous circulation in dΓdt
    dΓdt .= .-Γ

    # solve for the new circulation
    if derivatives
        circulation_derivatives!(Γ, dΓ, AIC, w, dw)
    else
        circulation!(Γ, AIC, w)
    end

    # solve for dΓdt using finite difference `dΓdt = (Γ - Γp)/dt`
    dΓdt .+= Γ # add newly computed circulation
    dΓdt ./= dt # divide by corresponding time step

    # compute transient forces on each panel (if necessary)
    if near_field_analysis
        if derivatives
            near_field_forces_derivatives!(properties, dproperties,
                current_surfaces, wakes, ref, fs, Γ, dΓ; dΓdt,
                additional_velocity, Vh, Vv, symmetric, nwake,
                surface_id, wake_finite_core, wake_shedding_locations,
                trailing_vortices, xhat)
        else
            near_field_forces!(properties, current_surfaces, wakes,
                ref, fs, Γ; dΓdt, additional_velocity, Vh, Vv,
                symmetric, nwake, surface_id, wake_finite_core,
                wake_shedding_locations, trailing_vortices, xhat)
        end

        # save flag indicating that a near-field analysis has been performed
        system.near_field_analysis[] = near_field_analysis
    end

    # save flag indicating that derivatives wrt freestream variables have been obtained
    system.derivatives[] = derivatives

    # update wake velocities
    get_wake_velocities!(wake_velocities, current_surfaces,
        wakes, ref, fs, Γ, additional_velocity, Vte, symmetric,
        repeated_points, nwake, surface_id, wake_finite_core,
        wake_shedding_locations, trailing_vortices, xhat)

    # shed additional wake panel (and translate existing wake panels)
    shed_wake!(wakes, wake_shedding_locations, wake_velocities,
        dt, current_surfaces, Γ, nwake)

    return system
end


text:
include("FoggyMountainObj.jl")

function main()
# One BFGS run with poor starting values
theta = [8.0, -8.0]
thetahat, obj_value, converged = fminunc(FoggyMountainObj, theta)
println()
println("The result with poor start values")
println("objective function value: ", round(obj_value,digits=5))
println("local minimizer: ", thetahat)
println()

# Now try simulated annealing
lb = [-20.0,-20.0]
ub = -lb
xopt = samin(FoggyMountainObj, theta, lb, ub, verbosity=1)

# Now try BFGS with multiple start values
bestobj = 1000
thetahat = zeros(2)
for i = 1:50
    theta = (ub-lb).*rand(2) + lb # use same bounds as SA
    theta, obj_value, converged = fminunc(FoggyMountainObj, theta)
    if obj_value < bestobj
        thetahat = theta
        bestobj = obj_value
    end
end
println("the estimate using BFGS with multiple start values: ", thetahat)
println("the best objective function value: ", bestobj)
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
login(; kwargs...) = wandb.login(; kwargs...)

mutable struct WandbLogger <: AbstractLogger
    wrun::PyObject
    step_increment::Int
    global_step::Int
    min_level::LogLevel
end

function WandbLogger(; project, name=nothing, min_level=Info, step_increment=1, start_step=0, kwargs...)
    wrun = nothing
    @static if Sys.iswindows()
        if :settings in keys(kwargs)
            wrun = wandb.init(; project=project, name=name, kwargs...)
        else
            wrun = wandb.init(; project=project, name=name, settings=wandb.Settings(; start_method="thread"), kwargs...)
        end
    else
        wrun = wandb.init(; project=project, name=name, kwargs...)
    end
    if !isnothing(name) && wrun.name != name
        @warn "There is an ongoing wandb run. Please `close` the run before initializing a new one."
    end
    return WandbLogger(wrun, step_increment, start_step, min_level)
end

function Base.show(io::IO, lg::WandbLogger)
    str = "WandbLogger(\"$(lg.wrun.project)\", \"$(lg.wrun.name)\", " *
          "id=$(lg.wrun.id), min_level=$(lg.min_level), " *
          "current_step=$(lg.global_step))"
    Base.print(io, str)
    return Base.printstyled(io, " @ $(lg.wrun.url)"; color=:yellow)
end

increment_step!(lg::WandbLogger, Δ_Step) = lg.global_step += Δ_Step

# https://docs.wandb.ai/guides/track/log
log(lg::WandbLogger, logs::Dict; kwargs...) = lg.wrun.log(logs; kwargs...)

# https://docs.wandb.ai/guides/track/config
update_config!(lg::WandbLogger, dict::Dict; kwargs...) = lg.wrun.config.update(dict; kwargs...)

# https://docs.wandb.ai/ref/python/finish
Base.close(lg::WandbLogger; kwargs...) = lg.wrun.finish(; kwargs...)

function finish(lg::WandbLogger; kwargs...)
    @warn "The desired way to end the run is `close`. Please update your code to use the same" maxlog = 1
    return close(lg; kwargs...)
end

# https://docs.wandb.ai/ref/python/save
save(lg::WandbLogger, file::String; kwargs...) = lg.wrun.save(file; kwargs...)

get_config(lg::WandbLogger) = lg.wrun.config
get_config(lg::WandbLogger, key::String) = get(lg.wrun.config, key)

# Logging Types: Image, Video, Histogram, 3D Objects, PR Curve
## We assume the images to be ordered as per Flux conventions
Image(img::AbstractArray{T,3}; kwargs...) where {T} = wandb.Image(permutedims(img, (3, 1, 2)); kwargs...)
Image(img::AbstractMatrix; kwargs...) = wandb.Image(img; kwargs...)
Image(img::String; kwargs...) = wandb.Image(img; kwargs...)

Video(vid::AbstractArray{T,4}; kwargs...) where {T} = wandb.Video(permutedims(vid, (4, 3, 2, 1)); kwargs...)
Video(vid::AbstractArray{T,3}; kwargs...) where {T} = wandb.Video(permutedims(vid, (3, 2, 1)); kwargs...)
Video(vid::String; kwargs...) = wandb.Video(vid; kwargs...)

Histogram(x::AbstractArray; kwargs...) = wandb.Histogram(x; kwargs...)

Object3D(path::String) = wandb.Object3D(open(path, 'r'))

function precision_recall(y_test::AbstractVector, y_probs::AbstractVector, labels::AbstractVector)
    return wandb.plots.precision_recall(y_test, y_probs, labels)
end


text:
"""
	SMMOptions

SMMOptions is a struct that contains options related to the optimization.
"""
struct SMMOptions
	globalOptimizer::Symbol #algorithm for finding a global maximum
	localOptimizer::Symbol 	#algorithm for finding a local maximum
	maxFuncEvals::Int64			#maximum number of evaluations
	saveSteps::Int64				#maximum number of steps
	saveName::String				#name under which the optimization should be saved
	showDistance::Bool			#show the distance, everytime the objective function is calculated?
	minBox::Bool						#When looking for a local maximum, use Fminbox ?
	populationSize::Int64		#When using BlackBoxOptim, set the population size
	penaltyValue::Float64 	#Objective function's value when the model fails
	gridType::Symbol				#sampling procedure to use (latin hypercube by default)
	saveStartingValues::Bool #whether or not saving the starting values when using local_to_global
	maxTrialsStartingValues::Int64 #maximum number of attempts when searching for valid starting values
	thresholdStartingValue::Float64 #value under which a point is considered as a valid starting value
end

function SMMOptions( ;globalOptimizer::Symbol=:dxnes,
											localOptimizer::Symbol=:LBFGS,
											maxFuncEvals::Int64=1000,
											saveSteps::Int64 = 100,
											saveName::String = get_now(),
											showDistance::Bool = false,
											minBox::Bool = false,
											populationSize::Int64 = 50,
											penaltyValue::Float64 = 999999.0,
											gridType::Symbol = :latin,
											saveStartingValues::Bool = true,
											maxTrialsStartingValues::Int64 = 1000,
											thresholdStartingValue::Float64 = 99999.0)

	# Safety Checks
	#--------------
	if saveSteps == 0
		error("You cannot set saveSteps equal to 0.")
	end

	if saveSteps > maxFuncEvals
		error("Error in the constructor for SMMOptions. \n saveSteps = $(saveSteps) > maxFuncEvals = $(maxFuncEvals)")
	end

	if mod(maxFuncEvals, saveSteps) != 0
		error("Error in the constructor for SMMOptions. \n maxFuncEvals should be a multiple of saveSteps")
	end

	if thresholdStartingValue > penaltyValue
		error("Please set thresholdStartingValue < penaltyValue.")
	end

	# Warnings
	#----------
	# Saving too often (relative to population size) leads to a poor performance if using BlackBoxOptim
	# see here: https://github.com/robertfeldt/BlackBoxOptim.jl/issues/99
	#--------------------------------------------------------------------
	# [TODO] Give a warning relative to the Population Size instead.
	if saveSteps/populationSize < 6 && is_bb_optimizer(globalOptimizer) == true
		info("WARNING. When saveSteps/populationSize < 6 using BlackBoxOptim, the performance of the global maximizer may deteriorate.")
	end

	SMMOptions(globalOptimizer,
						localOptimizer,
						maxFuncEvals,
						saveSteps,
						saveName,
						showDistance,
						minBox,
						populationSize,
						penaltyValue,
						gridType,
						saveStartingValues,
						maxTrialsStartingValues,
						thresholdStartingValue)

end

"""
	SMMProblem

SMMProblem is a mutable struct that caries all the information needed to
perform the optimization and display the results.
"""
mutable struct SMMProblem
	iter::Int64
	priors::OrderedDict{String,Array{Float64,1}}
	empiricalMoments::OrderedDict{String,Array{Float64,1}}
	simulatedMoments::OrderedDict{String, Float64}
	distanceEmpSimMoments::Float64
	simulate_empirical_moments::Function					#returns an ordered dict
	simulate_empirical_moments_array::Function		#returns an Array
	objective_function::Function
	options::SMMOptions
	bbSetup::BlackBoxOptim.OptController					#set up when using BlackBoxOptim (global minimum)
	bbResults::BlackBoxOptim.OptimizationResults	#results when using BlackBoxOptim (global minimum)
	optimResults::Optim.OptimizationResults				#results when using Optim (local minimum)
	Sigma0::Array{Float64,2}											#distance matrix (in the terminology of Duffie and Singleton (1993))
	Avar::Array{Float64,2}											  #asymptotic variance of the SMM estimate
end

# Constructor for SMMProblem
#------------------------------------------------------------------------------
function SMMProblem(  ;iter::Int64 = 0,
											priors::OrderedDict{String,Array{Float64,1}} = OrderedDict{String,Array{Float64,1}}(),
											empiricalMoments::OrderedDict{String,Array{Float64,1}} = OrderedDict{String,Array{Float64,1}}(),
											simulatedMoments::OrderedDict{String, Float64} = OrderedDict{String,Float64}(),
											distanceEmpSimMoments::Float64 = 0.,
											simulate_empirical_moments::Function = default_function,       #returns an ordered dict
											simulate_empirical_moments_array::Function = default_function, #returns an Array
											objective_function::Function = default_function,
											options::SMMOptions = SMMOptions(),
											bbSetup::BlackBoxOptim.OptController = defaultbbOptimOptController,
											bbResults::BlackBoxOptim.OptimizationResults = defaultbbOptimOptimizationResults,
											optimResults::Optim.OptimizationResults = defaultOptimResults,
											Sigma0::Array{Float64,2} = Array{Float64}(0,0),
											Avar::Array{Float64,2} = Array{Float64}(0,0))

	SMMProblem(iter,
						priors,
						empiricalMoments,
						simulatedMoments,
						distanceEmpSimMoments,
						simulate_empirical_moments,
						simulate_empirical_moments_array,
						objective_function,
						options,
						bbSetup,
						bbResults,
						optimResults,
						Sigma0,
						Avar)

end

"""
	default_function(x)

Function x->x. Used to initialize functions.
"""
function default_function(x)
	println("default_function, returns input")
	x
end

"""
	rosenbrock2d(x)

Rosenbrock function. Used to initialize BlackBoxOptim.OptController().
"""
function rosenbrock2d(x)
  return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

# It is quite useful to have "default" BlackBoxOptim.OptController and
# BlackBoxOptim.OptimizationResults objects, since
# BlackBoxOptim.OptController() and BlackBoxOptim.OptimizationResults()
# do not work
#-------------------------------------------------------------------------------
defaultbbOptimOptController = bbsetup(x -> rosenbrock2d(x);
																		Method=:dxnes,
																		SearchRange = (-5.0, 5.0),
             												NumDimensions = 2, MaxFuncEvals = 2,
																		TraceMode = :silent)

defaultbbOptimOptimizationResults = bboptimize(defaultbbOptimOptController)

defaultOptimResults = optimize(rosenbrock2d, [0.0, 0.0], LBFGS())

"""
	is_global_optimizer(s::Symbol)

function to check that the global optimizer chosen is available.
"""
function is_global_optimizer(s::Symbol)

	# Global Optimizers using BlackBoxOptim
	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	#------------------------------------------------------------------------------
	listValidGlobalOptimizers = [:dxnes, :adaptive_de_rand_1_bin_radiuslimited, :xnes,
															 :de_rand_1_bin_radiuslimited, :adaptive_de_rand_1_bin,
															 :generating_set_search, :de_rand_1_bin,
															 :separable_nes, :resampling_inheritance_memetic_search,
															 :probabilistic_descent, :resampling_memetic_search,
															 :de_rand_2_bin_radiuslimited, :de_rand_2_bin,
															 :random_search, :simultaneous_perturbation_stochastic_approximation]

	in(s, listValidGlobalOptimizers)

end

"""
	is_bb_optimizer(s::Symbol)

function to check whether the global optimizer is using BlackBoxOptim
"""
function is_bb_optimizer(s::Symbol)

	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	listbbOptimizers = [:dxnes, :adaptive_de_rand_1_bin_radiuslimited, :xnes,
															 :de_rand_1_bin_radiuslimited, :adaptive_de_rand_1_bin,
															 :generating_set_search, :de_rand_1_bin,
															 :separable_nes, :resampling_inheritance_memetic_search,
															 :probabilistic_descent, :resampling_memetic_search,
															 :de_rand_2_bin_radiuslimited, :de_rand_2_bin,
															 :random_search, :simultaneous_perturbation_stochastic_approximation]

	in(s, listbbOptimizers)

end

"""
	is_local_optimizer(s::Symbol)

function to check that the local optimizer chosen is available.
"""
function is_local_optimizer(s::Symbol)

	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	listValidLocalOptimizers = [:NelderMead, :SimulatedAnnealing, :ParticleSwarm,
															:BFGS, :LBFGS, :ConjugateGradient, :GradientDescent,
															:MomentumGradientDescent, :AcceleratedGradientDescent]

	in(s, listValidLocalOptimizers)

end

"""
	is_optim_optimizer(s::Symbol)

function to check whether the local minimizer uses the package Optim.
"""
function is_optim_optimizer(s::Symbol)

	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	listOptimOptimizers = [:NelderMead, :SimulatedAnnealing, :ParticleSwarm,
															:BFGS, :LBFGS, :ConjugateGradient, :GradientDescent,
															:MomentumGradientDescent, :AcceleratedGradientDescent]

	in(s, listOptimOptimizers)

end

"""
	convert_to_optim_algo(s::Symbol)

function to convert local optimizer (of type Symbol) to an Optim algo.
"""
function convert_to_optim_algo(s::Symbol)
	eval(Meta.parse("$(s)()"))
end

"""
	convert_to_fminbox(s::Symbol)

function to convert local optimizer (of type Symbol) to a Fminbox usable
by Optim.
"""
function convert_to_fminbox(s::Symbol)

	# Old API (before v0.15.0)
	# To be changed when switching to Julia v0.7
	eval(Meta.parse("Fminbox{$(s)}()"))

end


text:
if lowercase(get(ENV, "CONCURRENTCOLLECTIONS_JL_ASSERTION", "false")) == "true"
    import ConcurrentCollections
    ConcurrentCollections.Implementations.enable_assertion()
    @assert ConcurrentCollections.Implementations.assertion_enabled()
    @info "ConcurrentCollections: Assertion enabled"
else
    @info "ConcurrentCollections: Assertion disenabled (default)"
end

using TestFunctionRunner
TestFunctionRunner.@run(paths = ["../benchmark/ConcurrentCollectionsBenchmarks"])


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterfaceTapConfiguration <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function NetworkInterfaceTapConfiguration(;id=nothing)
        o = new()
        validate_property(NetworkInterfaceTapConfiguration, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type NetworkInterfaceTapConfiguration

const _property_map_NetworkInterfaceTapConfiguration = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_NetworkInterfaceTapConfiguration = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ NetworkInterfaceTapConfiguration }) = collect(keys(_property_map_NetworkInterfaceTapConfiguration))
Swagger.property_type(::Type{ NetworkInterfaceTapConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkInterfaceTapConfiguration[name]))}
Swagger.field_name(::Type{ NetworkInterfaceTapConfiguration }, property_name::Symbol) =  _property_map_NetworkInterfaceTapConfiguration[property_name]

function check_required(o::NetworkInterfaceTapConfiguration)
    true
end

function validate_property(::Type{ NetworkInterfaceTapConfiguration }, name::Symbol, val)
end


text:
using Documenter, DICE

DocMeta.setdocmeta!(DICE, :DocTestSetup, :(using DICE); recursive=true)

makedocs(
    modules = [DICE],
    doctest = true,
    format = Documenter.HTML(prettyurls=!("local" in ARGS)),
    sitename = "DICE.jl",
    authors = "Tim DuBois",
    linkcheck = false,
    pages = [
        "Home" => "index.md",
        "API" => "api.md",
    ],
)

deploydocs(
    repo   = "github.com/Libbum/DICE.jl.git",
    target = "build",
)


text:
module Experimental
reload("OnlineStats")
using OnlineStats
O = OnlineStats

using DataGenerator, GLM

# For testing LogRegMM
x, y, β = logregdata(1_000_000, 5, false)
@time o = O.LogRegMM(x, y, LearningRate(.7))
o2 = StatLearn(x, y, LearningRate(.7), LogisticRegression(), AdaGrad(); intercept = false)
o3 = glm(x, y, Binomial())
@show coef(o)
@show coef(o2)
@show coef(o3)

@show loss(o, x, y)
@show loss(o2, x, y)
# # For Testing QuantRegMM
# x, y, β = linregdat(100_000, 5)
# o = QuantRegMM(x, y, LearningRate(.5))
# o2 = StatLearn(x, y, LearningRate(.5), QuantileRegression(.8))
# @show coef(o)
end


text:
@testset "Sphere" begin

    # Objective function
    sphere(x::AbstractVector) = sum(x.*x)

    # Parameters
    N = 30
    P = 25
    initial = ones(N)

    function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::ES, options)
        record["fitpop"] = state.fitness
        record["σ"] = strategy(state).σ
    end

    function Evolutionary.terminate(state::Evolutionary.ESState)
        strategy(state).σ < 1e-10
    end

    # Testing: (μ/μ_I, λ)-σ-Self-Adaptation-ES
    # with isotropic mutation operator y' := y + σ(N_1(0, 1), ..., N_N(0, 1))
    result = Evolutionary.optimize(
        sphere,
        initial,
        ES(
            initStrategy = IsotropicStrategy(N),
            recombination = average, srecombination = average,
            mutation = gaussian, smutation = gaussian,
            selection=:comma,
            μ = 3, λ = P
        ), Evolutionary.Options(show_trace=false,iterations=1000));
    # show(result)
    println("(3/3,$(P))-σ-SA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test minimum(result) ≈ 0.0 atol=1e-3
    @test sum(x->x.^2, Evolutionary.minimizer(result)) ≈ 0.0 atol=1e-3
    @test length(Evolutionary.minimizer(result)) == N

    # disable custom functions
    Evolutionary.terminate(state::Evolutionary.ESState) = false
    Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::ES, options) = ()

    # Testing: GA
    result = Evolutionary.optimize(
        sphere,
        initial,
        GA(
            populationSize = 4P,
            mutationRate = 0.15,
            ɛ = 0.1,
            selection = susinv,
            crossover = intermediate(0.25),
            mutation = domainrange(fill(0.5,N)),
        ));
    # show(result)
    println("GA:INTER:DOMRNG:(N=$(N), P=$(P)) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test minimum(result) ≈ 0.0 atol=1e-3
    @test sum(x->x.^2, Evolutionary.minimizer(result)) ≈ 0.0 atol=1e-3
    @test length(Evolutionary.minimizer(result)) == N

end


text:
ishidden = s -> startswith(basename(s),'.')
isimg = s -> !ishidden(s) && any(endswith.(lowercase(s),[".tif",".tiff",".png",".jpg",".jpeg"]))

"""
	interpimage (img)
Returns a callable function that performs interpolation/extrapolation of a filtered image.
"""
function interpimage(img)
	X = interpolate(img, BSpline(Linear()))
	Z = extrapolate(X,0)  # extend by zero
	return (i,j) -> float(Z(i,j))
end

"Return data relevant for cornea detection to be used on the given image."
function detectiondata(img,m,n;options=get_defaults())
	# select the angles for optimization of the circle residual
	θ = 2π*(-180:180)/360
	select = falses(size(θ))
	for rng in options.angle_ranges
		select[ rng[1] .<= θ .<= rng[2] ] .= true
	end
	θ = θ[select]

	X = imresize(img,m,n)
	purk = findpurkinje(X)

	# optimization initializations
	u_init = [ (m/size(X,1)).*i for i in initvals(X,purk) ]  # "smart"
	append!(u_init,options.initializer(m,n))            # dumb

	G = green.(X)

	# try to screen out the eyelid lines, using bright blue values as the indicator.
	C = blue.(X)
	Cmax = mapwindow(maximum,C,(15,15))  # windowed max
	thresh = 0.75*maximum(C)
	G[@. (Cmax > thresh) & (G > 0.25) ] .= 0.33

	# Screen out purkinje, using a generous overestimate
	if !isempty(purk)
		iran,jran = [ extrema(i[k] for i in purk) for k in 1:2 ]
		h,w = iran[end]-iran[1],jran[end]-jran[1]
		cen = round.(Int,(median(iran),median(jran)))
		rows = clamp.(cen[1]-2h:cen[1]+2h,1,m)
		cols = clamp.(cen[2]-2w:cen[2]+2w,1,n)
		G[rows,cols] .= median(G[rows,cols])
		end

	# smooth and interpolate green channel
	k = m/options.blur_width
	ker = KernelFactors.gaussian((k,k))
	Z = interpimage(imfilter(G,ker))

	return Z,θ,u_init
end

function detect(sz,Z,θ,u_init;options=get_defaults())
	# Given data and multiple intitializations and options, try them all and find the best.
	u,fmin,best = [],Inf,[]
	for ui in u_init, method in options.optimizers
		unew,fnew = fitcircle(Z,sz...,ui...,method,θ;options)
		@debug "latest: $unew,$fnew,$ui,$method"
		if fnew < fmin
			u,fmin = unew,fnew
			best = (ui,method)
		end
	end
	return u,fmin,best
end

function detect(img::AbstractMatrix{T} where T <: Colorant,sz=size(img);options=get_defaults())
	detect(sz,detectiondata(img,sz...)...;options)
end

"""
	detect(files[,size])
	detect(img[,size])
Perform cornea detection on all the image files named in the vector `files`, or on the image `img`. If `size` is given, all images to be resized to match it.

Returns a named tuple with the following fields:
* `cenrow` row location of the cornea center, normalized by the image height
* `cencol` column location of the cornea center, normalized by the image height
* `radius` radius of the cornea center, normalized by the image height
* `fmin` optimized value of the optimization objective function
* `init` initial condition that led to the optimal result
* `method` info about the optimizer that got the optimal result
"""
function detect(folder::AbstractVector{String},sz=[];options=get_defaults())
	result = (cenrow = Float32[], cencol = Float32[], radius = Float32[], fmin = Float32[], init = [], method = [] )
	if isempty(sz)
		sz = size(load(folder[1]))
	end
	@showprogress for fname in folder
		img = load(fname)
		Z,θ,u_init = detectiondata(img,sz...;options)
		if !isempty(result.radius)
			new_ui = sz[1].*(result.cenrow[end],result.cencol[end],result.radius[end])
			push!(u_init,new_ui)
		end
		u,fmin,best = detect(sz,Z,θ,u_init;options)

		push!(result.cenrow,u[1]/sz[1])
		push!(result.cencol,u[2]/sz[1])
		push!(result.radius,u[3]/sz[1])
		push!(result.fmin,fmin)
		push!(result.init,best[1])
		push!(result.method,best[2])
	end
	return result
end

max_stars_repo_name: tobydriscoll/cornea-detection

max_issues_repo_name: tobydriscoll/cornea-detection

max_forks_repo_name: tobydriscoll/cornea-detection


text:
const JSON_RESP_HDRS = Dict{String,String}("Content-Type" => "application/json; charset=utf-8")
const configuration = Dict(
    :virtualhost => "/",
    :host => "localhost",
    :port => AMQPClient.AMQP_DEFAULT_PORT,
    :auth_params => AMQPClient.DEFAULT_AUTH_PARAMS,
    :inactivity_timeout_secs => 240,
    :debug_level => Logging.WARNING
)

function configure(d::Dict)
    for (n,v) in d
        configuration[Symbol(n)] = v
    end
    api_clients.timeout = configuration[:inactivity_timeout_secs]
    Logging.configure(; level=configuration[:debug_level])
    nothing
end

amqp_params() = filter((n,v)->(n in (:virtualhost, :host, :port, :auth_params)), configuration)

"""
A message queue based client.
Can call multiple services.
Caches connections.
"""
type APIClient
    srvrs::Dict{String,Tuple}
    timeout::Int
    purged_time::Float64

    function APIClient(timeout::Int=typemax(Int))
        new(Dict{String,Tuple}(), timeout, time())
    end
end

"""Singleton client instance used by exported APIs"""
const api_clients = APIClient(configuration[:inactivity_timeout_secs])

function srvrcall(srvrname::String, command::String, args...; kwargs...)
    invoker,transport = getclient(srvrname)

    info("invoking ", srvrname, ": ", command)
    result = apicall(invoker, command, args...; kwargs...)
    info("got result ", result)

    if "data" in keys(result)
        data = result["data"]
        if "code" in keys(data)
            code = data["code"]
            if code == 0
                return data["data"]
            end
        end
    end
    error("Server error: ", result)
end

function getclient(name::String)
    clnt = api_clients
    if hasclient(name)
        invoker,transport,lastused = clnt.srvrs[name]
    else
        info("creating a conenction to ", name, "...")
        transport = JuliaWebAPIPlugins.AMQPTransport(name, :client; amqp_params()...)
        msgformat = JuliaWebAPI.SerializedMsgFormat()
        invoker = APIInvoker(transport, msgformat)
    end
    clnt.srvrs[name] = (invoker,transport,time())
    invoker,transport
end

function hasclient(name::String)
    clnt = api_clients
    !(name in keys(clnt.srvrs)) && (return false)
    invoker,transport,lastused = clnt.srvrs[name]
    if (time() - lastused) > clnt.timeout
        closeclient(name)
        return false
    end
    if (time() - clnt.purged_time) > (clnt.timeout * 5)
        purgeclients()
    end
    true
end

function purgeclients()
    clnt = api_clients
    for name in collect(keys(clnt.srvrs))
        invoker,transport,lastused = clnt.srvrs[name]
        if (time() - lastused) > clnt.timeout
            closeclient(name)
        end
    end
    clnt.purged_time = time()
    nothing
end

function closeclient(name::String)
    clnt = api_clients
    if name in keys(clnt.srvrs)
        invoker,transport,lastused = clnt.srvrs[name]
        close(transport)
        delete!(clnt.srvrs, name)
    end
end

"""
A message queue based service.
Receives a request to process from an upstream server.
Exits on inactivity, or upon receiving an exit command from upstream.
"""
type APIServer
    name::String
    srvr::APIResponder
    srvr_task::Union{Void,Task}
    last_busy::Float64
    command::Function
    atexit::Function
    inactivity_timeout::Int

    function APIServer(command, atexit, name::String)
        transport = JuliaWebAPIPlugins.AMQPTransport(name, :server; amqp_params()...)
        info("APIServer ", name, " created transport")
        msgformat = JuliaWebAPI.SerializedMsgFormat()

        srvr = APIResponder(transport, msgformat)
        endpt = JuliaWebAPI.default_endpoint(command)
        register(srvr, command_wrapper; endpt=endpt, resp_json=true, resp_headers=JSON_RESP_HDRS)
        register(srvr, atexit; endpt="exit", resp_json=true, resp_headers=JSON_RESP_HDRS)

        new(name, srvr, nothing, time(), command, atexit, configuration[:inactivity_timeout_secs])
    end
end

function run(srvr::APIServer)
    global srvr_instance = srvr
    info("APIServer $(srvr.name) starting...")
    srvr.srvr_task = @async process(srvr.srvr)

    # wait till we are receiving requests
    while !istaskdone(srvr.srvr_task) && ((time() - srvr.last_busy) < srvr.inactivity_timeout)
        sleep(5)
    end
    reason = istaskdone(srvr.srvr_task) ? "upstream command" : "inactivity"
    info("APIServer $(srvr.name) exiting due to $reason")
    close(srvr.srvr.transport)
    srvr.atexit()
    nothing
end

function command_wrapper(args...; kwargs...)
    global srvr_instance
    done = false
    keepalive_task = @async begin
        while !done
            srvr_instance.last_busy = time()
            sleep(30)
        end
    end

    try
        srvr_instance.command(args...; kwargs...)
    finally
        done = true
    end
end

max_stars_repo_name: tanmaykm/JuliaWebAPIPlugins

max_issues_repo_name: tanmaykm/JuliaWebAPIPlugins

max_forks_repo_name: tanmaykm/JuliaWebAPIPlugins


text:
function pow(a,b)
    k = b
    t = 1
    p = a
    while k > 0
        if k%2 == 0
            k ÷= 2
            p *= p
        else
            k -= 1
            t *= p
        end
    end
    return t
end


text:
using GaussianInterpolation, CMBLensing
using Test

@test size(coords(zeros(3,2))) == (6,2)
@test size(coords(FlatMap(rand(3,3)))) == (9, 2)


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_m5_opt_rho1_68(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m5_opt_rho1_68!(A_copy)
end

@inline function exp_sastre_m5_opt_rho1_68!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=-0.2940691417404044
    coeff2=0.48301970211963724
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-0.2940691417404044
    coeff2=0.48301970211963724
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.3371652952117715
    coeff2=0.4456299027272144
    coeff3=0.08435039764895197
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.04182430699533409
    coeff2=-0.19033210762699687
    coeff3=0.9365591932640563
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.0007801504959772192
    coeff2=0.03693404820491456
    coeff3=-0.3358998734592709
    coeff4=0.9668380398527895
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.008399213058468373
    coeff2=0.0017935941807159792
    coeff3=0.0005218334570420487
    coeff4=7.480755260551737e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.11409149134407987
    coeff2=0.6872926560470709
    coeff3=0.37829567639306994
    coeff4=0.03390388537056082
    coeff5=0.36476063896687405
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.13828449809707125
    coeff2=0.13000362502266968
    coeff3=0.0628565346339682
    coeff4=0.018703010057796603
    coeff5=0.6395504919933279
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.968291126022034
    coeff2=0.5165851751813066
    coeff3=0.4739135345011139
    coeff4=0.11560763975481164
    coeff5=5.160367463472579
    coeff6=0.6818767988800976
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.968291126022034
    coeff2=0.5165851751813066
    coeff3=0.4739135345011139
    coeff4=0.11560763975481164
    coeff5=5.160367463472579
    coeff6=0.6818767988800976
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.055527011448075834
    coeff2=-0.01959281862440215
    coeff3=-0.01516878804519402
    coeff4=-0.002839711497135245
    coeff5=0.01435796676858081
    coeff6=0.009741431366483796
    coeff7=1.02397276682243
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
points(H::Hermite,n)=gausshermite(n)[1]

plan_transform(H::Hermite,v::Vector) = gausshermite(length(v))
plan_itransform(H::Hermite,cfs::Vector) = points(H,length(cfs))
function transform(H::Hermite,vals,plan::Tuple{Vector,Vector})
    x,w = plan
    V=hermitep(0:length(vals)-1,x)'
    nrm=(V.^2)*w

    V*(w.*vals)./nrm
end
itransform(H::Hermite,cfs,plan::Vector) = hermitep(0:length(cfs)-1,tocanonical(H,plan))*cfs


text:
using Test, ModelConstructors
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

# TEST USING SLM OBJECT THAT IT MATCHES THE DESIRED OUTPUT
rp = joinpath(dirname(@__FILE__), "../reference/SLM")
in_inc = load(joinpath(rp, "eval.jld2"))
in_dec = load(joinpath(rp, "eval_decrease.jld2"))
in_sin = load(joinpath(rp, "eval_sine.jld2"))
in_sinnoineq = load(joinpath(rp, "eval_sinnoineq.jld2"))

slms = Vector{AbstractSLM}(undef, 0)

for (i, in_data) in zip(1:4, [in_inc, in_dec, in_sin, in_sinnoineq])
    if i == 1
        push!(slms, SLM(vec(in_data["w"]), vec(in_data["p_sol"]); knots = Int(in_data["knots"]),
                        increasing = true, concave_down = true, left_value = in_data["p0_norun"],
                        right_value = in_data["p1_norun"], min_value = in_data["p0_norun"] - 1e-3,
                        max_value = in_data["p1_norun"] + 1e-3))
    elseif i == 2
        push!(slms, SLM(vec(in_data["w"]), vec(in_data["rev_p_sol"]); knots = Int(in_data["knots"]),
                        decreasing = true, concave_up = true, right_value = in_data["p0_norun"],
                        left_value = in_data["p1_norun"], min_value = in_data["p0_norun"] - 1e-3,
                        max_value = in_data["p1_norun"] + 1e-3))
    elseif i == 3
        push!(slms, SLM(vec(in_data["x"]), vec(in_data["y"]); knots = Int(in_data["knots"]),
                        increasing_intervals = in_data["inc_int"], decreasing_intervals = in_data["dec_int"],
                        concave_up_intervals = in_data["cu_int"], concave_down_intervals = in_data["cd_int"],
                        right_value = in_data["rightvalue"], left_value = in_data["leftvalue"]))
    else
        push!(slms, SLM(vec(in_data["x"]), vec(in_data["y"]); knots = Int(in_data["knots"]),
                        right_value = in_data["rightvalue"], left_value = in_data["leftvalue"]))
    end
end

@testset "eval of SLM" begin
    for (i, in_data) in zip(1:4, [in_inc, in_dec, in_sin, in_sinnoineq])
        if i == 1
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["w"])) in_data["phat"]
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["w"]), 0) in_data["phat"]
        elseif i == 2
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["w"])) in_data["phat"]
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["w"]), 0) in_data["phat"]
        elseif i == 3
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["xhat"])) in_data["yhat"]
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["xhat"]), 0) in_data["yhat"]
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["xhat"]), 1) in_data["yphat"]
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["xhat"]), 2) in_data["ypphat"]
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["xhat"]), 3) in_data["yppphat"]
        else
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["xhat"])) in_data["yhat"]
            @test @test_matrix_approx_eq eval(slms[i], vec(in_data["xhat"]), 0) in_data["yhat"]
        end
    end
end


text:
# Recreate plot with CMF vs. orbital period:

include("../../../src/regress.jl")
using Statistics
using PyPlot
using DelimitedFiles
using Printf

period = [1.510826, 2.421937,4.049219, 6.101013, 9.207540,12.352446,18.772866]

dir = "../../../data/POSTERIOR_CMF/"
fname = ["T1b_IRON_FRAC.out","T1c_IRON_FRAC.out","T1d_IRON_FRAC.out","T1e_IRON_FRAC.out","T1f_IRON_FRAC.out","T1g_IRON_FRAC.out","T1h_IRON_FRAC.out"]

np = 7; nsamp = 10000

feonmg(cmf) = 1.778/(100/cmf-1)

cmf_norm = zeros(np,nsamp)
cmf_mean = zeros(np)
cmf_median = zeros(np)
cmf_sig = zeros(np)
cmf_sig1 = zeros(np)
cmf_sig2 = zeros(np)
feonmg_median = zeros(np)
feonmg_sig1 = zeros(np)
feonmg_sig2 = zeros(np)

for ip=1:np
  data = readdlm(string(dir,fname[ip]))
  cmf_norm[ip,:] .= data[:,2]
  cmf_mean[ip] = mean(cmf_norm[ip,:])
  cmf_median[ip] = median(cmf_norm[ip,:])
  cmf_sig[ip] = std(cmf_norm[ip,:])
  cmf_sig1[ip] = cmf_median[ip]-sort(cmf_norm[ip,:])[1587]
  cmf_sig2[ip] = sort(cmf_norm[ip,:])[8413]-cmf_median[ip]
  feonmg_median[ip] = median(feonmg.(cmf_norm[ip,:]))
  feonmg_sig1[ip] = feonmg_median[ip]-sort(feonmg.(cmf_norm[ip,:]))[1587]
  feonmg_sig2[ip] = sort(feonmg.(cmf_norm[ip,:]))[8413]-feonmg_median[ip]
end

# Now, carry out regressions for each of these:
coeff_samp = zeros(2,nsamp)
fn = zeros(2,np)
fn[1,:] .= 1.0
fn[2,:] .= period ./6
sig = ones(np)
cmf_avg = 0.0
for isamp=1:nsamp
  global cmf_avg += mean(cmf_norm[:,isamp])
  coeff,cov = regress(fn,cmf_norm[:,isamp],sig)
  coeff_samp[:,isamp] = coeff
end

cmf_avg /= nsamp

# Now, compute standard deviation:
cmf_var = 0.0
for isamp=1:nsamp
  global cmf_var += (mean(cmf_norm[:,isamp]) - cmf_avg)^2
end
cmf_rms = sqrt(cmf_var/nsamp)

# Now, draw 1-sigma and 2-sigma confidence intervals:
period_trial = collect(0.0:0.1:20.0)
conf = zeros(5,length(period_trial))
for (i,p) in enumerate(period_trial)
  # Predict samples at this period
  csamp = coeff_samp[1,:] .+ coeff_samp[2,:] .* p /6
  # Now, sort these and pick out 68.3% & 95% confidence limits:
  csamp = sort(csamp)
  conf[:,i] = [csamp[250],csamp[1415],csamp[5000],csamp[8415],csamp[9750]]
end
# Now, plot median, 1-sigma and 2-sigma confidence intervals:

clf()
#plot(period_trial,conf[1,:],color="C1",alpha=0.5,linestyle="--")
errorbar(period,cmf_mean,yerr=cmf_sig,fmt="o")
#errorbar(period,cmf_median,yerr=[cmf_sig1';cmf_sig2'],fmt="o")
#plot(period_trial,conf[5,:],color="C1",alpha=0.5,linestyle="--")
plot(period_trial,conf[3,:],color="C1",linewidth=2)
plot(period_trial,conf[2,:],color="C1",linewidth=2,linestyle=":")
plot(period_trial,conf[4,:],color="C1",linewidth=2,linestyle=":")
#plot(period_trial,cmf_avg  .+ 0.0 .* period_trial,color="C0",linewidth=2)
fill_between([0,20],[cmf_avg-cmf_rms,cmf_avg-cmf_rms],[cmf_avg+cmf_rms,cmf_avg+cmf_rms],alpha=0.2,color="C0")
xlabel("Orbital Period [d]")
ylabel("CMF/Iron mass fraction (%)")
# Overplot solar system bodies:
text([16],[28],"Mars")
plot([0,20],[30.0,30.0],"--",color="k")
text([1],[8],"Moon")
plot([0,20],[10.0,10.0],"--",color="k")
text([18],[33],"Earth")
plot([0,20],[32.0,32.0],"--",color="k")
#text([17],[31.5],"Venus")
#plot([0,20],[31.,31.],"--",color="k")
axis([0,20,0,40])
text([15],[20],"Chondrites (L/CI)")
#plot([0,20],[22.0,22.0],"--",color="k")
plot([0,20],[18.7,18.7],"--",color="k")
text([1],[20.6],"b")
text([2.7],[22.6],"c")
text([3.5],[13],"d")
text([6.3],[20],"e")
text([8.7],[15],"f")
text([11.8],[12],"g")
text([18.1],[11.8],"h")

tight_layout()
savefig("../Norm_cmf_vs_period.pdf",bbox_inches = "tight")

println("Coefficients of fit: ",mean(coeff_samp[1,:]),"+-",std(coeff_samp[1,:])," ",mean(coeff_samp[2,:]),"+-",std(coeff_samp[2,:]))

# Now, print out first line of Table 9:
tab09_01 = "CMF [wt\\%] & "
for i=1:7
  global tab09_01 = string(tab09_01," \$",@sprintf("%4.1f",cmf_median[i]),"_{-",@sprintf("%4.1f",cmf_sig1[i]),"}^{+",@sprintf("%4.1f",cmf_sig2[i]),"}\$ &")
end

tab09_01 = string(tab09_01," \$",@sprintf("%4.1f",cmf_avg),"\\pm",@sprintf("%4.1f",cmf_rms),"\$ \\cr")
println(tab09_01)

# Now, print out second line of Table 9:
tab09_02 = "Fe/Mg molar ratio & "
for i=1:7
  global tab09_02 = string(tab09_02," \$",@sprintf("%4.2f",feonmg_median[i]),"_{-",@sprintf("%4.2f",feonmg_sig1[i]),"}^{+",@sprintf("%4.2f",feonmg_sig2[i]),"}\$ &")
end

tab09_02 = string(tab09_02," \$",@sprintf("%4.2f",feonmg.(cmf_avg)),"\\pm",@sprintf("%4.2f",feonmg.(cmf_rms)),"\$ \\cr")
println(tab09_02)


text:
export BasisFunction,
    index

abstract type TypedFunction{S,T} end

domaintype(::Type{<:TypedFunction{S,T}}) where {S,T} = S
domaintype(φ::TypedFunction) = domaintype(typeof(φ))

codomaintype(::Type{<:TypedFunction{S,T}}) where {S,T} = T
codomaintype(φ::TypedFunction) = codomaintype(typeof(φ))

prectype(F::Type{<:TypedFunction}) = prectype(domaintype(F),codomaintype(F))

eltype(::Type{<:Dictionary{S,T}}) where {S,T} = TypedFunction{S,T}

promote_rule(::Type{<:TypedFunction{S1,T1}}, ::Type{<:TypedFunction{S2,T2}}) where {S1,T1,S2,T2} =
    TypedFunction{promote_type(S1,S2),promote_type(T1,T2)}

"The supertype of functions that can be associated with a dictionary or a family of basis functions."
abstract type AbstractBasisFunction{S,T} <: TypedFunction{S,T} end

"A `BasisFunction` is one element of a dictionary."
struct BasisFunction{S,T,D<:Dictionary{S,T},I} <: AbstractBasisFunction{S,T}
    dictionary  ::  D
    index       ::  I
end

dictionary(φ::BasisFunction) = φ.dictionary
index(φ::BasisFunction) = φ.index

support(φ::BasisFunction) = support(dictionary(φ), index(φ))
measure(φ::BasisFunction) = measure(dictionary(φ))

show(io::IO, mime::MIME"text/plain", φ::BasisFunction) = composite_show(io, mime, φ)
Display.stencil_parentheses(φ::BasisFunction) = true
Display.displaystencil(φ::BasisFunction) = [dictionary(φ), '[', index(φ), ']']

(φ::BasisFunction)(x) = unsafe_eval_element1(dictionary(φ), index(φ), x)
(φ::BasisFunction)(x, y...) = φ(SVector(x, y...))

basisfunction(dict::Dictionary, idx) = BasisFunction(dict, idx)

function getindex(dict::Dictionary, idx)
    @boundscheck checkbounds(dict, idx)
    basisfunction(dict, idx)
end
function getindex(dict::Dictionary, i, j, indices...)
    @boundscheck checkbounds(dict, i, j, indices...)
    basisfunction(dict, (i,j,indices...))
end

# Inner product with a basis function: we choose the measure associated with the dictionary
innerproduct(φ::AbstractBasisFunction, g; options...) =
    innerproduct(φ, g, measure(φ); options...)
innerproduct(f, ψ::AbstractBasisFunction; options...) =
    innerproduct(f, ψ, measure(ψ); options...)
innerproduct(ϕ::AbstractBasisFunction, ψ::AbstractBasisFunction; options...) =
    innerproduct(ϕ, ψ, defaultmeasure(dictionary(ϕ), dictionary(ψ)); options...)

# The inner product between two basis functions: invoke the implementation of the dictionary
innerproduct(φ::AbstractBasisFunction, ψ::AbstractBasisFunction, measure; options...) =
    innerproduct(dictionary(φ), index(φ), dictionary(ψ), index(ψ), measure; options...)

# The inner product of a basis function with another function: this is an analysis integral
# We introduce a separate function name for this for easier dispatch.
innerproduct(φ::AbstractBasisFunction, g, measure; options...) =
    analysis_integral(dictionary(φ), index(φ), g, measure; options...)
innerproduct(f, ψ::AbstractBasisFunction, measure; options...) =
    conj(analysis_integral(dictionary(ψ), index(ψ), f, measure; options...))

# We want to check whether the supports of the basis function and the measure differ.
# The integral may be easier to evaluate by restricting to the intersection of these
# supports. However, we only perform this optimization if the intersection does not
# lead to an IntersectDomain (i.e., the intersection is not recognized) since
# that leads to incomputable integrals.
function analysis_integral(dict::Dictionary, idx, g, measure::Measure; options...)
    @boundscheck checkbounds(dict, idx)
    domain1 = support(dict, idx)
    domain2 = support(measure)
    qs = quadstrategy(promote_type(prectype(dict), prectype(measure)); options...)
    unsafe_analysis_integral1(dict, idx, g, measure, domain1, domain2, domain1 ∩ domain2, qs)
end

function analysis_integral(dict::Dictionary, idx, g, measure::DiscreteWeight; options...)
    @boundscheck checkbounds(dict, idx)
    unsafe_analysis_integral2(dict, idx, g, measure, support(dict, idx))
end

# unsafe for indexing
function unsafe_analysis_integral1(dict, idx, g, measure::Weight, d1, d2, domain, qs)
    if d1 == d2
        integral(qs, x->conj(unsafe_eval_element(dict, idx, x))*g(x), measure)
    else
        # -> do compute on the smaller domain and convert the measure
        integral(qs, x->conj(unsafe_eval_element(dict, idx, x))*g(x), domain, measure)
    end
end

unsafe_analysis_integral1(dict, idx, g, measure, d1, d2, domain::IntersectDomain, qs) =
    # -> disregard the intersection domain, but use eval_element to guarantee correctness
    integral(qs, x->conj(eval_element(dict, idx, x))*g(x), measure)

# unsafe for indexing and for support of integral
unsafe_analysis_integral2(dict::Dictionary, idx, g, measure::DiscreteWeight, domain) =
    integral(x->conj(unsafe_eval_element(dict, idx, x))*g(x), domain, measure)


text:
"""
    NaiveDecoder <: Decoder

    NaiveDecoder([max_qubits=10])

Construct a naive decoder that iterates through all possible errors, in ascending weight,
and resolves to the first error that matches the syndrome.

!!! note

    This decoder is slow for even moderate numbers of qubits. By default, it is restricted
    to codes with a maximum of 10 qubits.
"""
struct NaiveDecoder <: Decoder
    max_qubits::Int
    NaiveDecoder(max_qubits=10) = new(max_qubits)
end
Model.label(::NaiveDecoder) = "Naive"
function Model.decode(
    decoder::NaiveDecoder, code, syndrome::AbstractVector{Bool}; kwargs...
)
    if 0 <= decoder.max_qubits < nkd(code)[1]
        throw(ArgumentError("NaiveDecoder restricted to $(decoder.max_qubits) qubits"))
    end
    recovery = _minimum_weight_recovery(code, syndrome)
    return DecodeResult(recovery=recovery)
end

function _minimum_weight_recovery(code, syndrome)
    n_qubits = nkd(code)[1]
    min_weight = 0
    max_weight = n_qubits
    # iterate recovery operations of increasing weight
    for weight in min_weight:max_weight
        for qubit_indices in combinations(1:n_qubits, weight)
            for xyzs in Iterators.product(fill(1:3, weight)...)
                # pauli encoded as 0123 for IXZY
                v = zeros(Int, n_qubits)
                v[qubit_indices] = collect(xyzs)
                # pauli encoded in binary symplectic form
                recovery = BitVector(vcat(v .% 2, v .÷ 2))
                if bsp(stabilizers(code), recovery) == syndrome
                    return recovery
                end
            end
        end
    end
    return falses(n_qubits * 2)
end


text:
using Flux: Chain, train!, params
using Flux.Losses:mse
using Flux.Data:DataLoader
using Flux.Optimise:ADAM
using IterTools:ncycle
using Statistics:mean

"""
    AEDetector(encoder= Chain(),
               decoder = Chain(),
               batchsize= 32,
               epochs = 1,
               shuffle = false,
               partial = true,
               opt = ADAM(),
               loss = mse)

Calculate the anomaly score of an instance based on the reconstruction loss of an autoencoder, see [1] for an
explanation of auto encoders.

Parameters
----------
$AE_PARAMS

$AE_LOSS

Examples
--------
$(SCORE_UNSUPERVISED("AEDetector"))

References
----------
[1] Aggarwal, Charu C. (2017): Outlier Analysis.
"""
OD.@detector struct AEDetector <: UnsupervisedDetector
    encoder::Chain = Chain()
    decoder::Chain = Chain()
    batchsize::Integer = 32
    epochs::Integer = 1
    shuffle::Bool = false
    partial::Bool = true
    opt::Any = ADAM()
    loss::Function = mse
end

struct AEModel <: DetectorModel
    chain::Chain
end

function OD.fit(detector::AEDetector, X::Data; verbosity)::Fit
    loader = DataLoader(X, batchsize=detector.batchsize, shuffle=detector.shuffle, partial=detector.partial)

    # Create the autoencoder
    model = Chain(detector.encoder, detector.decoder)

    # train the neural network model
    train!(x -> detector.loss(model(x), x), params(model), ncycle(loader, detector.epochs), detector.opt)

    scores = detector.loss(model(X), X, agg=instance_mean)
    return AEModel(model), scores
end

function OD.transform(detector::AEDetector, model::AEModel, X::Data)::Scores
    detector.loss(model.chain(X), X, agg=instance_mean)
end

max_stars_repo_name: davnn/OutlierDetectionNetworks

max_issues_repo_name: davnn/OutlierDetectionNetworks

max_forks_repo_name: davnn/OutlierDetectionNetworks


text:
module q5

import Evolve

steps = 8000                # number of simulation steps
samplesteps = 50            # steps between measurements
L = 2                       # width of the box
xmin = -1                   # leftmost point in the box
N = 2000                    # number of subintervals
ϵ = L/N                     # width of each subinterval
h = 1e5                     # energy is 12336.8 for w = 0.02
b = 1
s = 0.005                    # wall thickness
w = 0.02         # pocket widths to try
ks = π/w .* [2.4:0.2:3.2]                  # wavenumbers
σ = 0.05                    # wavepacket width
V0s = [0]         # step potentials
x0 = -0.5                   # wavepacket starting position
Ts = Float64[]
Rs = Float64[]

# Run the trials
for k in ks
    for V0 in V0s
        print("Initializing for k = $k... ")
        V(xs) = map(x -> Evolve.pocket(L, h, b, s, w, xmin, x), xs)      # pocket potential
        AV = Evolve.A(V)            # evolution matrix constructor for square well
        ϕV = Evolve.Φ(k, σ, x0)     # initial state vector constructor for Gaussian packet

        A = AV(L, N, xmin)
        ϕ0 = ϕV(L, N, xmin)
        x = (L/N .* [0:N] + xmin)   # x positions

        # Evolve the system and record positions every so many steps
        print("evolving... ")
        test = x -> (abs(x[2]) > 1e-2 || abs(x[end-1]) > 1e-2)
        ϕ = Evolve.evolve(A, ϕ0, steps, test)
        phis = [ϕ0 ϕ[:,samplesteps:samplesteps:end]]

        # calculate transmission and reflection probabilities
        T = mean(abs(ϕ[N/2+1:end,end]).^2)
        R = mean(abs(ϕ[1:N/2,end]).^2)
        I = T + R
        T /= I
        R /= I
        push!(Ts, T)
        push!(Rs, R)
        print("T = $T, R = $R...")
        ### writedlm("trans_step_$V0", [T,R,I])

        println("done after $(size(ϕ)[2]) steps. Writing to file.")
        writedlm("phis_scatter_im_$k.dat", imag(phis))
        writedlm("phis_scatter_re_$k.dat", real(phis))
    end
end

writedlm("scatter1.dat", [ks Ts Rs])
end


text:
const NaNCF64 = (NaN+NaN*im) :: ComplexF64

function findconnected(xs::Vector, cmp=isequal)
    # NOTE: `cmp` must be symmetric.
    components = Vector{Int}[]
    visited = fill(false, length(xs))
    for i = 1:length(xs)
        visited[i] && continue # the value is already part of a component
        component = Int[]
        _find_component!(component, xs, visited, cmp, i)
        push!(components, component)
    end
    return components
end

function _find_component!(component, xs, visited, cmp, i)
    push!(component, i)
    visited[i] = true
    for j = (i+1):length(xs)
        visited[j] && continue # the value is already part of a component
        cmp(xs[i], xs[j]) && _find_component!(component, xs, visited, cmp, j)
    end
end

function findsimilar(xs::Vector{T}; atol=sqrt(eps(T))) where T <: Number
    ccs = findconnected(xs, (x,y) -> isapprox(x, y; atol=atol))
    map(ccs) do idxs
        idxs => mean(xs[idxs])
    end
end

countunique(xs) = [uq => count(isequal(uq), xs) for uq in unique(xs)]
findunique(xs) = [uq => findall(isequal(uq), xs) for uq in unique(xs)]

function subspace_matrix(M::AbstractMatrix, vectors::AbstractMatrix, subspace::Vector{Int})
    vs = vectors[:, subspace]
    # NOTE: also computing the transpose seems to give some numeric stability
    (vs' * M * vs + (vs' * M' * vs)')/2
end

function simeigen(A::AbstractMatrix, Bs...; atol=1e-10)
    for B in Bs
        @assert typeof(B) <: AbstractMatrix
        @assert size(A) == size(B)
    end

    e = eigen(A)
    if length(Bs) == 0
        return (values = e.values, vectors = e.vectors)
    else
        values = fill(NaNCF64, (length(Bs) + 1, length(e.values)))
        vectors = fill(NaNCF64, size(e.vectors))
        subspaces = findsimilar(e.values; atol=atol)
        for (ids, ev) in subspaces
            Bs_ss = [subspace_matrix(B, e.vectors, ids) for B in Bs]
            simeig = simeigen(Bs_ss...; atol=atol)
            values[1, ids] .= ev
            values[2:end, ids] = simeig.values
            vectors[:,ids] = e.vectors[:,ids] * simeig.vectors
        end
        return (values = values, vectors = vectors)
    end
end


text:
using ComplexSIMD
using Base.Test
using StaticArrays
using SIMD
using BenchmarkTools

a = @SVector rand(Complex128, 4)
b = @SVector rand(Complex128, 4)

d0, d1 = ComplexSIMD.cmul(a, b)

@btime ComplexSIMD.cmul($a, $b)

function naive(x, y)
    x .* y
end

@btime naive($a, $b)

@code_native ComplexSIMD.cmul(a, b)

# write your own tests here
@test 1 == 2


text:
# Autogenerated wrapper script for Gama_jll for aarch64-linux-gnu-cxx11
export gama_g3, gama_local, gama_local_gkf2yaml, gama_local_xml2sql, gama_local_xml2txt, gama_local_yaml2gkf

using Expat_jll
using yaml_cpp_jll
using SQLite_jll
JLLWrappers.@generate_wrapper_header("Gama")
JLLWrappers.@declare_executable_product(gama_g3)
JLLWrappers.@declare_executable_product(gama_local)
JLLWrappers.@declare_executable_product(gama_local_gkf2yaml)
JLLWrappers.@declare_executable_product(gama_local_xml2sql)
JLLWrappers.@declare_executable_product(gama_local_xml2txt)
JLLWrappers.@declare_executable_product(gama_local_yaml2gkf)
function __init__()
    JLLWrappers.@generate_init_header(Expat_jll, yaml_cpp_jll, SQLite_jll)
    JLLWrappers.@init_executable_product(
        gama_g3,
        "bin/gama-g3",
    )

    JLLWrappers.@init_executable_product(
        gama_local,
        "bin/gama-local",
    )

    JLLWrappers.@init_executable_product(
        gama_local_gkf2yaml,
        "bin/gama-local-gkf2yaml",
    )

    JLLWrappers.@init_executable_product(
        gama_local_xml2sql,
        "bin/gama-local-xml2sql",
    )

    JLLWrappers.@init_executable_product(
        gama_local_xml2txt,
        "bin/gama-local-xml2txt",
    )

    JLLWrappers.@init_executable_product(
        gama_local_yaml2gkf,
        "bin/gama-local-yaml2gkf",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ---
# title: 322. Coin Change
# id: problem322
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/coin-change/description/>
# hidden: true
# ---
#
# You are given coins of different denominations and a total amount of money
# _amount_. Write a function to compute the fewest number of coins that you need
# to make up that amount. If that amount of money cannot be made up by any
# combination of the coins, return `-1`.
#
# You may assume that you have an infinite number of each kind of coin.
#
#
#
# **Example 1:**
#
#
#
#     Input: coins = [1,2,5], amount = 11
#     Output: 3
#     Explanation: 11 = 5 + 5 + 1
#
#
# **Example 2:**
#
#
#
#     Input: coins = [2], amount = 3
#     Output: -1
#
#
# **Example 3:**
#
#
#
#     Input: coins = [1], amount = 0
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: coins = [1], amount = 1
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: coins = [1], amount = 2
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= coins.length <= 12`
#   * `1 <= coins[i] <= 231 - 1`
#   * `0 <= amount <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2b78a920-070d-11eb-194e-bd44cc52fd06
using Pkg, DrWatson

# ╔═╡ 0fc3846a-070e-11eb-1144-4918bcd157f3
begin
	@quickactivate "SR2TuringPluto"
	using Turing
	using StatisticalRethinking
end;

# ╔═╡ dfe57196-070c-11eb-1a81-fdf2e39f38ae
md"## Clip-04-48-49t.jl"

# ╔═╡ 2d8739d4-070d-11eb-1bd2-c9a5fd864e4d
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	df = df[df.age .>= 18, :]
	x̄ = mean(df.weight)
	x = range(minimum(df.weight), maximum(df.weight), length = 100)     # either
end;

# ╔═╡ 37566494-070d-11eb-0c3b-694670d5f2ea
@model function m4_3(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* (weights .- x̄)
    heights ~ MvNormal(μ, σ)
end

# ╔═╡ d987cc24-070c-11eb-12ba-4b04b7f9c218
md"### snippets 4.48, 4.49"

# ╔═╡ 04032cd8-070d-11eb-3610-0b5e5c16c8ec
begin
	figs = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3)
		for (indx, N) in enumerate([50, 150, 352])

		dN = df[1:N, :]
		quapN = quap(m4_3(dN.weight, dN.height))

		post = rand(quapN.distr, 20)
		post = DataFrame(post', quapN.params)

		figs[indx] = scatter(dN.weight, dN.height)
		for p in eachrow(post)
			plot!(x, p.a .+ p.b .* (x .- mean(dN.weight)), color = "black", alpha = 0.3)
		end
		plot!(legend = false, xlabel = "weight", ylabel = "height")
	end
	plot(figs..., layout=(1, 3))
end

# ╔═╡ b8f032bc-070f-11eb-2ad9-c33312592d49
md"## End of clip-04-48-49t.jl"

# ╔═╡ Cell order:
# ╟─dfe57196-070c-11eb-1a81-fdf2e39f38ae
# ╠═2b78a920-070d-11eb-194e-bd44cc52fd06
# ╠═0fc3846a-070e-11eb-1144-4918bcd157f3
# ╠═2d8739d4-070d-11eb-1bd2-c9a5fd864e4d
# ╠═37566494-070d-11eb-0c3b-694670d5f2ea
# ╟─d987cc24-070c-11eb-12ba-4b04b7f9c218
# ╠═04032cd8-070d-11eb-3610-0b5e5c16c8ec
# ╟─b8f032bc-070f-11eb-2ad9-c33312592d49


text:
# ------------------------------------------------------------
# Reconstruction of impulse responses, i.e., a single coefficient in the Haar
# basis from equidistant frequency samples

# Define scaling function and mother wavelet

function Haar_scaling(x::Real)
	if x >= 0 && x <= 1.0
		return 1.0
	else
		return 0.0
	end
end

function Haar_wavelet(x::Real)
	if x >= 0 && x < 0.5
		return 1.0
	elseif x >= 0.5 && x <= 1.0
		return -1.0
	else
		return 0.0
	end
end

# ------------------------------------------------------------
# Reconstruction of scaling function

# Input:
# M: The number of frequency samples
# epsislon: The distance between two consecutive samples
# N: The number of wavelet coefficients recovered

M = 32
epsilon = 1
#J = floor(Int, log2(M))
J = 2

# Sampling points
xi = [-M:M;]*epsilon

# Change of basis matrix
Tx, Ty = freq2Haar( xi, xi, J )
T = kron(Ty, Tx)

# Fourier observations
bx = FourHaarScaling(xi, J, 3)
by = FourHaarScaling(xi, J, 3)
b = bx * by.'

# Scaling function coefficients: 1 in the last (16th) entry, 0 elsewhere
y = pinv(T)*b[:] # Brute force
yy = pinv(Tx)*b*pinv(Ty).' # With Kronecker products


text:
module InteractiveDynamics

using Makie, Observables

const MARKER = Circle(Point2f0(0, 0), Float32(1)) # allows pixel size (zoom independent)
const DEFAULT_BG = RGBf0(1.0, 1.0, 1.0)
using Makie: px

# JULIADYNAMICS_COLORS = to_color.(("#7a60bb", "#202020", "#1ba5aa"))
JULIADYNAMICS_COLORS = to_color.(["#7d53e7", "#202020", "#17c7cd"])
JULIADYNAMICS_CMAP = to_color.(["#2c1633", "#744cd7", "#00ffc1"])
export JULIADYNAMICS_COLORS

include("utils.jl")
using Requires

function __init__()
    @require DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" begin
        # include("numericdata/plot_dataset.jl")
        # include("numericdata/trajectory_highlighter.jl")
        include("chaos/orbitdiagram.jl")
        include("chaos/poincareclick.jl")
        include("chaos/trajanim.jl")
        include("chaos/cobweb.jl")
        include("chaos/brainscan.jl")
    end
    @require DynamicalBilliards = "4986ee89-4ee5-5cef-b6b8-e49ba721d7a5" begin
        include("billiards/defs_plotting.jl")
        include("billiards/defs_animating.jl")
        include("billiards/interactive_billiard.jl")
    end
    @require Agents = "46ada45e-f475-11e8-01d0-f70cc89e6671" begin
        include("agents/stepping.jl")
        include("agents/plots_videos.jl")
        include("agents/interactive_parameters.jl")
    end
end

end