text:
export GeneralizedLinearRegression, GLR,
        LinearRegression, RidgeRegression,
        LassoRegression, ElasticNetRegression,
        LADRegression, LogisticRegression,
        MultinomialRegression, RobustRegression,
        HuberRegression, QuantileRegression

"""
    GeneralizedLinearRegression{L<:Loss, P<:Penalty}

Generalized Linear Regression (GLR) model with objective function:

``L(y, Xθ) + P(θ)``

where `L` is a loss function, `P` a penalty, `y` is the vector of observed
response, `X` is the feature matrix and `θ` the vector of parameters.

Special cases include:

* **OLS regression**:      L2 loss, no penalty.
* **Ridge regression**:    L2 loss, L2 penalty.
* **Lasso regression**:    L2 loss, L1 penalty.
* **Logistic regression**: Logit loss, [no,L1,L2] penalty.
"""
@with_kw mutable struct GeneralizedLinearRegression{L<:Loss, P<:Penalty}
    # Parameters that can be tuned
    loss::L                  = L2Loss()    # L(y, ŷ=Xθ)
    penalty::P               = NoPenalty() # P(θ)
    fit_intercept::Bool      = true        # add intercept ? def=true
    penalize_intercept::Bool = false
end

const GLR = GeneralizedLinearRegression

getc(g::GLR) = getc(g.loss)
getc(g::GLR, y) = getc(g.loss, y)

## Specific constructors

"""
$SIGNATURES

Objective function: ``|Xθ - y|₂²/2``.
"""
LinearRegression(; fit_intercept::Bool=true) = GLR(fit_intercept=fit_intercept)

"""
$SIGNATURES

Objective function: ``|Xθ - y|₂²/2 + λ|θ|₂²/2``.
"""
function RidgeRegression(λ::Real=1.0; lambda::Real=λ, fit_intercept::Bool=true,
                         penalize_intercept::Bool=false)
    check_pos(lambda)
    GLR(penalty=lambda*L2Penalty(),
        fit_intercept=fit_intercept,
        penalize_intercept=penalize_intercept)
end

"""
$SIGNATURES

Objective function: ``|Xθ - y|₂²/2 + λ|θ|₁``.
"""
function LassoRegression(λ::Real=1.0; lambda::Real=λ, fit_intercept::Bool=true,
                         penalize_intercept::Bool=false)
    check_pos(lambda)
    GLR(penalty=lambda*L1Penalty(),
        fit_intercept=fit_intercept,
        penalize_intercept=penalize_intercept)
end

"""
$SIGNATURES

Objective function: ``|Xθ - y|₂²/2 + λ|θ|₂²/2 + γ|θ|₁``.
"""
function ElasticNetRegression(λ::Real=1.0, γ::Real=1.0;
                              lambda::Real=λ, gamma::Real=γ,
                              fit_intercept::Bool=true,
                              penalize_intercept::Bool=false)
    check_pos.((lambda, gamma))
    GLR(penalty=lambda*L2Penalty()+gamma*L1Penalty(),
        fit_intercept=fit_intercept,
        penalize_intercept=penalize_intercept)
end

"""
$SIGNATURES

Helper function for objectives which can have 0/L1/L2 penalties.
"""
function _l1l2en(lambda, gamma, penalty, r)
    check_pos.((lambda, gamma))
    penalty ∈ (:l1, :l2, :en, :none) ||
        throw(ArgumentError("Unrecognised penalty for the $r: '$penalty' " *
                            "(expected none/l1/l2/en)"))
    penalty = if penalty == :none
       NoPenalty()
    elseif penalty == :l1
        lambda * L1Penalty()
    elseif penalty == :l2
        lambda * L2Penalty()
    else
        lambda * L2Penalty() + gamma * L1Penalty()
    end
    return penalty
end

