text:
#=
Functions that automatically create Goppa codes.

Author: Ivan A. Moreno Soto
Last updated: 2019/July/15
=#

module GoppaCode

export
    generateGoppaCode,
    decode!

include("./Polynomials.jl")
using .Polynomials

include("./ExtensionField.jl")
using .ExtensionField

include("./FiniteFieldMatrix.jl")
using .FiniteFieldMatrix

include("./Misc.jl")
using .Misc

using StatsBase
using LinearAlgebra

using DelimitedFiles

"""
    selectCodeSupport(n::Int, m::Int)

Randomly selects and returns n distinct elements from the
extension field 2^m.
"""
function selectCodeSupport(n::Int, m::Int)
    return sample(Array{UInt32, 1}(0:(2^m - 1)), n, replace=false)
end

"""
    getY!(codeSupport, t, g)

Creates the Y matrix needed to compute the parity check matrix H. Basically
computes a Vandermonde matrix with the code support.
"""
function getY!(codeSupport, t, g)
    n = length(codeSupport)
    Y = zeros(UInt32, (t, n))

    for i in 1:n
        Y[1, i] = 1
    end

    for i in 1:n
        for j in 2:t
            Y[j, i] = power(codeSupport[i], (j-1), g)
        end
    end

    return Y
end

"""
    getZ!(goppaPolynomial, moduloPolynomial, codeSupport)

Creates the Z matrix needed to compute the parity check matrix H.
"""
function getZ!(goppaPolynomial, moduloPolynomial, codeSupport)
    n = length(codeSupport)
    Z = zeros(UInt32, (n, n))

    for i in 1:n
        Z[i, i] = inverse(evaluate!(codeSupport[i], goppaPolynomial, moduloPolynomial), moduloPolynomial)
    end

    return Z
end

"""
    toBits!(matrix, t, n, m)

Takes a matrix of dimension txn and returns a binary matrix
with m-bits representations of every element spread vertically.
"""
function toBits!(matrix, t, n, m)
    bits = falses(m*t, n)

    for i in 1:t
        for j in 1:n
            bits[((i - 1) * m + 1):(i * m), j] = digits(matrix[i, j], base=2, pad=m)
        end
    end

    return bits
end

"""
    generateGoppaCode(n::Int, t::Int, m::Int)

Generates a (n, k) irreducible linear Goppa code with minimum
distance d >= 2t + 1. Returns the generator matrix G, the goppa
polynomial, the code support, and the number of rows k of G.
"""
function generateGoppaCode(n::Int, t::Int, m::Int)
    # Get polynomials
    moduloPolynomial = primitivePolynomials[m]
    goppaPolynomial = randomPoly(t, m)
    # Get code support
    codeSupport = selectCodeSupport(n, m)

    # Compute H = YZ
    Y = getY!(codeSupport, t, moduloPolynomial)
    Z = getZ!(goppaPolynomial, moduloPolynomial, codeSupport)
    T = multiplyPolyMatrices!(Y, Z, moduloPolynomial)

    H = toBits!(T, t, n, m)
    gaussianEliminationColumnPivoting!(H, m*t, codeSupport)

    # Compute G
    k = n - m*t
    ATransposed = transpose(H[:, (m*t + 1):end])
    G = hcat(ATransposed, Matrix{UInt32}(I, k, k))

    return G, goppaPolynomial, codeSupport, k
end

"""
    syndrome!(msg, codeSupport, goppaPolynomial, m)

Computes the syndrome of a vector using a binary irreducible goppa code.
"""
function syndrome!(msg, codeSupport, goppaPolynomial, m)
    syndrome = Array{UInt32, 1}([0])
    for (c, γ) in zip(msg, codeSupport)
        if c == 1
            currentInverse = inverseRing!(Array{UInt32, 1}([γ, 1]), goppaPolynomial, m)
            syndrome = addRing!(syndrome, currentInverse)
        end
    end
    return syndrome
end

"""
    syndromeSquareRoot!(syndrome, goppaPolynomial, m)

Computes the square root of syndrome + x mod goppaPolynomial.
"""
function syndromeSquareRoot!(syndrome, goppaPolynomial, m)
    # Computing T(x) + x.
    syndromeInverse = inverseRing!(syndrome, goppaPolynomial, m)
    h = addRing!(syndromeInverse, PolynomialConstants['x'])

    # Computing the square root of x mod goppaPolynomial.
    goppaEven, goppaOdd = splitPolynomial!(goppaPolynomial, m)
    sx = multiplyRing!(goppaEven, inverseRing!(goppaOdd, goppaPolynomial, m), m)

    # Computing the odd and even parts of the square root of T(x) + x.
    evenh = [multiplyRing!(power(coeff, 2^(m-1), primitivePolynomials[m]), monomial(i-1), m) for (i, coeff) in enumerate(h[1:2:end])]
    oddh = [multiplyRing!(multiplyRing!([power(coeff, 2^(m-1), primitivePolynomials[m])], monomial(i-1), m), sx, m) for (i, coeff) in enumerate(h[2:2:end])]
    if length(oddh) == 0
        oddh = [copy(PolynomialConstants['0'])]
    end

    # Computing the square root by adding the even and odd parts.
    τ = copy(PolynomialConstants['0'])
    for term in vcat(evenh, oddh)
        τ = addRing!(τ, term)
    end

    return τ
