text:
#=##############################################################################
# DESCRIPTION
    Types defining maneuvers of flight vehicles.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

################################################################################
# ABSTRACT MANEUVER TYPE
################################################################################
"""
    `AbstractManeuver{N, M}`

`N` indicates the number of tilting systems in this maneuver, while and `M`
indicates the number of rotor systems.

Every implementation of `AbstractManeuver` must have the properties:

 * `angle::NTuple{N, Function}` where `angle[i](t)` returns the angle of the
        i-th tilting system at time `t` (t is nondimensionalized by the total
        time of the maneuver, from 0 to 1, beginning to end).
 * `RPM::NTuple{M, Function}` where `RPM[i](t)` returns the normalized RPM of
        the i-th rotor system at time `t`. This RPM values are normalized by the
        an arbitrary RPM value (usually RPM in hover or cruise).
"""
abstract type AbstractManeuver{N, M} end

##### FUNCTIONS REQUIRED IN IMPLEMENTATIONS ####################################
"""
    `calc_dV(maneuver::AbstractManeuver, vehicle::Vehicle, t, dt, ttot, Vref)`

Returns the change in velocity `dV=[dVx, dVy, dVz]` (m/s) of `vehicle` performing
`maneuver` at time `t` (s) after a time step `dt` (s).  `Vref` and `tot` are the
reference velocity and the total time at which this maneuver is being performed,
respectively. `dV` is in the global reference system.
"""
function calc_dV(self::AbstractManeuver, vehicle::AbstractVehicle, t::Real,
                                            dt::Real, ttot::Real, Vref::Real)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `calc_dw(maneuver::AbstractManeuver, vehicle::Vehicle, t, dt, Vref, ttot)`

Returns the change in angular velocity `dW=[dWx, dWy, dWz]` (about global
axes, in radians) of `vehicle` performing `maneuver` at time `t` (s) after a
time step `dt` (s). `ttot` is the total time at which this maneuver is to be
performed.
"""
function calc_dW(self::AbstractManeuver, vehicle::AbstractVehicle, t::Real,
                                                        dt::Real, ttot::Real)
    error("$(typeof(self)) has no implementation yet!")
end

##### COMMON FUNCTIONS  ########################################################
"""
    `get_ntltsys(self::AbstractManeuver)`
Return number of tilting systems.
"""
get_ntltsys(self::AbstractManeuver) = typeof(self).parameters[1]

"""
    `get_nrtrsys(self::AbstractManeuver)`
Return number of rotor systems.
"""
get_nrtrsys(self::AbstractManeuver) = typeof(self).parameters[2]

"""
    `get_angle(maneuver::AbstractManeuver, i::Int, t::Real)`

Returns the angle (in degrees) of the i-th tilting system at the non-dimensional
time t.
"""
function get_angle(self::AbstractManeuver, i::Int, t::Real)
    if i<=0 || i>get_ntltsys(self)
        error("Invalid tilting system #$i (max is $(get_ntltsys(self))).")
    end
    if t<0 || t>1
        @warn("Got non-dimensionalized time $(t).")
    end
    return self.angle[i](t)
end

"""
    `get_angles(maneuver::AbstractManeuver, t::Real)`

Returns the angle (in degrees) of every tilting systems at the non-dimensional
time t.
"""
get_angles(self::AbstractManeuver, t::Real) = Tuple(a(t) for a in self.angle)

"""
    `get_RPM(maneuver::AbstractManeuver, i::Int, t::Real)`

Returns the normalized RPM of the i-th rotor system at the non-dimensional time
t.
"""
function get_RPM(self::AbstractManeuver, i::Int, t::Real)
    if i<=0 || i>get_nrtrsys(self)
        error("Invalid rotor system #$i (max is $(get_nrtrsys(self))).")
    end
    if t<0 || t>1
        @warn("Got non-dimensionalized time $(t).")
    end
    return self.RPM[i](t)
end

"""
    `get_RPMs(maneuver::AbstractManeuver, t::Real)`

Returns the normalized RPM of every rotor systems at the non-dimensional time
t.
"""
get_RPMs(self::AbstractManeuver, t::Real) = Tuple(rpm(t) for rpm in self.RPM)

##### COMMON INTERNAL FUNCTIONS  ###############################################

##### END OF ABSTRACT MANEUVER #################################################

################################################################################
# KINEMATIC MANEUVER TYPE
################################################################################
"""
    `KinematicManeuver(angle, RPM, Vvehicle, anglevehicle)`

A vehicle maneuver where the kinematic are prescribed.

# ARGUMENTS
* `angle::NTuple{N, Function}` where `angle[i](t)` returns the angles
        `[Ax, Ay, Az]` (in degrees)of the i-th tilting system at time `t` (t is
        nondimensionalized by the total time of the maneuver, from 0 to 1,
        beginning to end).