"""
$SIGNATURES

Objective function: ``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is either the
logistic loss in the binary case or the multinomial loss otherwise.
"""
function LogisticRegression(λ::Real=1.0, γ::Real=0.0;
                            lambda::Real=λ, gamma::Real=γ,
                            penalty::Symbol=iszero(gamma) ? :l2 : :en,
                            fit_intercept::Bool=true,
                            penalize_intercept::Bool=false,
                            multi_class::Bool=false,
                            nclasses::Integer=0)
    penalty = _l1l2en(lambda, gamma, penalty, "Logistic regression")
    loss = LogisticLoss()
    if nclasses > 2     # number of classes is explicitly specified
        loss = MultinomialLoss(nclasses)
    elseif multi_class  # number of classes will be inferred from data
        loss = MultinomialLoss()
    end
    GLR(loss=loss,
        penalty=penalty,
        fit_intercept=fit_intercept,
        penalize_intercept=penalize_intercept)
end

"""
$SIGNATURES

Objective function: ``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is the
multinomial loss.
"""
MultinomialRegression(a...; kwa...) =
    LogisticRegression(a...; multi_class=true, kwa...)

# ========

"""
$SIGNATURES

Objective function: ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` where ρ is a given function
on the residuals.
"""
function RobustRegression(ρ::RobustRho=HuberRho(0.1), λ::Real=1.0, γ::Real=0.0;
                          rho::RobustRho=ρ, lambda::Real=λ, gamma::Real=γ,
                          penalty::Symbol=iszero(gamma) ? :l2 : :en,
                          fit_intercept::Bool=true,
                          penalize_intercept::Bool=false)
    penalty = _l1l2en(lambda, gamma, penalty, "Robust regression")
    GLR(loss=RobustLoss(rho),
        penalty=penalty,
        fit_intercept=fit_intercept,
        penalize_intercept=penalize_intercept)
end

"""
$SIGNATURES

Huber Regression with objective:

``∑ρ(Xθ - y) + λ|θ|₂²/2 + γ|θ|₁``

Where `ρ` is the Huber function `ρ(r) = r²/2``  if `|r|≤δ` and
`ρ(r)=δ(|r|-δ/2)` otherwise.
"""
function HuberRegression(δ::Real=0.5, λ::Real=1.0, γ::Real=0.0;
                         delta::Real=δ, lambda::Real=λ, gamma::Real=γ,
                         penalty::Symbol=iszero(gamma) ? :l2 : :en,
                         fit_intercept::Bool=true,
                         penalize_intercept::Bool=false)
    return RobustRegression(HuberRho(delta), lambda, gamma;
                            penalty=penalty, fit_intercept=fit_intercept,
                            penalize_intercept=penalize_intercept)
end

"""
$SIGNATURES

Quantile Regression with objective:

``∑ρ(Xθ - y) + λ|θ|₂²/2 + γ|θ|₁``

Where `ρ` is the check function `ρ(r) = r(δ - 1(r < 0))`.
"""
function QuantileRegression(δ::Real=0.5, λ::Real=1.0, γ::Real=0.0;
                            delta::Real=δ, lambda::Real=λ, gamma::Real=γ,
                            penalty::Symbol=iszero(gamma) ? :l2 : :en,
                            fit_intercept::Bool=true,
                            penalize_intercept::Bool=false)
    return RobustRegression(QuantileRho(delta), lambda, gamma;
                            penalty=penalty, fit_intercept=fit_intercept,
                            penalize_intercept=penalize_intercept)
end

"""
$SIGNATURES

Least Absolute Deviation regression with objective:

``|Xθ - y|₁ + λ|θ|₂²/2 + γ|θ|₁``

This is a specific type of Quantile Regression with `δ=0.5` (median).
"""
function LADRegression(λ::Real=1.0, γ::Real=0.0;
                       lambda::Real=λ, gamma::Real=γ,
                       penalty::Symbol=iszero(gamma) ? :l2 : :en,
                       fit_intercept::Bool=true, penalize_intercept::Bool=false)
    return QuantileRegression(0.5, lambda, gamma;
                              penalty=penalty, fit_intercept=fit_intercept,
                              penalize_intercept=penalize_intercept)
end


text:
function f(x::AbstractIntVar,alpha::Float64)
    return 1/(assignedValue(x)+1)^alpha
end

"""
    struct CPReward <: AbstractReward end
This is the default reward, that will be used if no custom reward is specified when constructing the `LearnedHeuristic`.
"""
mutable struct CPReward <: AbstractReward
    value::Float32