end

"""
    halfwayXGCD!(goppaPolynomial, tau, t, m)

Runs the extended eucliedean algorithm halfway to compute the least
degree polynomials that can define the error locator polynomial for
a codeword with errors.
"""
function halfwayXGCD!(goppaPolynomial, tau, t, m)
    r = copy(goppaPolynomial)
    R = copy(tau)
    s = Array{UInt32, 1}([1])
    S = Array{UInt32, 1}([0])
    u = Array{UInt32, 1}([0])
    U = Array{UInt32, 1}([1])

    while length(R) - 1 > div(t, 2)
        q = divRing!(r, R, m)

        x = R; R = addRing!(r, multiplyRing!(q, R, m)); r = x
        x = S; S = addRing!(s, multiplyRing!(q, S, m)); s = x
        x = U; U = addRing!(u, multiplyRing!(q, U, m)); u = x
    end

    return R, U
end

"""
    findRootsErrorLocatorPolynomial!(sigma, codeSupport, m)

Finds the roots of sigma by evaluating every element of the code support.
"""
function findRootsErrorLocatorPolynomial!(sigma, codeSupport, m)
    totalErrors = 0
    errors = []
    for (i, gamma) in enumerate(codeSupport)
        if evaluate!(gamma, sigma, primitivePolynomials[m]) == 0
            push!(errors, i)
            totalErrors += 1
        end
    end

    return errors, totalErrors
end

"""
    decode!(msg, goppaPolynomial, codeSupport, m, t)

Given a msg with noise, finds a codeword in the goppa code defined
by the goppaPolynomial and codeSupport in GF(2^m).
"""
function decode!(msg, goppaPolynomial, codeSupport, m, t)
    # Compute syndrome
    msgSyndrome = syndrome!(msg, codeSupport, goppaPolynomial, m)
    if msgSyndrome == PolynomialConstants['0']
        throw(InvalidCodewordError())
    end

    tau = syndromeSquareRoot!(msgSyndrome, goppaPolynomial, m)

    # Using the extended euclidean algorithm to find the error locator polynomial.
    a, b = halfwayXGCD!(goppaPolynomial, tau, t, m)
    sigma = addRing!(multiplyRing!(a, a, m), multiplyRing!(PolynomialConstants['x'], multiplyRing!(b, b, m), m))
    sigma = makePolynomialMonic!(sigma, m)

    # Finding the roots of the error locator polynomial.
    errors, totalErrors = findRootsErrorLocatorPolynomial!(sigma, codeSupport, m)

    # Checking that it was a valid cyphertext.
    if totalErrors < t
        throw(InvalidCodewordError())
    end

    # Correcting the t errors.
    for position in errors
        msg[position] = (msg[position] + 1) % 2
    end

    return msg
end

end # GoppaCode

max_stars_repo_name: rexemin/McEliece

max_issues_repo_name: rexemin/McEliece

max_forks_repo_name: rexemin/McEliece


text:
module DistributedSolver

export dlinsolve, prbgs

using Distributed
@everywhere using LinearAlgebra, SparseArrays

include("./pgmres.jl")
include("./prbgs.jl")

function dlinsolve(A::SparseMatrixCSC, x0::Vector{Float64}, b::Vector{Float64})

end

end

ext: jl


text:
import SearchBall, ParticleSwarmOptimizer, DelimitedFiles, Profile, ProfileView
const sb = SearchBall
const pso = ParticleSwarmOptimizer
const df = DelimitedFiles

import Random
Random.seed!(1234);

myconfig = sb.get_config(["-f", "/home/tobias/doc/SearchBall/SearchBall.jl/configs/2dlinNN.config"])
sb.get_strategy(myconfig).init_strategies(myconfig)
@code_warntype sb.get_strategy(myconfig).calc_myPlayers(myconfig, 0.1, :default, nothing)
@code_warntype sb.update_game!(myconfig.ws, 0.1, myconfig)


text:
# Map classificator for Unique values
struct Unique <: AbstractUniqueMapClassificator
    uniq::Union{Vector{String}, Vector{Symbol}, Nothing}
    Unique(uniq::Union{Vector{String}, Vector{Symbol}, Nothing} = nothing) = return new(uniq);
end

# Structure for Unique Values  Map Classification
struct UniqueMapClassification <: AbstractMapClassification
    mcr::AbstractMapClassificator
    k::Int
    group::Vector{Int}
    ngroup::Vector{Int}
    grouplabs::Vector{String}
end

function mapclassify(mcr::AbstractUniqueMapClassificator, x::AbstractVector)::UniqueMapClassification

    if isnothing(mcr.uniq)
        uniq = unique(x)
    else
        uniq = mcr.uniq
    end
    k = length(uniq)
    group = zeros(Int, length(x))
    ngroup = zeros(Int, k)
    grouplabs = string.(uniq)

    for i in 1:k
        subset = x .== uniq[i]
        group[subset] .= i
        ngroup[i] = sum(subset)
    end

    return UniqueMapClassification(mcr, k, group, ngroup, grouplabs)
end