* `RPM::NTuple{M, Function}` where `RPM[i](t)` returns the normalized RPM of
        the i-th rotor system at time `t`. This RPM values are normalized by the
        an arbitrary RPM value (usually RPM in hover or cruise).
* `Vvehicle::Function` where `Vvehicle(t)` returns the normalized vehicle
        velocity `[Vx, Vy, Vz]` at the normalized time `t`. Velocity is
        normalized by a reference velocity (typically, cruise velocity).
* `anglevehicle::Function` where `anglevehicle(t)` returns the angles
        `[Ax, Ay, Az]` (in degrees) of the vehicle relative to the global
        coordinate system at the normalized time `t`.
"""
struct KinematicManeuver{N, M} <: AbstractManeuver{N, M}
    angle::NTuple{N, Function}
    RPM::NTuple{M, Function}
    Vvehicle::Function
    anglevehicle::Function
end

# # Implicit N and M constructor
# KinematicManeuver(a::NTuple{N, Function}, b::NTuple{M, Function},
#                     c::Function, d::Function
#                  ) where {N, M} = KinematicManeuver{N, M}(a, b, c, d)

##### FUNCTIONS  ###############################################################
function calc_dV(self::KinematicManeuver, vehicle::AbstractVehicle, t::Real,
                                            dt::Real, ttot::Real, Vref::Real)
    return Vref * (self.Vvehicle((t+dt)/ttot) - self.Vvehicle(t/ttot))
end

function calc_dW(self::KinematicManeuver, vehicle::AbstractVehicle, t::Real,
                                                         dt::Real, ttot::Real)
    prev_W = (self.anglevehicle(t/ttot) - self.anglevehicle((t-dt)/ttot)) / dt
    cur_W = (self.anglevehicle((t+dt)/ttot) - self.anglevehicle(t/ttot)) / dt
    return pi/180 * (cur_W - prev_W)
end

##### INTERNAL FUNCTIONS  ######################################################

##### END OF KINEMATICMANEUVER  ################################################

################################################################################
# KINEMATIC MANEUVER TYPE
################################################################################
struct DynamicManeuver{N, M} <: AbstractManeuver{N, M}
    angle::NTuple{N, Function}
    RPM::NTuple{M, Function}
end

# # Implicit N and M constructor
# DynamicManeuver(a::NTuple{N, Function}, b::NTuple{M, Function}
#                                     ) where {N, M} = DynamicManeuver{N, M}(a, b)

##### FUNCTIONS  ###############################################################
function calc_dV(self::DynamicManeuver, vehicle::AbstractVehicle, t::Real,
                                            dt::Real, ttot::Real, Vref::Real)
    # Here calculate change in velocity of the vehicle based on current
    # aerodynamic forces
    error("Implementation pending")
end

function calc_dw(self::KinematicManeuver, vehicle::AbstractVehicle, t::Real,
                                                         dt::Real, ttot::Real)
     # Here calculate change in angular velocity of the vehicle based on current
     # aerodynamics forces
    error("Implementation pending")
end

##### INTERNAL FUNCTIONS  ######################################################

##### END OF KINEMATICMANEUVER  ################################################

#

max_stars_repo_name: rachealerhard/FLOWUnsteady

max_issues_repo_name: rachealerhard/FLOWUnsteady

max_forks_repo_name: rachealerhard/FLOWUnsteady


text:
infile = "day19/input.txt"
rules, text = split(read(infile, String), "\n\n")
rules = split(rules, "\n")
text = split(text, "\n")

r_dict = Dict()
for rule in rules
    number, matchlist = match(r"(\d+): (.+)", rule).captures
    number = parse(Int, number)
    if contains(matchlist, "\"")
        r_dict[number] = matchlist[2]
    elseif contains(matchlist, "|")
        matchlist1, matchlist2 = split.(strip.(split(matchlist, "|")))
        matchlist1 = parse.(Int, matchlist1)
        matchlist2 = parse.(Int, matchlist2)
        r_dict[number] = (matchlist1, matchlist2)
    else
        r_dict[number] = parse.(Int, split(matchlist))
    end
end

function matches_rule(rule, text)
    global r_dict
    if text === nothing || text == ""
        return nothing
    end
    if isa(rule, Tuple)
        if (firstmatch = matches_rule(rule[1], text)) === nothing
            return matches_rule(rule[2], text)
        end
        return firstmatch
    elseif isa(rule, Array)
        rem_text = text
        for r in rule
            rem_text = matches_rule(r_dict[r], rem_text)
        end
        return rem_text
    elseif isa(rule, Char)
        if text[1] == rule
            return text[2:end]
        end
    end
end

println(count(t -> matches_rule(r_dict[0], t) == "", text))

ext: jl


text:
"""
    AutoBody(sdf,map=(x,t)->x) <: AbstractBody

    - sdf(x::AbstractVector,t::Real)::Real: signed distance function
    - map(x::AbstractVector,t::Real)::AbstractVector: coordinate mapping function