end

CPReward(model::CPModel) = CPReward(0)

"""
    set_reward!(::StepPhase, lh::LearnedHeuristic{DefaultReward, A}, model::CPModel, symbol::Union{Nothing, Symbol})
Change the "current_reward" attribute of the LearnedHeuristic at the StepPhase.
"""
function set_reward!(::Type{StepPhase}, lh::LearnedHeuristic{SR, CPReward, A}, model::CPModel, symbol::Union{Nothing, Symbol}) where {
    SR <: AbstractStateRepresentation,
    A <: ActionOutput
}
    nothing
end

"""
    set_reward!(::DecisionPhase, lh::LearnedHeuristic{DefaultReward, O}, model::CPModel)
Change the current reward at the DecisionPhase. This is called right before making the next decision, so you know you have the very last state before the new decision
and every computation like fixPoints and backtracking has been done.
"""
function set_reward!(::Type{DecisionPhase}, lh::LearnedHeuristic{SR, CPReward, A}, model::CPModel) where {
    SR <: AbstractStateRepresentation,
    A <: ActionOutput
}
    nothing
end

"""
    set_reward!(::EndingPhase, lh::LearnedHeuristic{DefaultReward, A}, model::CPModel, symbol::Union{Nothing, Symbol})
Increment the current reward at the EndingPhase. Called when the search is finished by an optimality proof or by a limit in term of nodes or
in terms of number of solution. This is useful to add some general results to the reward like the number of ndoes visited during the episode for instance.
"""
function set_reward!(::Type{EndingPhase}, lh::LearnedHeuristic{SR, CPReward, A}, model::CPModel, symbol::Union{Nothing, Symbol}) where {
    SR <: AbstractStateRepresentation,
    A <: ActionOutput
}
    alpha = 0.5
    if symbol == :FoundSolution
        if isnothing(model.objective)
            lh.reward.value += 1
        else
            lh.reward.value += f(model.objective, alpha)
        end
    else
        lh.reward.value -= length(branchable_variables(model))/nb_boundvariables(model) - 1
    end
end

ext: jl


text:
module TestCore

include("preamble.jl")

@testset begin
    @test implements(push!, [])
    @test implements(push!, Dict(:a=>1))
    @test !implements(push!, ())
    @test !implements(push!, (a=1,))
    @test !implements(push!, SVector(0))
end

end  # module


text:
function viterbi(
        a::AbstractVector,
        A::AbstractArray,
        L::AbstractArray,
        track_length::AbstractArray;
        logl = false
    )
    T1 = Array{Float64}(undef, size(L))
    T2 = Array{Int}(undef, size(L))
    z = Matrix{Int}(undef, size(L, 1), size(L, 3))
    c = Matrix{Float64}(undef, size(z))
    if logl
        error("logviterbi! not implemented yet.")
#         viterbilog!(T1, T2, z, a, A, L)
    else
        warn_logl(L)
        viterbi!(T1, T2, z, a, A, L, c, track_length)
    end
    z
end