function maplabels(mc::UniqueMapClassification; counts::Bool = true)::Vector{String}

    ngroup = mc.ngroup
    grouplabs = mc.grouplabs
    labels = grouplabs

    if counts
        labels = labels .* " (" .* string.(ngroup) .* ")"
    end

    labels
end

levels(mc::UniqueMapClassification)::Vector{String} = return mc.grouplabs;


text:
# Autogenerated wrapper script for TetGen_jll for armv6l-linux-musleabihf
export libtet

JLLWrappers.@generate_wrapper_header("TetGen")
JLLWrappers.@declare_library_product(libtet, "libtet.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtet,
        "lib/libtet.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
newton(f,dfdx,x₁)

Use Newton's method to find a root of `f` starting from `x₁`, where
`dfdx` is the derivative of `f`. Returns a vector of root estimates.
"""
function newton(f,dfdx,x₁)
    # Operating parameters.
    funtol = 100*eps();  xtol = 100*eps();  maxiter = 40;

    x = [x₁]
    y = f(x₁)
    dx = Inf   # for initial pass below
    k = 1

    while (abs(dx) > xtol) && (abs(y) > funtol) && (k < maxiter)
        dydx = dfdx(x[k])
        dx = -y/dydx            # Newton step
        push!(x,x[k]+dx)        # append new estimate

        k += 1
        y = f(x[k])
    end

    if k==maxiter
        @warn "Maximum number of iterations reached."
    end

    return x
end

"""
secant(f,x₁,x₂)

Use the secant method to find a root of `f` starting from `x₁` and
`x₂`. Returns a vector of root estimates.
"""
function secant(f,x₁,x₂)
    # Operating parameters.
    funtol = 100*eps();  xtol = 100*eps();  maxiter = 40;

    x = [x₁,x₂]
    y₁ = f(x₁); y₂ = 100;
    dx = Inf   # for initial pass below
    k = 2

    while (abs(dx) > xtol) && (abs(y₂) > funtol) && (k < maxiter)
        y₂ = f(x[k])
        dx = -y₂ * (x[k]-x[k-1]) / (y₂-y₁)   # secant step
        push!(x,x[k]+dx)        # append new estimate

        k += 1
        y₁ = y₂    # current f-value becomes the old one next time
    end

    if k==maxiter
        @warn "Maximum number of iterations reached."
    end

    return x
end

"""
newtonsys(f,jac,x₁)

Use Newton's method to find a root of a system of equations,
starting from `x₁`. The functions `f` and `jac should return the
residual vector and the Jacobian matrix, respectively. Returns
history of root estimates as a vector of vectors.
"""
function newtonsys(f,jac,x₁)
    # Operating parameters.
    funtol = 1000*eps();  xtol = 1000*eps();  maxiter = 40;

    x = [float(x₁)]
    y,J = f(x₁),jac(x₁)
    dx = Inf   # for initial pass below
    k = 1

    while (norm(dx) > xtol) && (norm(y) > funtol) && (k < maxiter)
        dx = -(J\y)             # Newton step
        push!(x,x[k] + dx)    # append to history
        k += 1
        y,J = f(x[k]),jac(x[k])
    end

    if k==maxiter
        @warn "Maximum number of iterations reached."
    end

    return x
end

"""
fdjac(f,x0,y0)

Compute a finite-difference approximation of the Jacobian matrix for
`f` at `x0`, where `y0`=`f(x0)` is given.
"""
function fdjac(f,x0,y0)

δ = sqrt(eps())   # FD step size
m,n = length(y0),length(x0)
if n==1
    J = (f(x0+δ) - y0) / δ
else
    J = zeros(m,n)
    Iₙ = I(n)
    for j = 1:n
        J[:,j] = (f(x0 + δ*Iₙ[:,j]) - y0) / δ
    end
end

return J
end

"""
levenberg(f,x₁,tol)

Use Levenberg's quasi-Newton iteration to find a root of the system
`f`, starting from `x₁`, with `tol` as the stopping tolerance in
both step size and residual norm. Returns root estimates as a
matrix, one estimate per column.
"""
function levenberg(f,x₁,tol=1e-12)

# Operating parameters.
ftol = tol;  xtol = tol;  maxiter = 40;

x = zeros(length(x₁),maxiter)
x = [float(x₁)]
fₖ = f(x₁)
k = 1;  s = Inf;
Aₖ = fdjac(f,x₁,fₖ)   # start with FD Jacobian
jac_is_new = true

λ = 10;
while (norm(s) > xtol) && (norm(fₖ) > ftol) && (k < maxiter)
    # Compute the proposed step.
    B = Aₖ'*Aₖ + λ*I
    z = Aₖ'*fₖ
    s = -(B\z)

    xnew = x[k] + s
    fnew = f(xnew)

    # Do we accept the result?
    if norm(fnew) < norm(fₖ)    # accept
        y = fnew - fₖ
        push!(x,xnew)
        fₖ = fnew
        k += 1

        λ = λ/10   # get closer to Newton
        # Broyden update of the Jacobian.
        Aₖ = Aₖ + (y-Aₖ*s)*(s'/(s'*s))
        jac_is_new = false
    else                       # don't accept
        # Get closer to steepest descent.
        λ = 4λ
        # Re-initialize the Jacobian if it's out of date.
        if !jac_is_new
            Aₖ = fdjac(f,x[k],fₖ)
            jac_is_new = true
        end
    end
end

if (norm(fₖ) > 1e-3)
    @warn "Iteration did not find a root."
end

return x
end

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
using Gateway
using Test

@testset "GetGateway.jl" begin
	if Sys.islinux()
		result = Gateway.parse_linux_proc_net_route(open("/proc/net/route"))
	    @test Gateway.getgateway() == result
	end
end


text:
@testset "NAC" begin
    nac = NAC(2, 3)
    x = rand(2, 10)
    z = nac(x)

    @test size(z) == (3,10)

    ps = params(nac)
    @test length(ps) == 2

    loss() = Flux.mse(nac(x), reshape(repeat(x[1,:], inner=(3,1)),3,10))
    opt = RMSProp()
    train_data = Iterators.repeated((), 10000)
    Flux.train!(loss, ps, train_data, opt)

    @test all(isapprox.(weights(nac), [1.0 0.0; 1.0 0.0; 1.0 0.0], atol=1e-2))
end


text:
package java.lang;

public final class Integer extends Number implements Comparable {

	private int value;

	// Obviously ridiculous way to get the input to flow to the field
	public Integer(String integerConstructorArgString) throws NumberFormatException {
		if(integerConstructorArgString == null){
			throw new NumberFormatException("null");
		}
		value = integerConstructorArgString.length();
	}

	public Integer(int integerConstructorArgInt){
		value = integerConstructorArgInt;
	}

	public int intValue(){
		return value;
	}

	public int compareTo(Object another){
		Integer anotherInteger = (Integer)another;
		int thisVal = this.value;
		int anotherVal = anotherInteger.value;
		return (thisVal<anotherVal ? -1 : (thisVal==anotherVal ? 0 : 1));
    }

    public long longValue(){
    	return (long)value;
    }

    public float floatValue(){
    	return (float)value;
    }

    public double doubleValue(){
    	return (double)value;
    }

    public static Integer valueOf(int integerValue){
    	return new Integer(integerValue);
    }

    public static String toHexString(int i) {
	return "x: " + i;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
include("building_blocks.jl")
include("../src/physics.jl")

# Run experiments in parallel in order to compare behaviour
function run_experiment(contraptions, settings_list;
                        Δt::Float64 = 0.01, t_total::Float64 = 10)

    @assert(length(contraptions) == length(settings_list))
    n = length(contraptions)

    results = Array()
    for i in 1:n
        outcome = engine(contraptions[i], settings_list[i], Δt, t_total)
        push!(results, outcome)
    end

    return results
end

# Physics behaviour tests ####

# Contraptions should not accelerate in absence of external forces
function newton_first_law_test(contraption)
    results = run_experiment([contraption],
                             [null_settings])
    return false
end

@test newton_first_law_test(one_contraption)
@test newton_first_law_test(two_contraption)
@test newton_first_law_test(slow_contraption)

# Springs should pull points together when they're stretched
function stretch_test()
    results = run_experiment([stretch_contraption],
                             [null_settings],
                             t_total = 1.)

    return false
end

@test stretch_test()

# Springs should push points away when they're compressed
function compression_test()
    results = run_experiment([squash_contraption],
                             [null_settings],
                             t_total = 1.)
    return false
end

@test compression_test()

# Contraptions should fall when there's gravity
function falling_test(contraption)
    results = run_experiment([contraption],
                             [default_settings],
                             t_total = 1.)

    return false
end

@test falling_test(one_contraption)
@test falling_test(two_contraption)
@test falling_test(three_contraption)

# Contraptions should fall faster in more gravity
function gravity_falling_test(contraption)
    results = run_experiment([contraption for _ in 1:5],
                             [antigrav_settings,
                              nograv_settings,
                              lowgrav_settings,
                              default_settings,
                              higrav_settings],
                             t_total = 1.)

    return false
end

@test gravity_falling_test(one_contraption)
@test gravity_falling_test(three_contraption)

# Contraptions should fall slower with more drag
function drag_falling_test()
    results = run_experiment([contraption for _ in 1:4],
                             [nodrag_settings,
                              lowdrag_settings,
                              default_settings,
                              hidrag_settings],
                             t_total = 1.)
    return false
end

@test drag_falling_test(one_contraption)
@test drag_falling_test(three_contraption)

# Contraptions should oscillate at about the same speed regardless of initial compression or extension
function oscillation_direction_test()
    results = run_experiment([stretch_contraption, squash_contraption],
                             [null_settings, null_settings])

    return false
end

@test oscillation_direction_test()

# Contraptions should oscillate at about the same speed regardless of drag
function oscillation_drag_test()
    results = run_experiment([stretch_contraption for _ in 1:4],
                             [nodrag_settings,
                              lowdrag_settings,
                              default_settings,
                              hidrag_settings])
    return false
end

@test oscillation_drag_test()

# Contraptions should have dampened oscillations, with faster damping in more drag
function oscillation_damping_test()
    results = run_experiment([stretch_contraption for _ in 1:4],
                             [nodrag_settings,
                              lowdrag_settings,
                              default_settings,
                              hidrag_settings])
    return false
end

@test oscillation_damping_test()

# Contraptions should bounce
function bounce_test(contraption)
    results = run_experiment([contraption],
                             [default_settings])

    return false
end

@test bounce_test(one_contraption)
@test bounce_test(two_contraption)
@test bounce_test(three_contraption)
@test bounce_test(seven_contraption)

# Contraptions should bounce higher with more elasticity
function elasticity_bounce_test()
    results = run_experiment([contraption for _ in 1:5],
                             [inelastic_settings,
                              lowelastic_settings,
                              default_settings,
                              hielastic_settings,
                              elastic_settings])
    return false
end

@test elasticity_bounce_test(one_contraption)
@test elasticity_bounce_test(two_contraption)
@test elasticity_bounce_test(three_contraption)
@test elasticity_bounce_test(seven_contraption)
@test elasticity_bounce_test(pogo_contraption)

# Energy should be conserved in the absence of drag or inelasticity
function energy_conservation_test(contraption)
    results = run_experiment([contraption],
                             [noloss_settings])

    return false
end

@test energy_conservation_test(one_contraption)
@test energy_conservation_test(seven_contraption)
@test energy_conservation_test(fast_contraption)
@test energy_conservation_test(strong_contraption)

# Contraptions should eventually come to rest
function come_to_rest_test(contraption)
    results = run_experiment([contraption],
                             [default_settings])
    return false
end

@test come_to_rest_test(one_contraption)
@test come_to_rest_test(two_contraption)
@test come_to_rest_test(three_contraption)

# Contraptions should spin after bouncing iff they're off-centre
function spin_test()
    results = run_experiment([seven_contraption, rotated],
                             [default_settings, default_settings])

    return false
end

@test spin_test()

# Larger contraptions should hit the ground sooner
function large_collision_test()
    results = run_experiment([small_contraption,
                              seven_contraption,
                              large_contraption],
                             [default_settings,
                              default_settings,
                              default_settings])
    return false
end

@test large_collision_test()

# Contraptions should take longer to hit the wall in larger environments
function large_env_test()
    results = run_experiment([seven_contraption,
                              seven_contraption,
                              seven_contraption],
                             [small_settings,
                              default_settings,
                              large_settings])

    return false
end

@test large_env_test()

# Larger contraptions should rest at a higher point
function large_rest_test()
    results = run_experiment([small_contraption,
                              seven_contraption,
                              large_contraption],
                             [default_settings,
                              default_settings,
                              default_settings])

    return false
end

@test large_rest_test()

# Light contraptions should rest at a higher point
function heavy_rest_test()
    results = run_experiment([light_contraption,
                              seven_contraption,
                              heavy_contraption],
                             [default_settings,
                              default_settings,
                              default_settings])

    return false
end

@test heavy_rest_test()

# Strong contraptions should rest at a higher point
function strong_rest_test()
    results = run_experiment([weak_contraption,
                              seven_contraption,
                              strong_contraption],
                             [default_settings,
                              default_settings,
                              default_settings])

    return false
end

@test heavy_rest_test()

# Fast contraptions should be approximately uniformly distributed over the bounding box at random times
function fast_distribution_test()
    results = run_experiment([fast_contraption],
                             [null_settings])

    return false
end

@test fast_distribution_test()


text:
################################################################################
#
#   fasta : Accelerated Forward Backward splitting with
#            adaptative step size and backtracking
#
################################################################################

const NLookback=10 # previous steps lookback to decide backtracking

function fasta!{T<:FloatingPoint}(
    rmodel::SupervisedRiskModel,     # model definition
    θ::StridedArray{T},              # initial parameter guess
    x::StridedArray{T},              # inputs
    y::StridedArray{T};              # responses
    reg::Regularizer = ZeroReg(),    # regularizer
    λ::T             = one(T),       # initial step size
    cbinterval::Int  = 0,            # callback interval
    callback         = fasta_trace,  # callback function
    maxsteps::Int    = 100,          # max number of steps
    maxtime::Int     = 0         )   # max calculation time (0 = no limit)

    g              = similar(θ) # holds gradient
    prevg          = similar(θ) # gradient of previous step
    deltag         = similar(θ) # g - predg

    prevθ          = similar(θ) # θ of previous step
    deltaθ         = similar(θ) # θ - prevθ

    dualθ          = copy(θ)    # dual θ
    candidatedualθ = similar(θ) # candidate dual θ
    deltadualθ     = similar(θ)

    vhist     = fill(-Inf, NLookback) # past values

    function fbstep!(θ1, g, θ0, λ)
        copy!(θ1, θ0)
        axpy!(-λ, g, θ1)
        prox!(reg, θ1, θ1, λ)
    end

    ### Loop initialization
    starttime = time()
    α, μ      = 1., 0.      # acceleration parameters
    v, _      = value_and_addgrad!(rmodel, 0., g, 1., dualθ, x, y)
    vhist[1]  = v

    cbinterval > 0 && callback(0, θ, v, value(reg, θ), λ, α, 0)
    for t in 1:maxsteps
        (maxtime != 0) && (time() - starttime > maxtime) && break

        # acceleration
        if t >= 2
            if dot(g, deltadualθ) >= 0.
                α = 1.
                μ = 0.
            else
                candα = (1. + sqrt(1+4α*α))/2.
                μ = (α - 1.) / candα
                α = candα
            end
        end

        # adaptative step size
        if t >= 2
            copy!(deltag, g) ; axpy!(-1., prevg, deltag)
            tmp = dot(deltadualθ, deltag)
            λs  = dot(deltadualθ, deltadualθ) / tmp
            λm  = tmp / dot(deltag, deltag)
            if λs > 0. && λm > 0.
                λ = (λm / λs > 0.5) ? λm : λs - 0.5λm
            end
        end

        # update, with backtracking
        btsteps = 0
        maxv = maximum(vhist)
        copy!(prevθ, θ)
        while btsteps < 10
            fbstep!(θ, g, dualθ, λ)
            copy!(deltaθ, θ) ; axpy!(-1., prevθ, deltaθ)
            copy!(candidatedualθ, θ)
            axpy!(μ, deltaθ, candidatedualθ)

            candidatev = value(rmodel, candidatedualθ, x, y)
            deltadualθ = candidatedualθ - dualθ # TODO : preallocate
            if candidatev < maxv + dot(deltadualθ, g) + sumabs2(deltadualθ)/2λ
                break
            end
            λ /= 2.
            btsteps += 1
        end
        copy!(dualθ     , candidatedualθ)

        # update variables for next step
        copy!(prevg     , g)
        v, _ = value_and_addgrad!(rmodel, 0., g, 1., dualθ, x, y)
        unshift!(vhist, v)
        pop!(vhist)

        # callback
        if cbinterval > 0 && t % cbinterval == 0
            callback(t, θ, v, value(reg, θ), λ, α, btsteps)
        end
    end

    θ
end

# Non mutating version
function fasta{T<:FloatingPoint}(
    rmodel::SupervisedRiskModel,     # model definition
    θ::StridedArray{T},              # initial parameter guess
    x::StridedArray{T},              # inputs
    y::StridedArray{T};              # responses
    vargs... )

    fasta!(rmodel, copy(θ), x, y; vargs...)
end


text:
using NorvigsTrieSpellchecker
using DataStructures
using Test

@testset "NorvigsTrieSpellchecker.jl" begin
    t = Trie{Tuple{String, Float64}}()
    t["abc"] = ("abc", 1.0)

    # delete
    @test correction(t, "abdc") == "abc"
    # insert
    @test correction(t, "ac") == "abc"
    # replace
    @test correction(t, "adc") == "abc"
    # transpose
    @test correction(t, "bac") =="abc"
    # replace 2
    @test correction(t, "dbe") == "abc"
    # insert 2
    @test correction(t, "a") == "abc"
    # transpose + delete
    @test correction(t, "ba") =="abc"
    # known
    @test correction(t, "abc") == "abc"
    # unknown
    @test correction(t, "quintessential") == "quintessential"
    # end insert
    @test correction(t, "ab") == "abc"
    # end replace
    @test correction(t, "abd") == "abc"
    # end transpose
    @test correction(t, "acb") == "abc"
    # end delete
    @test correction(t, "abcd") == "abc"
    # begin insert
    @test correction(t, "bc") == "abc"

    t = Trie{Tuple{String, Float64}}()
    t["abcd"] = ("abcd", 0.5)
    t["abc"] = ("abc", 1.0)
    @test correction(t, "abce") == "abc"

    t = Trie{Tuple{String, Float64}}()
    t["inst"] = ("inst", 0.5)
    t["nut"] = ("nut", 1.0)
    @test correction(t, "inut") == "nut"

    t = Trie{Tuple{String, Float64}}()
    t["trees"] = ("trees", 0.5)
    t["these"] = ("these", 1.0)
    @test correction(t, "thees") == "these"
end


text:
typealias ReComp Union{Real,Complex}

immutable Dual{T<:ReComp} <: Number
    value::T
    epsilon::T
end
Dual{S<:ReComp,T<:ReComp}(x::S, y::T) = Dual(promote(x,y)...)
Dual(x::ReComp) = Dual(x, zero(x))

const ɛ = Dual(false, true)
const imɛ = Dual(Complex(false, false), Complex(false, true))

typealias Dual128 Dual{Float64}
typealias Dual64  Dual{Float32}
typealias Dual32  Dual{Float16}
typealias DualComplex256 Dual{Complex128}
typealias DualComplex128 Dual{Complex64}
typealias DualComplex64  Dual{Complex32}

convert{T<:ReComp}(::Type{Dual{T}}, z::Dual{T}) = z
convert{T<:ReComp}(::Type{Dual{T}}, z::Dual) = Dual{T}(convert(T, value(z)), convert(T, epsilon(z)))
convert{T<:ReComp}(::Type{Dual{T}}, x::Number) = Dual{T}(convert(T, x), convert(T, 0))
convert{T<:ReComp}(::Type{T}, z::Dual) = (epsilon(z)==0 ? convert(T, value(z)) : throw(InexactError()))

promote_rule{T<:ReComp, S<:ReComp}(::Type{Dual{T}}, ::Type{Dual{S}}) = Dual{promote_type(T, S)}
promote_rule{T<:ReComp, S<:ReComp}(::Type{Dual{T}}, ::Type{S}) = Dual{promote_type(T, S)}
promote_rule{T<:ReComp}(::Type{Dual{T}}, ::Type{T}) = Dual{T}

widen{T}(::Type{Dual{T}}) = Dual{widen(T)}

value(z::Dual) = z.value
epsilon(z::Dual) = z.epsilon

dual(x::ReComp, y::ReComp) = Dual(x, y)
dual(x::ReComp) = Dual(x)
dual(z::Dual) = z

Compat.@dep_vectorize_1arg ReComp dual
Compat.@dep_vectorize_2arg ReComp dual
Compat.@dep_vectorize_1arg Dual dual
Compat.@dep_vectorize_1arg Dual value
Compat.@dep_vectorize_1arg Dual epsilon

const realpart = value
const dualpart = epsilon

isnan(z::Dual) = isnan(value(z))
isinf(z::Dual) = isinf(value(z))
isfinite(z::Dual) = isfinite(value(z))
isdual(x::Dual) = true
isdual(x::Number) = false
eps(z::Dual) = eps(value(z))
eps{T}(::Type{Dual{T}}) = eps(T)

function dual_show{T<:Real}(io::IO, z::Dual{T}, compact::Bool)
    x, y = value(z), epsilon(z)
    if isnan(x) || isfinite(y)
        compact ? showcompact(io,x) : show(io,x)
        if signbit(y)==1 && !isnan(y)
            y = -y
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
        compact ? showcompact(io, y) : show(io, y)
        printtimes(io, y)
        print(io, "ɛ")
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function dual_show{T<:Complex}(io::IO, z::Dual{T}, compact::Bool)
    x, y = value(z), epsilon(z)
    xr, xi = reim(x)
    yr, yi = reim(y)
    if isnan(x) || isfinite(y)
        compact ? showcompact(io,x) : show(io,x)
        if signbit(yr)==1 && !isnan(y)
            yr = -yr
            print(io, " - ")
        else
            print(io, " + ")
        end
        if compact
            if signbit(yi)==1 && !isnan(y)
                yi = -yi
                showcompact(io, yr)
                printtimes(io, yr)
                print(io, "ɛ-")
                showcompact(io, yi)
            else
                showcompact(io, yr)
                print(io, "ɛ+")
                showcompact(io, yi)
            end
        else
            if signbit(yi)==1 && !isnan(y)
                yi = -yi
                show(io, yr)
                printtimes(io, yr)
                print(io, "ɛ - ")
                show(io, yi)
            else
                show(io, yr)
                print(io, "ɛ + ")
                show(io, yi)
            end
        end
        printtimes(io, yi)
        print(io, "imɛ")
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function dual_show{T<:Bool}(io::IO, z::Dual{T}, compact::Bool)
    x, y = value(z), epsilon(z)
    if !value(z) && epsilon(z)
        print(io, "ɛ")
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function dual_show{T<:Bool}(io::IO, z::Dual{Complex{T}}, compact::Bool)
    x, y = value(z), epsilon(z)
    xr, xi = reim(x)
    yr, yi = reim(y)
    if !xr
        if xi*!yr*!yi
            print(io, "im")
        elseif !xi*yr*!yi
            print(io, "ɛ")
        elseif !xi*!yr*yi
            print(io, "imɛ")
        end
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function printtimes(io::IO, x::Real)
    if !(isa(x,Integer) || isa(x,Rational) ||
         isa(x,AbstractFloat) && isfinite(x))
        print(io, "*")
    end
end

show(io::IO, z::Dual) = dual_show(io, z, false)
showcompact(io::IO, z::Dual) = dual_show(io, z, true)

function read{T<:ReComp}(s::IO, ::Type{Dual{T}})
    x = read(s, T)
    y = read(s, T)
    Dual{T}(x, y)
end
function write(s::IO, z::Dual)
    write(s, value(z))
    write(s, epsilon(z))
end

## Generic functions of dual numbers ##

convert(::Type{Dual}, z::Dual) = z
convert(::Type{Dual}, x::Number) = Dual(x)

==(z::Dual, w::Dual) = value(z) == value(w)
==(z::Dual, x::Number) = value(z) == x
==(x::Number, z::Dual) = value(z) == x

isequal(z::Dual, w::Dual) = isequal(value(z),value(w)) && isequal(epsilon(z), epsilon(w))
isequal(z::Dual, x::Number) = isequal(value(z), x) && isequal(epsilon(z), zero(x))
isequal(x::Number, z::Dual) = isequal(z, x)

isless{T<:Real}(z::Dual{T},w::Dual{T}) = value(z) < value(w)
isless{T<:Real}(z::Real,w::Dual{T}) = z < value(w)
isless{T<:Real}(z::Dual{T},w::Real) = value(z) < w

hash(z::Dual) = (x = hash(value(z)); epsilon(z)==0 ? x : bitmix(x,hash(epsilon(z))))

float{T<:AbstractFloat}(z::Union{Dual{T},Dual{Complex{T}}})=z
complex{T<:Real}(z::Dual{Complex{T}})=z

floor{T<:Real}(::Type{T}, z::Dual) = floor(T, value(z))
ceil{ T<:Real}(::Type{T}, z::Dual) = ceil( T, value(z))
trunc{T<:Real}(::Type{T}, z::Dual) = trunc(T, value(z))
round{T<:Real}(::Type{T}, z::Dual) = round(T, value(z))

for op in (:real,:imag,:conj,:float,:complex)
    @eval begin
        $op(z::Dual) = Dual($op(value(z)),$op(epsilon(z)))
    end
end

abs(z::Dual) = sqrt(abs2(z))
abs2(z::Dual) = real(conj(z)*z)

real{T<:Real}(z::Dual{T}) = z
abs{T<:Real}(z::Dual{T}) = z ≥ 0 ? z : -z

angle{T<:Real}(z::Dual{T}) = z ≥ 0 ? zero(z) : one(z)*π
angle{T<:Real}(z::Dual{Complex{T}}) = z == 0 ? (imag(epsilon(z)) == 0 ? Dual(zero(T),zero(T)) : Dual(zero(T),convert(T, Inf))) : real(log(sign(z))/im)

# algebraic definitions
conjdual(z::Dual) = Dual(value(z),-epsilon(z))
absdual(z::Dual) = abs(value(z))
abs2dual(z::Dual) = abs2(value(z))

# algebra

+(z::Dual, w::Dual) = Dual(value(z)+value(w), epsilon(z)+epsilon(w))
+(z::Number, w::Dual) = Dual(z+value(w), epsilon(w))
+(z::Dual, w::Number) = Dual(value(z)+w, epsilon(z))

-(z::Dual) = Dual(-value(z), -epsilon(z))
-(z::Dual, w::Dual) = Dual(value(z)-value(w), epsilon(z)-epsilon(w))
-(z::Number, w::Dual) = Dual(z-value(w), -epsilon(w))
-(z::Dual, w::Number) = Dual(value(z)-w, epsilon(z))

# avoid ambiguous definition with Bool*Number
*(x::Bool, z::Dual) = ifelse(x, z, ifelse(signbit(real(value(z)))==0, zero(z), -zero(z)))
*(x::Dual, z::Bool) = z*x

*(z::Dual, w::Dual) = Dual(value(z)*value(w), epsilon(z)*value(w)+value(z)*epsilon(w))
*(x::Number, z::Dual) = Dual(x*value(z), x*epsilon(z))
*(z::Dual, x::Number) = Dual(x*value(z), x*epsilon(z))

/(z::Dual, w::Dual) = Dual(value(z)/value(w), (epsilon(z)*value(w)-value(z)*epsilon(w))/(value(w)*value(w)))
/(z::Number, w::Dual) = Dual(z/value(w), -z*epsilon(w)/value(w)^2)
/(z::Dual, x::Number) = Dual(value(z)/x, epsilon(z)/x)

for f in [:^, :(NaNMath.pow)]
    @eval function ($f)(z::Dual, w::Dual)
        if epsilon(w) == 0.0
            return $f(z,value(w))
        end
        val = $f(value(z),value(w))

        du =
        epsilon(z)*value(w)*(($f)(value(z),value(w)-1))+epsilon(w)*($f)(value(z),value(w))*log(value(z))

        Dual(val, du)
    end
end

mod(z::Dual, n::Number) = Dual(mod(value(z), n), epsilon(z))

# these two definitions are needed to fix ambiguity warnings
^(z::Dual, n::Integer) = Dual(value(z)^n, epsilon(z)*n*value(z)^(n-1))
^(z::Dual, n::Rational) = Dual(value(z)^n, epsilon(z)*n*value(z)^(n-1))

^(z::Dual, n::Number) = Dual(value(z)^n, epsilon(z)*n*value(z)^(n-1))
NaNMath.pow(z::Dual, n::Number) = Dual(NaNMath.pow(value(z),n), epsilon(z)*n*NaNMath.pow(value(z),n-1))
NaNMath.pow(z::Number, w::Dual) = Dual(NaNMath.pow(z,value(w)), epsilon(w)*NaNMath.pow(z,value(w))*log(z))

# force use of NaNMath functions in derivative calculations
function to_nanmath(x::Expr)
    if x.head == :call
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(x.args[1]))
        return Expr(:call,funsym,[to_nanmath(z) for z in x.args[2:end]]...)
    else
        return Expr(:call,[to_nanmath(z) for z in x.args]...)
    end
end
to_nanmath(x) = x

for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
    funsym == :exp && continue
    funsym == :abs2 && continue
    @eval function $(funsym)(z::Dual)
        x = value(z)
        xp = epsilon(z)
        Dual($(funsym)(x),xp*$exp)
    end
    # extend corresponding NaNMath methods
    if funsym in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10,
          :lgamma, :log1p)
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(funsym))
        @eval function $(funsym)(z::Dual)
            x = value(z)
            xp = epsilon(z)
            Dual($(funsym)(x),xp*$(to_nanmath(exp)))
        end
    end
end

# only need to compute exp/cis once
exp(z::Dual) = (expval = exp(value(z)); Dual(expval, epsilon(z)*expval))
cis(z::Dual) = (cisval = cis(value(z)); Dual(cisval, im*epsilon(z)*cisval))

## TODO: should be generated in Calculus
sinpi(z::Dual) = Dual(sinpi(value(z)),epsilon(z)*cospi(value(z))*π)
cospi(z::Dual) = Dual(cospi(value(z)),-epsilon(z)*sinpi(value(z))*π)