Define a geometry by its `sdf` and optional coordinate `map`. All other
properties are determined using Automatic Differentiation. Note: the `map`
is composed automatically if provided, ie `sdf(x,t) = sdf(map(x,t),t)`.
"""
struct AutoBody{F1<:Function,F2<:Function} <: AbstractBody
    sdf::F1
    map::F2
    function AutoBody(sdf,map=(x,t)->x)
        comp(x,t) = sdf(map(x,t),t)
        new{typeof(comp),typeof(map)}(comp, map)
    end
end

"""
    measure!(flow::Flow, body::AutoBody; t=0, ϵ=1)

Uses `body.sdf` and `body.map` to fill the arrays:

    `flow.μ₀`, Zeroth kernel moment
    `flow.μ₁`, First kernel moment scaled by the body normal
    `flow.V`,  Body velocity
    `flow.σᵥ`,  Body velocity divergence scaled by `μ₀-1`

at time `t` using an immersion kernel of size `ϵ`.
See [Maertens & Weymouth](https://eprints.soton.ac.uk/369635/)
"""
function measure!(a::Flow{N},body::AutoBody;t=0,ϵ=1) where N
    a.V .= 0; a.μ₀ .= 1; a.μ₁ .= 0; a.σᵥ .= 0
    sdf = a.σ
    apply_sdf!(x->body.sdf(x,t), sdf)   # distance to cell center
    for I ∈ inside(a.p)
        if abs(sdf[I])<ϵ+1   # near interface
            a.σᵥ[I] = μ₀(sdf[I],ϵ)-1
            # measure properties and fill arrays at face (i,I)
            for i ∈ 1:N
                dᵢ,n,V = measure(body,loc(i,I),t)
                m = √sum(abs2,n); dᵢ /= m; n /= m
                a.V[I,i] = V[i]
                a.μ₀[I,i] = μ₀(dᵢ,ϵ)
                a.μ₁[I,i,:] = μ₁(dᵢ,ϵ).*n
            end
        elseif sdf[I]<0      # completely inside body
            a.μ₀[I,:] .= 0
            a.σᵥ[I] = -1
        end
    end
    @inside a.σᵥ[I] = a.σᵥ[I]*div(I,a.V) # scaled divergence
    BC!(a.μ₀,zeros(SVector{N}))
end

function apply_sdf!(f,a,margin=2,stride=1)
    # strided index and signed distance function
    @inline J(I) = stride*I+oneunit(I)
    @inline sdf(I) = f(WaterLily.loc(0,J(I)) .- (stride-1)/2)

    # if the strided array is indivisible, fill it using the sdf
    dims = (size(a) .-2 ) .÷ stride
    if sum(mod.(dims,2)) != 0
        @loop a[J(I)] = sdf(I) over I ∈ CartesianIndices(dims)

    # if not, fill an array with twice the stride first
    else
        apply_sdf!(f,a,margin,2stride)

    # and only improve the values within a margin of sdf=0
        tol = stride*(√length(dims)+margin)
        @loop (d = @inbounds a[J(I)+stride*CI(mod.(I.I,2))];
               a[J(I)] = abs(d)<tol ? sdf(I) : d) over I ∈ CartesianIndices(dims)
    end
end

using ForwardDiff
"""
    measure(body::AutoBody,x,t)

ForwardDiff is used to determine the geometric properties from the `sdf`.
Note: The velocity is determined _soley_ from the optional `map` function.
"""
measure(body,x,t) = (body.sdf(x,t),                               # d
                    ForwardDiff.gradient(x->body.sdf(x,t), x),    # n=∇d
                    -ForwardDiff.derivative(t->body.map(x,t), t)) # V = -dot(X)

using LinearAlgebra: tr
"""
    curvature(A::AbstractMatrix)

Return `H,K` the mean and Gaussian curvature from `A=hessian(sdf)`.
K=tr(minor(A)) in 3D and K=0 in 2D.
"""
function curvature(A::AbstractMatrix)
    H,K = 0.5*tr(A),0
    if size(A)==(3,3)
        K = A[1,1]*A[2,2]+A[1,1]*A[3,3]+A[2,2]*A[3,3]-A[1,2]^2-A[1,3]^2-A[2,3]^2
    end
    H,K
end


text:
module Photon

using Statistics, Reexport
import Knet

include("core.jl")
include("layers/core.jl")
include("layers/conv.jl")
include("layers/recurrent.jl")
include("layers/container.jl")
include("train.jl")
include("losses.jl")

include("metrics/Metrics.jl")
@reexport using .Metrics

include("data/Data.jl")
@reexport using .Data

include("utils.jl")

export Dense, BatchNorm, Sequential, Flatten, Residual, Concurrent, Dropout,
	   Conv2D, Conv3D, Conv2DTranspose, Conv3DTranspose, output_size,
	   MaxPool2D, AvgPool2D, AdaptiveAvgPool, AdaptiveMaxPool,
	   RNN, LSTM, GRU,
	   Activation, add, forward, ContextSwitch

dir(path...) = joinpath(dirname(@__DIR__),path...)

# Re-export some of the Knet features
module K
	using Knet
	export relu, elu, selu, sigm
	export softmax, nll
	export Adam, SGD, Momentum, Nesterov, Adagrad, Rmsprop, Adadelta
	export xavier, gaussian, bilinear
end

using .K

export relu, elu, selu, sigm
export softmax, nll
export Adam, SGD, Momentum, Nesterov, Adagrad, Rmsprop, Adadelta
export xavier, gaussian, bilinear

@debug "Loaded Photon"

end


text:
function interp2D( x1Vec0, x2Vec0, yMat, x1Val, x2Val )
#UNTITLED Summary of this function goes here
#   Detailed explanation goes here

x1Vec=vec(x1Vec0)
x2Vec=vec(x2Vec0)

#Check that both x1Vec and x2Vec are n_i * 1 vectors
sizex1Vec = size(x1Vec);
sizex2Vec = size(x2Vec);
sizeyMat = size(yMat);

#println("Enter $sizex1Vec and $sizex2Vec and $sizeyMat")

# if (sizex1Vec[1]<=1) || (sizex2Vec[1]<=1) || (sizex1Vec[2]!=1) || (sizex2Vec[2]!=1)
#     error("Both x1Vec and x2Vec need to be n_i by 1 vectors. The length of the vectors may differ")
# end

#Check whether x1Vec and x2Vec are monotonic
# isx1VecMonotone = all(diff[x1Vec]>=0);
# isx2VecMonotone = all(diff[x2Vec]>=0);

# if (isx1VecMonotone*isx2VecMonotone == 0)
#     error("One of x1Vec or x2Vec is non-monotonic");
# end

n1 = sizex1Vec[1];

function nonNeg(args)
    return args>=0
end

#Find where x1val fits into x1Vec and pick out the vectors in yMat
#conditional on these values
if (x1Val<=x1Vec[1])
    x1LowerIndex = 1;
    x1UpperIndex = 2;
elseif (x1Val>=x1Vec[n1])
    x1UpperIndex = n1;
    x1LowerIndex = n1 - 1;
else #x1Val is between the lower and the upper components of x1Vec
    #Get a vect
    x1VecLessX1Val = x1Vec - x1Val;

    # Get the smallest non-negative
    posI=find(nonNeg, x1VecLessX1Val)
    (a,b)=findmin(x1VecLessX1Val[posI],1)
    x1UpperIndex=posI[b];
    x1LowerIndex = x1UpperIndex - 1.;
end

newYMatAfterInterpOnX2 = Array(Float64,2, 1);

yMatCondOnx1Lower = yMat[x1LowerIndex, :];
yMatCondOnx1Upper = yMat[x1UpperIndex, :];

intfL  = InterpIrregular(vec(x2Vec),vec(yMatCondOnx1Lower), BCnil, InterpLinear);
intfH  = InterpIrregular(vec(x2Vec),vec(yMatCondOnx1Upper), BCnil, InterpLinear);

# The last bit of interpolation!

newYMatAfterInterpOnX2[1]  =interp1extra(vec(x2Vec),vec(yMatCondOnx1Lower),x2Val,intfL);
newYMatAfterInterpOnX2[2]  =interp1extra(vec(x2Vec),vec(yMatCondOnx1Upper),x2Val,intfH);

newX1VecShort=[x1Vec[x1LowerIndex] x1Vec[x1UpperIndex] ]'

intfF  = InterpIrregular(vec(newX1VecShort),vec(newYMatAfterInterpOnX2), BCnil, InterpLinear);
yVal   = interp1extra(vec(newX1VecShort),vec(newYMatAfterInterpOnX2),x1Val,intfF); #  intfF[x1Val];

return yVal

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using Checkers

import Base.Test: record, finish, Result, Pass, Fail, Error, AbstractTestSet

function is_tiny(x)
    abs(x) < 10.0^(-6)
end

function is_huge(x)
    x > 10.0^(6)
end

function is_increasing(f::Function)
            eval(:(@test_formany -Inf<x::Float64<Inf, -Inf<y::Float64<Inf, x < y --> $f(x) < $f(y)))
end

@test isa( ( @test_exists 0<x::Float64<<1, is_tiny(x) ) , Pass )

# Borrowed from julia/test/test.jl:
# Fail/Error results do not throw errors when tests are performed,
# so *this* set of tests can test for expected fails / errors and pass.
type NoThrowTestSet <: AbstractTestSet
    results::Vector
    NoThrowTestSet(desc) = new([])
end
record(ts::NoThrowTestSet, t::Result) = (push!(ts.results, t); t)
finish(ts::NoThrowTestSet) = ts.results

results = @testset NoThrowTestSet begin
	#1 Fail
	@test_formany 0<x::Float64<10, x<9 ntests = 1000
	#2 Fail
	@test_exists -Inf<x::Int64<Inf, x>10000
	#3 Fail
	is_increasing(x -> x^2)
	#4 Error: generator for y receives incorrect lower and upper bounds
	@test_formany Inf > x > 0, x<y<10, x>0
	#5 Pass
	@test_forall i in 1:20, 1==1
	#6 Pass
	is_increasing(x -> x)
	#7 Pass
	is_increasing(x -> x^3)
end

#Now the actual expectations:
expected = [Fail, Fail, Fail, Error, Pass, Pass, Pass]
@testset "@test_... macros behave as expected" begin
    for i in eachindex(expected)
        @test isa(results[i],expected[i])
    end
end

results_2 = @testset NoThrowTestSet begin
    @test_forall 100==100
    @test_formany maxtests=100 0 < x < 10, x < 11
    @test_exists maxtests=100 0 < x < 10, x < 10
    @test_exists ntests=100 0 < x < 10, x < 10
    @test_exists ntests=100 0 < x < 10, x < 9 --> x < 10
    @test_exists maxtests=100 0 < x < 10, x < 9 --> x < 10
    @test_exists 0 < x < 10, x < 9 --> x < 10
end
expected_2 = [Error, Pass, Pass, Pass, Pass, Pass, Pass]

results_3 = @testset NoThrowTestSet begin
    #Workaround from @TotalVerb on Julia Gitter
    temp_file = tempname()
    @eval (@test_formany logto=$temp_file 0 < x < 10, x < 11)
    rm(temp_file)

    #Should pass
    @test_formany 0<<x<1, x < 2

    #Should fail...formany wants all of them
    @test_formany 0<<x<1, x < 0.5

    #Should pass
    @test_exists ntests=1000 0<<x<1, x < 2

    #This should fail, because x will veer towards 1
    @test_exists ntests=1000 0<<x<1, x < 0.5

    #This should pass. Great!
    @test_exists ntests=1000 0<x<<1, x < 0.5

    #This should pass too.
    @test_exists ntests=1000 0<x<<1, x < 0.01

end

expected_3 = [Pass, Pass, Fail, Pass, Fail, Pass, Pass]
@testset "@test_... macros behave as expected" begin
    for i in eachindex(expected_3)
        @test isa(results_3[i],expected_3[i])
    end
end


text:
"""
    gsheet_json(color::Colorant)

Returns a dictionary of values to describe a color.

# Arguments
- `color::Colorant`: color
"""
function gsheet_json(color::Colorant)
    c = convert(RGBA, color)
    return Dict(
            "red" => red(c),
            "green" => green(c),
            "blue" => blue(c),
            "alpha" => alpha(c),
        )
end

"""
    gsheet_json(range::CellIndexRange1D, dim::AbstractString)

Returns a dictionary of values to describe a CellIndexRange1D.

# Arguments
- `range::CellIndexRange1D`: cell index range
- `dim::AbstractString`: dimension
"""
function gsheet_json(range::CellIndexRange1D, dim::AbstractString)
    return Dict(
        "sheetId" => range.sheet.id,
        "dimension" => dim,
        "startIndex" => range.start_index,
        "endIndex" => range.end_index,
    )
end

"""
    gsheet_json(range::CellIndexRange2D)

Returns a dictionary of values to describe a CellIndexRange2D.

# Arguments
- `range::CellIndexRange2D`: cell index range
"""
function gsheet_json(range::CellIndexRange2D)
    return Dict(
            "sheetId" => range.sheet.id,
            "startRowIndex" => range.start_row_index,
            "endRowIndex" => range.end_row_index+1,
            "startColumnIndex" => range.start_col_index,
            "endColumnIndex" => range.end_col_index+1,
        )
end

"""
    gsheet_json(format::CellFormat)

Returns a dictionary of values to describe a Format.

# Arguments
- `format::CellFormat`: cell format
"""
function gsheet_json(format::CellFormat)
    rst = Dict{String,Any}()

    number_format = Dict{String,Any}()
    if !isnothing(format.number_format_type)  number_format["type"] = gsheet_string(format.number_format_type) end
    if !isnothing(format.number_format_pattern)  number_format["pattern"] = format.number_format_pattern end
    if length(number_format) > 0  rst["numberFormat"] = number_format end

    text_format = Dict{String,Any}()
    if !isnothing(format.text_bold)  text_format["bold"] = true end
    if !isnothing(format.text_italic)  text_format["italic"] = true end
    if !isnothing(format.text_color)  text_format["foreground_color"] = gsheet_json(format.text_color) end
    if !isnothing(format.text_font_size)  text_format["fontSize"] = format.text_font_size end
    if length(text_format) > 0  rst["textFormat"] = text_format end

    if !isnothing(format.background_color)  rst["backgroundColor"] = gsheet_json(format.background_color) end

    return rst
end


text:
function _safe_kill(gw::GitWorker, phint)
    pdat = _find_proc_reg(gw, phint)
    pid = get!(pdat, "pid", "-1")
    (pid == "-1") && return false
    if _validate_proc(gw, pid)
        return _force_kill(pid)
    end
    return false
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
# -----------------------------------------------
# Non-Linear Regression routines
# -----------------------------------------------

# #########   NON-LINEAR REGRESSION   #########

function nlinfit{T<:Number}(X::Array{T}, Y::Array{T}, regtype=Optimize(); fast=false)
    if isdefined(regtype, :intercept)
        X = regtype.intercept ? addintercept(X) : nointercept(X)
    end
    fitted_beta, converged, iter = regtype(X, Y)
    residuals = Y - regtype.modelfun(X, fitted_beta)
    # mse = Optim.minimum(optimresult)
    if regtype.robust == false || typeof(regtype.robust) == BasicVariance
        variance = BasicVariance(X, residuals, length(fitted_beta))
    elseif typeof(regtype.robust) == HCEVariance
        variance = HCEVariance(X, residuals)
    end
    if fast
        return (fitted_beta, variance)
    else
        return nlinfitresult(regtype, X, Y, fitted_beta, residuals, variance, regtype.robust, converged, iter)
    end
end

# #########   OPTIMIZE   #########

function (optimize_nlinfit::Optimize)(X::Array{T}, Y::Array{T}) where {T<:Number}
    # modelfun, beta0, method::AbsOptim=NewtonTrustRegion()
    # test__modelfun(modelfun, X[1:2,:], beta0)

    # Residual function given the model:
    residual(_beta_) = sum((Y - optimize_nlinfit.modelfun(X, _beta_)).^2)
    optimresult = optimize(residual, optimize_nlinfit.beta0, optimize_nlinfit.method)
    return Optim.minimizer(optimresult), Optim.converged(optimresult), Optim.iterations(optimresult)
end

function test__modelfun(modelfun, x, b)
    try
        modelfun(x, b)
        return true
    catch
        error("modelfun is not correctly defined")
    end
end


text:
const pixelstring =
"00000000000000000000000000000000" *
"01111111110000000111111110000000" *
"01110001111000001111001111000000" *
"01110000111000001110000111000000" *
"01110001111000001110000000000000" *
"01111111110000001110000000000000" *
"01110111100000001110000111000000" *
"01110011110011101111001111011100" *
"01110001111011100111111110011100" *
"00000000000000000000000000000000"
const pixels = reshape([UInt8(c- 48) for c in pixelstring], (32,10))'

function surroundtesting(px, i, j, step)
    if px[i,j] == 0
        return false
    end
    isize, jsize = size(px)
    if i < 1 || j < 1 || i == isize || j == jsize                         # criteria 0.both
        return false
    end
    s = Array{Int,1}(9)
    s[1] = s[9] = px[i-1,j]; s[2] = px[i-1,j+1]; s[3] = px[i,j+1]; s[4] = px[i+1,j+1]
    s[5] = px[i+1,j]; s[6] = px[i+1,j-1]; s[7] = px[i,j-1]; s[8] = px[i-1,j-1]
    b = sum(s[1:8])
    if b < 2 || b > 6                                                     # criteria 1.both
        return false
    end
    if sum([(s[i] == 0 && s[i+1] == 1) for i in 1:length(s)-1]) != 1      # criteria 2.both
        return false
    end
    if step == 1
        rightwhite = s[1] == 0 || s[3] == 0 || s[5] == 0                  # 1.3
        downwhite = s[3] == 0 || s[5] == 0 || s[7] == 0                   # 1.4
        return rightwhite && downwhite
    end
    upwhite = s[1] == 0 || s[3] == 0 || s[7] == 0                         # 2.3
    leftwhite = s[1] == 0 || s[5] == 0 || s[7] == 0                       # 2.4
    return upwhite && leftwhite
end

function zsthinning(mat)
    retmat = copy(mat)
    testmat = zeros(Int, size(mat))
    isize, jsize = size(testmat)
    needredo = true
    loops = 0
    while(needredo)
        loops += 1
        println("loop number $loops")
        needredo = false
        for n in 1:2
            for i in 1:isize, j in 1:jsize
                testmat[i,j] = surroundtesting(retmat, i, j, n) ? 1 : 0
            end
            for i in 1:isize, j in 1:jsize
                if testmat[i,j] == 1
                    retmat[i,j] = 0
                    needredo = true
                end
            end
        end
    end
    retmat
end

function asciiprint(mat)
    for i in 1:size(mat)[1]
        println(join(map(i -> i == 1 ? '#' : ' ', mat[i,:])))
    end
end

asciiprint(zsthinning(pixels))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
using DynamicHMCExamples: rel_path
using Random
Random.seed!(UInt32[0x57a97f0d, 0x1a38664c, 0x0dddb228, 0x7dbba96f])

# includenbinclude(rel_path("Estimation of a covariance matrix.ipynb"))
# @test all(ESS .≥ 500)

@testset "independent bernoulli" begin
    include(rel_path("example_independent_bernoulli.jl"))
    @test ess_α ≥ 200
end

@testset "linear regression" begin
    include(rel_path("example_linear_regression.jl"))
    @test length(ess) == 4
end

@testset "logistic regression" begin
    include(rel_path("example_logistic_regression.jl"))
    @test all(ess .≥ 400)
end

@testset "multinomial logistic regression" begin
    include(rel_path("example_multinomial_logistic_regression.jl"))
    @test all(ess .≥ 200)
end


text:
let
    solver = NelderMead()
    T = typeof(solver)
    trace = OptimizationTrace{T}()
    push!(trace,OptimizationState{T}(1,1.0,1.0,Dict()))
    push!(trace,OptimizationState{T}(2,1.0,1.0,Dict()))
    @test length(trace) == 2
    @test trace[end].iteration == 2

    prob = Optim.UnconstrainedProblems.examples["Rosenbrock"]
    f_prob = prob.f
    for g_free in (NelderMead(), SimulatedAnnealing())
        res = Optim.optimize(f_prob, prob.initial_x, method=g_free)

        io = IOBuffer()
        show(io, res)
        s = takebuf_string(io)

        lines = split(s, '\n')
        @test lines[1] == "Results of Optimization Algorithm"
        @test startswith(lines[2], " * Algorithm: ")
        @test startswith(lines[3], " * Starting Point: ")
        @test startswith(lines[4], " * Minimizer: [")
        @test startswith(lines[5], " * Minimum: ")
        @test startswith(lines[6], " * Iterations: ")
        @test startswith(lines[7], " * Convergence: ")
        if res.method == "Nelder-Mead"
            @test startswith(lines[8], "   *  √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: ")
            @test startswith(lines[9], "   * Reached Maximum Number of Iterations: ")
            @test startswith(lines[10], " * Objective Function Calls: ")
        elseif res.method == "Simulated Annealing"
            @test startswith(lines[8], "   * |x - x'| < ")
            @test startswith(lines[9], "   * |f(x) - f(x')| / |f(x)| < ")
            @test startswith(lines[10], "   * |g(x)| < ")
            @test startswith(lines[11], "   * Reached Maximum Number of Iterations: ")
            @test startswith(lines[12], " * Objective Function Calls: ")
        end
    end

    res = Optim.optimize(prob.f, prob.g!, prob.initial_x, method=LBFGS())

    io = IOBuffer()
    show(io, res)
    s = takebuf_string(io)

    lines = split(s, '\n')
    @test lines[1] == "Results of Optimization Algorithm"
    @test startswith(lines[2], " * Algorithm: ")
    @test startswith(lines[3], " * Starting Point: ")
    @test startswith(lines[4], " * Minimizer: [")
    @test startswith(lines[5], " * Minimum: ")
    @test startswith(lines[6], " * Iterations: ")
    @test startswith(lines[7], " * Convergence: ")
    @test startswith(lines[8], "   * |x - x'| < ")
    @test startswith(lines[9], "   * |f(x) - f(x')| / |f(x)| < ")
    @test startswith(lines[10], "   * |g(x)| < ")
    @test startswith(lines[11], "   * Reached Maximum Number of Iterations: ")
    @test startswith(lines[12], " * Objective Function Calls: ")
    @test startswith(lines[13], " * Gradient Calls: ")
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Distributions, StatsBase, IntervalSets, ValueShapes, ArraysOfArrays

@testset "hierarchial_density" begin
    let
        parent_density = NamedTupleDist(
            foo = Exponential(3.5),
            bar = Normal(2.0, 1.0)
        )

        f = v -> NamedTupleDist(baz = fill(Normal(v.bar, v.foo), 3))

        @test typeof(@inferred HierarchicalDensity(f, parent_density)) <: HierarchicalDensity

        hd = HierarchicalDensity(f, parent_density)

        @test @inferred(sampler(hd)) == BAT.HierarchicalDensitySampler(hd)
        @test @inferred(rand(sampler(hd))) isa AbstractVector{<:Real}
        @test @inferred(varshape(hd)) == NamedTupleShape(foo = ScalarShape{Real}(), bar = ScalarShape{Real}(), baz = ArrayShape{Real}(3))

        @test @inferred(BAT.eval_logval_unchecked(hd, [2.7, 4.3, 8.7, 8.7, 8.7])) ≈ logpdf(parent_density.foo, 2.7) + logpdf(parent_density.bar, 4.3) + 3 * logpdf(Normal(4.3, 2.7), 8.7)
        @test @inferred(BAT.eval_logval_unchecked(hd, (foo = 2.7, bar = 4.3, baz = fill(8.7, 3)))) ≈ logpdf(parent_density.foo, 2.7) + logpdf(parent_density.bar, 4.3) + 3 * logpdf(Normal(4.3, 2.7), 8.7)

        @test @inferred(BAT.var_bounds(hd)) == BAT.HierarchicalDensityBounds(hd)

        hd_bounds = @inferred(BAT.var_bounds(hd))
        @test all(in.(nestedview(@inferred rand(sampler(hd), 10^3)), Ref(hd_bounds)))
        @test @inferred(fill(-1.0, totalndof(hd)) in hd_bounds) == false

        @test (@inferred(BAT.renormalize_variate!(fill(NaN, totalndof(hd)), hd_bounds, fill(-1.0, totalndof(hd)))) in hd_bounds) == false

        posterior = PosteriorDensity(LogDVal(0), hd)
        samples = bat_sample(posterior, MCMCSampling(sampler = MetropolisHastings(), trafo = NoDensityTransform(), nsteps = 10^5)).result
        isapprox(cov(unshaped.(samples)), cov(hd), rtol = 0.05)
    end

    let
        parent_density = BAT.DistributionDensity(NamedTupleDist(
            foo = 2..4,
            bar = Normal(2.0, 1.0)
        ), bounds_type = BAT.reflective_bounds)

        f = v -> BAT.DistributionDensity(NamedTupleDist(baz = fill(Normal(v.bar, v.foo), 3)), bounds_type = BAT.reflective_bounds)

        hd = HierarchicalDensity(f, parent_density)

        hd_bounds = @inferred(BAT.var_bounds(hd))
        @test @inferred(fill(1.5, totalndof(hd)) in hd_bounds) == false
        @test @inferred(BAT.renormalize_variate!(fill(NaN, totalndof(hd)), hd_bounds, fill(1.5, totalndof(hd)))) == [2.5, 1.5, 1.5, 1.5, 1.5]
        @test (@inferred(BAT.renormalize_variate!(fill(NaN, totalndof(hd)), hd_bounds, fill(1.5, totalndof(hd)))) in hd_bounds) == true
    end

    let
        hd = HierarchicalDensity(
            v -> NamedTupleDist(b = Normal(v.a, 1.2)),
            NamedTupleDist(a = Normal(2.3, 1.9))
        )

        cov_expected = [1.9^2 1.9^2; 1.9^2 1.9^2 + 1.2^2]

        @test isapprox(cov(hd), cov_expected, rtol = 0.05)
        @test isapprox(mean(nestedview(rand(sampler(hd), 10^5))), [2.3, 2.3], rtol = 0.05)
    end
end


text:
function JacobiPreconditioner(arg0::Vector{jdouble}, arg1::jboolean)
    return JacobiPreconditioner((Vector{jdouble}, jboolean), arg0, arg1)
end

function create(::Type{JacobiPreconditioner}, arg0::RealLinearOperator)
    return jcall(JacobiPreconditioner, "create", JacobiPreconditioner, (RealLinearOperator,), arg0)
end

function get_column_dimension(obj::JacobiPreconditioner)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_row_dimension(obj::JacobiPreconditioner)
    return jcall(obj, "getRowDimension", jint, ())
end

function operate(obj::JacobiPreconditioner, arg0::RealVector)
    return jcall(obj, "operate", RealVector, (RealVector,), arg0)
end

function sqrt(obj::JacobiPreconditioner)
    return jcall(obj, "sqrt", RealLinearOperator, ())
end