function viterbi!(
        T1::AbstractArray,
        T2::AbstractArray,
        z::AbstractArray,
        a::AbstractVector,
        A::AbstractArray,
        L::AbstractArray,
        c::AbstractArray,
        track_length::AbstractArray
    )
    @argcheck size(T1, 1) == size(T2, 1) == size(L, 1) == size(z, 1) == size(c, 1)
    @argcheck size(T1, 2) == size(T2, 2) == size(L, 2) == size(a, 1) == size(A, 1) == size(A, 2)
    @argcheck size(T1, 3) == size(T2, 3) == size(L, 3) == size(z, 2) == size(c, 2)

    T, K, N = size(L)
    (T == 0) && return

    fill!(T1, 0.0)
    fill!(T2, 0)
    fill!(c, 0.0)

    for i in 1:N

        for j in 1:K
            T1[1, j, i] = a[j] * L[1, j, i]
            c[1, i] += T1[1, j, i]
        end

        for j in 1:K
            T1[1, j, i] /= c[1, i]
        end

        @inbounds for t = 2:track_length[i]
            for j1 in 1:K
                # TODO: If there is NaNs in T1 this may
                # stay to 0 (NaN > -Inf == false).
                # Hence it will crash when computing z[t].
                # Maybe we should check for NaNs beforehand ?
                amax = 0
                vmax = -Inf

                for j2 in 1:K
                    v = T1[t-1, j2, i] * A[j2, j1]
                    if v > vmax
                        amax = j2
                        vmax = v
                    end
                end

                T1[t, j1, i] = vmax * L[t, j1, i]
                T2[t, j1, i] = amax
                c[t, i] += T1[t, j1, i]
            end

            for j in 1:K
                T1[t, j, i] /= c[t, i]
            end
        end
    end

    z[track_length[i], i] = argmax(T1[track_length[i], :, i])
    for t in reverse(1:track_length[i]-1)
        z[t, i] = T2[t+1, z[t+1, i], i]
    end
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
module Taxize
  using Requests, JSON, LightXML

  export itisPing, gettaxonomicranknamefromtsn

  """
  itisPing - ping the ITIS server and get some metadata on the service

  Usage:

    Taxize.itisPing()
  """
  function itisPing()
    temp = Requests.get("http://www.itis.gov/ITISWebService/services/ITISService/getDescription")
    return temp.data
  end

  """
  Returns the kingdom and rank information for the TSN.

  :param tsn: TSN for a taxonomic group (numeric)

  Usage:
  Taxize.gettaxonomicranknamefromtsn(tsn = 202385)
  """
  function gettaxonomicranknamefromtsn(tsn)
    url = "http://www.itis.gov/ITISWebService/services/ITISService/getTaxonomicRankNameFromTSN"
    payload = {"tsn" => tsn}
    out = Requests.get(url; query = payload)
    check_status(out)
    dat = out.data
    return replace(replace(dat, "\r", ""), "\n", "")
  end

  """
  Uses the Global Names Index to parse scientific names

  :param names: List of scientific names.

  Usage:
  Taxize.gni_parse(names = ("Cyanistes caeruleus","Helianthus annuus"))
  """
  function gni_parse(names)
    url = "http://gni.globalnames.org/parsers.json"
    names = join(names, "|")
    out = Requests.get(url; query = {"names" => names})
    check_status(out)
    return out.data
  end

  function check_status(response_obj)
    if response_obj.status > 200
      @sprintf("Call error, http status code: %s", response_obj.status)
    end
  end

end


text:
#!/usr/bin/env julia

using VibronicToolkit

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "--conf"
        metavar = "FILE"
        help = "path to config file"
        required = true
    "--beta"
        metavar = "T"
        help = "reciprocal temperature"
        arg_type = Float64
        required = true
    "--basis-size"
        metavar = "N"
        help = "single-mode basis size"
        arg_type = Int
        required = true
end
c = parse_args(ARGS, s, as_symbols=true)

sys = read(c[:conf], DenseSystem)
beta = c[:beta]
basis_size = c[:basis_size]

sos = SumOverStates(sys, beta, basis_size)

println("Z: $(sos.Z)")
println("E: $(sos.E)")
println("Cv: $(sos.Cv)")


text:
const G = 9.8
const CARMAXSTEER = pi/6

struct Car
    # coordinates
    rlx::Float32
    rly::Float32
    rrx::Float32
    rry::Float32
    flx::Float32
    fly::Float32
    frx::Float32
    fry::Float32
    minx::Float32
    miny::Float32
    maxx::Float32
    maxy::Float32
    cx::Float32
    cy::Float32
    fcx::Float32 # front center
    fcy::Float32
    rcx::Float32 # rear center
    rcy::Float32

    # physical properties
    full_steering_dt::Float32
    maxsteer::Float32
    wheelbase::Float32
    fronttrack::Float32
    steering_speed::Float32
    wheel_static_friction_coefficient::Float32

    # state (other than coordinates)
    angle::Float32
    speed::Float32
    steer::Float32
    accel::Float32

    function Car(rcx,rcy,angle,steer,accel,speed)
        full_steering_dt = 3.0
        maxsteer = CARMAXSTEER
        wheelbase = 2.7
        fronttrack = 1.65
        steering_speed = 2.0*maxsteer / full_steering_dt
        wheel_static_friction_coefficient = 0.7

        ca9 = cos(pi / 2 - angle)
        sa9 = sin(pi / 2 - angle)
        ca = cos(angle)
        sa = sin(angle)

        fa2 = fronttrack / 2.0

        rlx = rcx - fa2 * ca9
        rly = rcy + fa2 * sa9
        rrx = rcx + fa2 * ca9
        rry = rcy - fa2 * sa9

        fx = rcx + wheelbase * ca
        fy = rcy + wheelbase * sa

        flx = fx - fa2 * ca9
        fly = fy + fa2 * sa9
        frx = fx + fa2 * ca9
        fry = fy - fa2 * sa9

        minx = min(rlx, rrx, flx, frx)
        maxx = max(rlx, rrx, flx, frx)
        miny = min(rly, rry, fly, fry)
        maxy = max(rly, rry, fly, fry)

        new(
            rlx, rly, rrx, rry, flx, fly, frx, fry,
            minx,
            miny,
            maxx,
            maxy,
            0.5 * (minx + maxx),
            0.5 * (miny + maxy),
            0.5 * (flx + frx),
            0.5 * (fly + fry),
            rcx, rcy,

            full_steering_dt,
            maxsteer,
            wheelbase,
            fronttrack,
            steering_speed,
            wheel_static_friction_coefficient,

            angle, speed, steer, accel,
        )
    end
end

function get_steer(c::Car, Rw::Float64)
    Rr = Rw - c.fronttrack / 2.0
    steer = pi/2.0 - atan(Rr/c.wheelbase)
    return steer
end

function get_max_safe_steer(c::Car)
    max_safe_steer = c.maxsteer
    if abs(c.speed) > 1.0e-2
        max_safe_steer = abs(atan(c.wheel_static_friction_coefficient*G*c.wheelbase/c.speed/c.speed))
    end
    return min(max_safe_steer, c.maxsteer)
end

function move(c::Car, dt::Float64, steering::Float64, accel::Float64)
    @assert dt > 1.0e-3

    maxsteer = get_max_safe_steer(c)
    if steering > maxsteer
        steering = maxsteer
    elseif steering < -maxsteer
        steering = -maxsteer
    end

    if steering - c.steer < -c.steering_speed*dt
        steering = c.steer - c.steering_speed*dt
    elseif steering - c.steer > c.steering_speed*dt
        steering = c.steer + c.steering_speed*dt
    end

    angle = c.angle
    rcx = c.rcx
    rcy = c.rcy

    arclen = c.speed * dt + 0.5 * accel * dt * dt
    if abs(steering) > 1.0e-4
        Rr = abs(c.wheelbase * tan(pi / 2.0 - steering))

        Rw = (Rr + c.fronttrack / 2.0)

        rx = c.rcx - Rw * cos(pi / 2.0 - c.angle) * sign(steering)
        ry = c.rcy + Rw * sin(pi / 2.0 - c.angle) * sign(steering)

        dangle = arclen / Rr * sign(steering)

        # translate to origin
        ox = c.rcx - rx
        oy = c.rcy - ry
        # rotate
        cda = cos(dangle)
        sda = sin(dangle)
        rcx = rx + ox * cda - oy * sda
        rcy = ry + ox * sda + oy * cda
        angle += dangle
    else
        # asumming steer of 0.0
        rcx += arclen * cos(c.angle)
        rcy += arclen * sin(c.angle)
    end

    speed = c.speed + accel * dt
    if angle < -pi
        angle += 2.0*pi
    elseif angle > pi
        angle -= 2.0*pi
    end

    Car(rcx,rcy,angle,steering,accel,speed)
end

function get_plot(c::Car)
    [c.flx,c.frx,c.rrx,c.rlx,c.flx],
    [c.fly,c.fry,c.rry,c.rly,c.fly]
end

getvxvy(c::Car) =
    cos(c.angle)*c.speed, sin(c.angle)*c.speed