text:
"""
    point foot quadruped with linear springs and dampers on each foot
    q = (p, r, f1, f2, f3, f4)
        p - body position
        r - body orientation (modified Rodriques parameters)
        f1 - foot 1 position
        f2 - foot 2 position
        f3 - foot 3 position
        f4 - foot 4 position
"""
mutable struct PointFootQuadruped{T} <: Model{T}
    # dimensions
	nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters
    nc::Int # contact points

	# parameters
	body_height::T
	foot_x::T
	foot_y::T

	# parameters
    mass_body::T
    inertia_body::Matrix{T}
    mass_foot::T

	# joint friction
	friction_joint::Vector{T}
	spring_stiffness_joint::Vector{T}

	# environment
    friction_foot_world::Vector{T}
    gravity::T
end

# function L_mult(x)
#     [x[1] -transpose(x[2:4]);
#      x[2:4] x[1] * I(3) + skew(x[2:4])]
# end
#
# # right quaternion multiply as matrix
# function R_mult(x)
#     [x[1] -transpose(x[2:4]); x[2:4] x[1] * I(3) - skew(x[2:4])]
# end

# rotation matrix
function quat_rot_mat(q)
    H = [zeros(1, 3); I(3)]
    transpose(H) * L_mult(q) * transpose(R_mult(q)) * H
end

function quat_from_mrp(p)
    """Quaternion (scalar first) from MRP"""
    return (1.0 / (1.0 + dot(p, p))) * [(1 - dot(p, p)); 2.0 * p]
end

function mrp_rot_mat(x)
    quat_rot_mat(quat_from_mrp(x))
end

# Kinematics
function kinematics(model::PointFootQuadruped, q)
	q[6 .+ (1:12)]
end

lagrangian(model::PointFootQuadruped, q, q̇) = 0.0

function mass_matrix(model::PointFootQuadruped, q)
    cat(
        model.mass_body * Diagonal(ones(3)),     # body position
        model.inertia_body,                      # body orienation
        model.mass_foot * Diagonal(ones(3 * 4)), # feet position
        dims=(1, 2)
        )
end

function dynamics_bias(model::PointFootQuadruped, q, q̇)
	offsets = [
		[+model.foot_x, +model.foot_y, -model.body_height],
		[+model.foot_x, -model.foot_y, -model.body_height],
		[-model.foot_x, +model.foot_y, -model.body_height],
		[-model.foot_x, -model.foot_y, -model.body_height],
		]

	spring = model.spring_stiffness_joint
    [
        model.mass_body * [0,0,model.gravity] + spring .* (+sum(offsets) + 4 * q[1:3] - q[7:9] - q[10:12] - q[13:15] - q[16:18]);     # body position
        skew(q̇[4:6]) * model.inertia_body * q̇[4:6]; # body orienation
        model.mass_foot * [0,0,model.gravity] + spring .* (-offsets[1] + q[7:9] - q[1:3]);
        model.mass_foot * [0,0,model.gravity] + spring .* (-offsets[2] + q[10:12] - q[1:3]);
        model.mass_foot * [0,0,model.gravity] + spring .* (-offsets[3] + q[13:15] - q[1:3]);
        model.mass_foot * [0,0,model.gravity] + spring .* (-offsets[4] + q[16:18] - q[1:3]);
    ]
end

function signed_distance(model::PointFootQuadruped, q)

    position_foot1 = q[6 .+ (1:3)]
    position_foot2 = q[9 .+ (1:3)]
    position_foot3 = q[12 .+ (1:3)]
	position_foot4 = q[15 .+ (1:3)]

	# return [q[3]; position_foot1[3]; position_foot2[3]; position_foot3[3]; position_foot4[3]]
	return [position_foot1[3]; position_foot2[3]; position_foot3[3]; position_foot4[3]]
end

function input_jacobian(model::PointFootQuadruped, q)
    position_body = q[1:3]
    orientation_body = q[3 .+ (1:3)]
	R = mrp_rot_mat(orientation_body)
	# R = euler_rotation_matrix(orientation_body)

	# kinematics in world frame
	r1 = q[6 .+ (1:3)] - position_body
	r2 = q[9 .+ (1:3)] - position_body
	r3 = q[12 .+ (1:3)] - position_body
	r4 = q[15 .+ (1:3)] - position_body

	z3 = zeros(3, 3)

	transpose([
        I(3) z3   I(3) I(3) I(3) I(3);
		z3   I(3) transpose(R) * skew(r1) transpose(R) * skew(r2) transpose(R) * skew(r3) transpose(R) * skew(r4);
        z3   z3   -I(3)    z3    z3   z3;
        z3   z3   z3    -I(3)    z3   z3;
        z3   z3   z3       z3 -I(3)   z3;
        z3   z3   z3       z3    z3 -I(3)
    ])
end

function contact_jacobian(model::PointFootQuadruped, q)
    z3 = zeros(3, 3)

    [
        z3   z3 I(3)   z3   z3   z3;
		z3   z3   z3 I(3)   z3   z3;
        z3   z3   z3   z3 I(3)   z3;
        z3   z3   z3   z3   z3 I(3);
    ]
end

function nominal_configuration(model::PointFootQuadruped)
    [
        0.0; 0.0; model.body_height;
        0.0; 0.0; 0.0;
        model.foot_x; model.foot_y; 0.0;
        model.foot_x;-model.foot_y; 0.0;
       -model.foot_x; model.foot_y; 0.0;
       -model.foot_x;-model.foot_y; 0.0;
    ]
end

function nominal_state(model::PointFootQuadruped)
	[nominal_configuration(model); zeros(model.nq)]
end

# friction coefficients
function friction_coefficients(model::PointFootQuadruped)
    model.friction_foot_world
end

# dimensions
nq = 3 + 3 + 3 * 4       # generalized coordinates
nu = 6 + 3 * 4           # controls
nw = 3                   # parameters
nc = 4                   # contact points

# inertial properties
scale = 10.0
mass_body = 13.5 / scale
i_xx = 0.0178533 / scale
i_xy = 0.0 / scale
i_xz = 0.0 / scale
i_yz = 0.0 / scale
i_yy = 0.0377999 / scale
i_zz = 0.0456542 / scale
inertia_body = Array(Diagonal([i_xx, i_yy, i_zz]))
mass_foot = 0.2

# parameters
body_height = 0.3
foot_x = 0.17
foot_y = 0.15

friction_joint = 0.15 * ones(3)
spring_stiffness_joint = [5,5,15.0]
friction_foot_world = 0.3*ones(nc)     # coefficient of friction
gravity = 9.81                   # gravity

point_foot_quadruped = PointFootQuadruped(nq, nu, nw, nc,
				body_height,
				foot_x,
				foot_y,
				mass_body,
				inertia_body,
				mass_foot,
				friction_joint,
				spring_stiffness_joint,
				friction_foot_world,
				gravity,
                )

point_foot_quadruped_contact_kinematics = [
    q -> q[6  .+ (1:3)],
    q -> q[9  .+ (1:3)],
    q -> q[12 .+ (1:3)],
    q -> q[15 .+ (1:3)],
]

point_foot_quadruped_contact_kinematics_jacobians = [
    q -> [zeros(3, 6) I(3) zeros(3, 9)],
    q -> [zeros(3, 9) I(3) zeros(3, 6)],
    q -> [zeros(3, 12) I(3) zeros(3, 3)],
    q -> [zeros(3, 15) I(3)],
]

floating_base_dim(::PointFootQuadruped) = 6
name(::PointFootQuadruped) = :point_foot_quadruped


text:
function solve()
    n = parse(Int, readline())
    a = [parse(Int, x) for x in split(readline())]
    b = [parse(Int, x) for x in split(readline())]

    a'b == 0 ? "Yes" : "No"
end

println(solve())

ext: jl

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
#=

  Global constraint inverse in Julia ConstraintSolver.jl

  From Global Constraint Catalog
  http://www.emn.fr/z-info/sdemasse/gccat/Cinverse.html
  """
  inverse(NODES)

  Enforce each vertex of a digraph to have exactly one predecessor and
  one successor. In addition the following two statements are equivalent:
    - The successor of the ith node is the jth node.
    - The predecessor of the jth node is the ith node.
  """

  Note: This is implemented as "self-assignment", i.e. that
     inverse(X) -> assignment(X,X)
  This means that for each element X[I] either
     - X[I] = I
     or
     - X[I] = J <=> X[J] = I

  Example N=4:

    [1,2,3,4]
    [1,2,4,3]    (X[3]=4, X[4]=3)
    [1,3,2,4]    (X[2]=3, X[3]=2)
    [1,4,3,2]
    [2,1,3,4]
    [2,1,4,3]    (X[1]=2,X[2]=1, X[3]=4,X[4]=3)
    [3,2,1,4]
    [3,4,1,2]
    [4,2,3,1]
    [4,3,2,1]

   Number of solutions for N=1..11:
     [1,2,4,10,26,76,232,764,2620,9496,35696]
   OEIS:
   http://oeis.org/search?q=1%2C2%2C4%2C10%2C26%2C76%2C232%2C764%2C2620%2C9496%2C35696&language=english&go=Search

   http://oeis.org/A000085
   """
   Number of self-inverse permutations on n letters, also known as
   involutions; number of Young tableaux with n cells.
   """

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

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function inverse_test(n,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    @variable(model, 1 <= x[1:n] <= n, Int)
    @constraint(model, x in CS.AllDifferent())

    inverse(model,x)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$xx")

            end
        end
        # com = JuMP.backend(model).optimizer.model.inner
        # println("info:$(com.info)")
        # println("start_time:$(com.start_time)")
        # println("solve_time:$(com.solve_time)")
        println()
    else
        println("status:$status")
    end

    return status
end

function test_inverse_range(from,to)
    print_solutions=true
    all_solutions=true
    for n in from:to
        println("n:$n")
        if n > 4
            print_solutions = false
        end
        if n > 9
            all_solutions=false
            print_solutions=true
        end
        @time status = inverse_test(n,print_solutions,all_solutions)
        if status != MOI.OPTIMAL
            break
        end
    end
end

test_inverse_range(1,21)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -20.345339284625858 <= q <= 16.479136897808104)

                     add_NL_constraint(m, :(sigmoid(0.24719801624822457 + 0.29438950671798914*sigmoid(0.7291821855933653 + 0.701331929085816*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.7293079760660532*sigmoid(0.05096633991424904 + 0.31501139336769945*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8601808654748551*sigmoid(0.18255055897839467 + -0.7759686413687503*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.9720568561360068*sigmoid(0.5199334546165302 + 0.4081309459110076*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + sigmoid(-0.49402133773229506 + 0.2465550492028581*sigmoid(0.7291821855933653 + 0.701331929085816*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.17626257051164496*sigmoid(0.05096633991424904 + 0.31501139336769945*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8688281780619751*sigmoid(0.18255055897839467 + -0.7759686413687503*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.3904178182819744*sigmoid(0.5199334546165302 + 0.4081309459110076*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + sigmoid(-0.6759660684871824 + 0.9311431633625067*sigmoid(0.7291821855933653 + 0.701331929085816*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.2015839225926066*sigmoid(0.05096633991424904 + 0.31501139336769945*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.01676614484829253*sigmoid(0.18255055897839467 + -0.7759686413687503*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.07861065087286123*sigmoid(0.5199334546165302 + 0.4081309459110076*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + sigmoid(-0.08101830092609674 + 0.5708385401734484*sigmoid(0.7291821855933653 + 0.701331929085816*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.9629021826429458*sigmoid(0.05096633991424904 + 0.31501139336769945*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.4780289946412233*sigmoid(0.18255055897839467 + -0.7759686413687503*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.25224290978348307*sigmoid(0.5199334546165302 + 0.4081309459110076*sigmoid(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*sigmoid(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*sigmoid(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*sigmoid(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#Code written by Skipper af Rosenborg and Esben Bager
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
println("Leeeeroooy Jenkins")

industry = "NoDur"
folder = "240"

# Load all 3 AALR data
nRows = 845
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/AALRTest/"*string(folder)*"-1/"
cd(path)

AALRMainData = zeros(3, nRows, 1472)
for i = 1:nRows
    f = open(string(i)*"AALRBestK_240.csv")
    readline(f)
    s = readline(f)
    AALRMainData[1, i, :] = parse.(Float64, split(s, ","))

    s = readline(f)
    AALRMainData[2, i, :] = parse.(Float64, split(s, ","))

    s = readline(f)
    AALRMainData[3, i, :] = parse.(Float64, split(s, ","))
end
AALRMainData = AALRMainData[:,:,9:end]

# Get column names
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
industry = "NoDur"
mainData = loadIndexDataLOGReturn(industry, path)
colNames = names(mainData)
featureNames = expandColNamesTimeToString(colNames[1:end-3])
s = split(expandedColNamesToString(featureNames,true,true),",")
s = Array{String}(s)

# Load all best RMSE bsolved
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/"
path = path*"Results/IndexData/LassoTest/NoDur/"
summary = CSV.read(path*"Summary "*folder*".csv",nullable = false)
minIndex = summary[end-1,end]

#Load best beta matrix
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/"
path = path*"Results/IndexData/LassoTest/NoDur/"*string(folder)*"-1/"
cd(path)
path
dType = [Float64]
for i = 1:1463
	dType = vcat(dType,Float64)
end
LassoMainData = CSV.read("240_VIX_MacroTimeExpTA_bMatrix240_"*string(minIndex)".csv", header=s,
    delim = ',', nullable=false, types = dType)

#Remove all columns that are never active for both AALR and LASSO
LassoMainDataActive = Array{Float64}(size(LassoMainData)[1],0)
LassoHeaderActive = Array{String}(0)
LassoParamUsage = Array{Float64}(0)
for i=1:size(LassoMainData)[2]
	paramUsage = countnz(LassoMainData[:,i])
	if paramUsage > 0
		LassoParamUsage = vcat(LassoParamUsage, paramUsage)
		LassoMainDataActive = hcat(LassoMainDataActive,LassoMainData[:,i])
		LassoHeaderActive = vcat(LassoHeaderActive,s[i])
	end
end
#Number of active parameters
LassoTotalNumActiveParam = size(LassoMainDataActive)[2]
println("Number of active parameters is ", LassoTotalNumActiveParam)
#Average parameter usage
LassoAvgParamUsage = mean(LassoParamUsage)/size(LassoMainData)[2]
println("Most used parameter ", maximum(LassoParamUsage)/size(LassoMainData)[1])
LassoParamUsageCopy = LassoParamUsage
test = zeros(20,2)
#Top 20 most used prediction terms
for i = 1:20
    maxind = findmax(LassoParamUsageCopy)[2]
    LassoParamUsageCopy[maxind] = 0
    test[i,1] = maxind
    println(i, " most used prediction term is index ", maxind)
end
println("An active parameter in LASSO is used ", mean(LassoParamUsage)/size(LassoMainData)[2]*100,"% of the time")

AALRMainData[1,:,:]
AALRMainDataActive = Array{Float64}(size(AALRMainData)[2],0)
AALRHeaderActive = Array{String}(0)
AALRParamUsage = Array{Float64}(0)
for i=1:size(AALRMainData[1,:,:])[2]
    paramUsage = length(find(x -> x>2e-6 || x<-2e-6, AALRMainData[1,:,i]))
	if paramUsage > 0
		AALRParamUsage = vcat(AALRParamUsage, paramUsage)
		AALRMainDataActive = hcat(AALRMainDataActive,AALRMainData[1,:,i])
		AALRHeaderActive = vcat(AALRHeaderActive,s[i])
	end
end
AALRTotalNumActiveParam = size(AALRMainDataActive)[2]
println("Number of active parameters is ",AALRTotalNumActiveParam)
AALRAvgParamUsage = mean(AALRParamUsage)/size(AALRMainData[1,:,:])[2]
println("Most used parameter ", maximum(AALRParamUsage)/size(LassoMainData)[1])
findmax(AALRParamUsage)
AALRParamUsageCopy = AALRParamUsage
for i = 1:20
    maxind = findmax(AALRParamUsageCopy)[2]
    AALRParamUsageCopy[maxind] = 0
    test[i,2] = maxind
    println(i, " most used prediction term is index ", maxind)
end
println("An active parameter in AALR is used ", mean(AALRParamUsage)/size(AALRMainData[1,:,:])[2]*100,"% of the time")

for i = 1:20
    println("At i = ", i, " we found ", find(x -> x == test[i,1], test[:,2]))
end

#Count average number of parameters in AALR and LASSO
LassoMainDataArr = Array{Float64}(LassoMainData)
LassoKArr = Array{Int64}(size(LassoMainData)[1])
for i=1:size(LassoMainData)[1]
	LassoKArr[i] = countnz(LassoMainDataArr[i,:])
end

lassoAvgK = mean(LassoKArr)
lassoStd = std(LassoKArr)

AALRMainDataArr = Array{Float64}(AALRMainData[1,:,:])
AALRKArr = Array{Int64}(size(AALRMainDataArr)[1])
for i=1:size(AALRMainDataArr)[1]
    AALRKArr[i] = size(find(x -> x >= 1e-6 || x <=-1e-6, AALRMainDataArr[i,:]))[1]
end
AALRAvgK = mean(AALRKArr)
AALRStd = std(AALRKArr)

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/LassoAALRComparison.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
module MatrixLMnet

using MatrixLM
import MatrixLM.calc_preds, MatrixLM.calc_preds!,
    MatrixLM.calc_resid, MatrixLM.calc_resid!

using DataFrames
using Distributed
using LinearAlgebra
using LinearAlgebra.BLAS
import LinearAlgebra.mul!, LinearAlgebra.norm, LinearAlgebra.dot,
     LinearAlgebra.diagm, LinearAlgebra.eigmin, LinearAlgebra.eigmax
using Random
using Statistics
using MLBase

export Response, Predictors, RawData, get_X, get_Z, get_Y, contr,
    add_intercept, remove_intercept, shuffle_rows, shuffle_cols,
    cd!, cd_active!, ista!, fista!, fista_bt!, admm!,
    Mlmnet, mlmnet, coef, coef_2d, predict, fitted, resid,
    mlmnet_perms,
    make_folds, make_folds_conds, Mlmnet_cv, mlmnet_cv,
    avg_mse, lambda_min, avg_prop_zero, mlmnet_cv_summary

# Helper functions
include("std_helpers.jl")
include("mlmnet_helpers.jl")

# L1 algorithms
include("cd.jl")
include("ista.jl")
include("fista.jl")
include("fista_bt.jl")
include("admm.jl")

# Top level functions that call L1 algorithms using warm starts
include("mlmnet.jl")
# Predictions and residuals
include("predict.jl")

# Permutations
include("mlmnet_perms.jl")

# Cross-validation
include("mlmnet_cv_helpers.jl")
include("mlmnet_cv.jl")
include("mlmnet_cv_summary.jl")

end


text:
"
  Name     : c4_16_plot_julia.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 1/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com

plot(rand(5,5),color="grep",linewidth=2,title="Simple plot")
#plot(rand(5,5),color="grep",linewidth=2,title=name)
"

# Pkg.add("Plots")

using Plots
srand(123)
name="Simple plot"
plot(rand(5,5),color="gray",linewidth=2,title=name)


text:
module IncrementalInference

info("Using multithreaded convolutions Threads.nthreads()=$(Threads.nthreads())")

import Base: convert
import HDF5: root
import Distributions: sample
import Base: rand, rand!
import KernelDensityEstimate: kde!

using
  Graphs,
  NLsolve,
  Optim,
  Distributions,
  StatsBase,
  KernelDensityEstimate,
  HDF5,
  JLD,
  ProgressMeter,
  DocStringExtensions,
  Compat

export
  # pass through from Graphs.jl
  # plot,

  # added methods to functions from KernelDensityEstimate
  kde!,
  getPoints,

  # pass through functions commonly used lower down
  Npoints,
  Ndim,
  getBW,

  evalLikelihood,

  # data layer variables
  dlapi,
  localapi,
  showcurrentdlapi,
  # and callback setting function
  setdatalayerAPI!,
  DataLayerAPI,

  # general types for softtyping of variable nodes
  InferenceVariable,
  ContinuousScalar,
  ContinuousMultivariate,

  # using either dictionary or cloudgraphs
  VariableNodeData,
  PackedVariableNodeData,
  FactorMetadata,
  encodePackedType,
  FunctionNodeData,
  PackedFunctionNodeData,
  encodePackedType,
  decodePackedType,
  normalfromstring,
  categoricalfromstring,
  extractdistribution,

  FactorGraph,
  addNode!,
  addFactor!,
  doautoinit!,
  resetData!,
  getVert,
  getData,
  setData!,
  getVarNode,
  getVal,
  setVal!,
  getNumPts,
  getBWVal,
  setBW!,
  setValKDE!,
  isInitialized,
  ensureAllInitialized!,
  updateFullVert!,
  getOutNeighbors,
  BayesTree,
  EasyMessage,
  NBPMessage,
  ExploreTreeType,
  emptyFactorGraph,
  subgraphShortestPath,
  subgraphFromVerts,
  getEliminationOrder,
  buildBayesNet!,
  emptyBayesTree,
  buildTree!,
  prepBatchTree!,
  wipeBuildNewTree!,
  whichCliq,
  getKDE,
  getVertKDE,
  initializeNode!,
  batchSolve!,

  #functors need
  getSample,
  freshSamples!,

  #Visualization
  writeGraphPdf,
  ls,
  lsf,
  ls2,
  hasOrphans,
  getfnctype,
  drawCopyFG,

  # Tree stuff
  # spyCliqMat,
  evalPotential,
  evalFactor2,
  approxConv,

  # weiged sampling
  AliasingScalarSampler,
  rand!,
  rand,

  # dev
  CommonConvWrapper, # new wrapper (experimental) -- not ready for use

  # is deprecated
  FastGenericRoot,
  FastRootGenericWrapParam,
  GenericWrapParam,

  # solve inference
  inferOverTree!,
  inferOverTreeR!,

  #development interface
  upMsgPassingRecursive,

  # Inference types
  InferenceType,
  PackedInferenceType,
  Singleton,
  Pairwise,
  # introduced for approximate convolution operations
  SingleThreaded,
  MultiThreaded,

  # functor abstracts
  FunctorInferenceType,
  FunctorPairwise,
  FunctorPairwiseMinimize,
  FunctorSingleton,
  # FunctorPartialSingleton,
  FunctorPairwiseNH,
  FunctorSingletonNH,

  # Solving utils
  findRelatedFromPotential,
  shuffleXAltD!,
  numericRoot,
  numericRootGenericRandomized,
  numericRootGenericRandomizedFnc,
  numericRootGenericRandomizedFnc!,

  # user functions
  proposalbeliefs,
  predictbelief,

  # generic marginal used during elimitation game
  GenericMarginal,
  PackedGenericMarginal,

  uppA,
  convert, # for protobuf stuff
  compare,
  extractdistribution,

  # factor graph operating system utils (fgos)
  convert2packedfunctionnode,
  encodefg,
  decodefg,
  savejld,
  loadjld,
  landmarks,

  # For 1D example, should be refactored and renamed
  Odo,
  odoAdd,
  PackedOdo,
  Obsv2,
  PackedObsv2,
  Ranged,
  PackedRanged,

  # dev exports
  addGraphsVert!,
  makeAddEdge!,

  # define evalPotential functions outside IIF
  registerCallback!,

  # development
  shuffleXAltD,
  reshapeVec2Mat

  # deprecated
  # VNDencoder,
  # VNDdecoder,
  # FNDencode,
  # FNDdecode

const VoidUnion{T} = Union{Void, T}

include("FactorGraphTypes.jl")
include("AliasScalarSampling.jl")
include("DataLayerAPI.jl")
include("FactorGraph01.jl")
include("DispatchPackedConversions.jl")
include("FGOSUtils.jl")

include("JunctionTree.jl")
include("GraphConstraintTypes.jl")
include("SolverUtilities.jl")
include("TreePotentials01.jl")
include("ExplicitDiscreteMarginalizations.jl")
include("ApproxConv.jl")
include("SolveTree01.jl")

include("deprecated.jl")

# function plot(fg::FactorGraph)
#   Graphs.plot(fg.g)
# end

end


text:
using SolveBio
using Base.Test

@test SolveBio.login()["is_active"]

# Depository tests

println("Running depository tests...")
@test SolveBio.Depository.list()["total"] > 240
@test SolveBio.Depository.get("ClinVar")["id"] == 223
@test SolveBio.Depository.get(223)["name"] == "ClinVar"
@test SolveBio.Depository.versions("ClinVar")["total"] > 24
println("Successfully ran all depository tests")

# Dataset tests

println("Running dataset tests...")
@test SolveBio.Dataset.list()["total"] > 200
ds = SolveBio.Dataset.get("ClinVar/ClinVar")
@test ds["depository"] == "ClinVar"
@test ds["name"] == "ClinVar"
@test ds["documents_count"] > 1e5
ds = SolveBio.Dataset.search("ClinVar/ClinVar")
@test ds["total"] > 1e5
@test ds["genome_build"] == "GRCh37"
println("Successfully ran all dataset tests")


text:
"""
    hartley_transformation(pts::AbstractArray{T}) where T<:HomogeneousPoint

Returns a matrix which can be used to map a set of ``d``-dimensional Cartesian
points  which are represented by ``\\text{length-}(d+1)`` homogeneous coordinates into a
data-dependent coordinate system. In the data-dependent coordinate system the
origin is the center of mass (centroid) of the points  and the root-mean-square
distance of the points to the origin is equal to ``\\sqrt{d}``.

# Details

A point in ``\\mathbb{R}^d`` with
Cartesian coordinates  ``\\left(m_1, m_2, \\ldots, m_d \\right)`` can also be
expressed in homogeneous coordinates with the vector ``\\mathbf{m} =
\\left[m_1, m_2, \\ldots, m_d , 1 \\right]^\\top``.

Suppose one has a set ``\\left \\{ \\mathbf{m}_n \\right \\}_{n = 1}^{N} `` of
Cartesian points which are represented by homogeneous coordinates.
Let
```math
\\overline{\\mathbf{m}} = \\frac{1}{N} \\sum_{n = 1}^{N} \\mathbf{m}_n
\\quad  \\text{and} \\quad
\\sigma = \\left( \\frac{1}{d \\times n} \\sum_{n = 1}^{N}  \\left \\| \\mathbf{m}_n -
\\overline{\\mathbf{m}} \\right \\|^{2} \\right)^{1/2}
```
represent the centroid of the points and the root-mean-square distance of the
points to the centroid, respectively.

This function returns the matrix
```math
\\mathbf{T} =
\\begin{bmatrix}
\\sigma^{-1} & 0              &   0     & \\ldots        & -\\sigma^{-1} \\overline{m}_1 \\\\
           0 & \\sigma^{-1}   &   0     & \\ldots        & -\\sigma^{-1} \\overline{m}_2 \\\\
           0 & 0              & \\ddots &  0             &  \\vdots                      \\\\
           \\vdots & \\vdots  &   0     &  \\sigma^{-1}  & -\\sigma^{-1} \\overline{m}_d \\\\
           0 & 0              &   0     &           0    &                             1
\\end{bmatrix}
```
such that a transformed point ``\\tilde{\\mathbf{m}}_n = \\mathbf{T} \\mathbf{m}_n``
has a root-mean-square distance to the origin of a new coordinate system
equal to ``\\sqrt{d}``.

"""
function hartley_transformation(ℳ::Vector{T})::SMatrix where T <:AbstractArray
    if isempty(ℳ)
        throw(ArgumentError("Array cannot be empty."))
    end
    npts = length(ℳ)
    ndim = length(ℳ[1])
    𝐜 = centroid(ℳ)
    σ = root_mean_square(ℳ, 𝐜)
    σ⁻¹ = 1 / σ
    𝐓 = SMatrix{ndim+1,ndim+1,Float64, (ndim+1)^2}([σ⁻¹*Matrix{Float64}(I,ndim,ndim) -σ⁻¹*𝐜 ; zeros(1,ndim) 1.0])
end

function centroid(positions::Vector{T}) where T <: AbstractArray
    x = zeros(T)
    for pos ∈ positions
        x = x + pos
    end
    return x / length(positions)
end

function root_mean_square(ℳ::Vector{T}, 𝐜::T ) where  T <: AbstractArray
    total = 0.0
    npts = length(ℳ)
    ndim = length(ℳ[1])
    for 𝐦 ∈ ℳ
         total  = total + ∑((𝐦-𝐜).^2)
    end
    σ = √( (1/(npts*ndim)) * total)
end

"""
    hartley_normalization(pts::AbstractArray{T}) where T<:HomogeneousPoint

Maps a set of ``d``-dimensional Cartesian points  which are represented by
``\\text{length-}(d+1)`` homogeneous coordinates into a data-dependent coordinate system.
In the data-dependent coordinate system the origin is the center of mass
(centroid) of the points  and the root-mean-square distance of the points to the
origin is equal to ``\\sqrt{d}``.

# Details

A point in ``\\mathbb{R}^d`` with
Cartesian coordinates  ``\\left(m_1, m_2, \\ldots, m_d \\right)`` can also be
expressed in homogeneous coordinates with the vector ``\\mathbf{m} =
\\left[m_1, m_2, \\ldots, m_d , 1 \\right]^\\top``.

Suppose one has a set ``\\left \\{ \\mathbf{m}_n \\right \\}_{n = 1}^{N} `` of
Cartesian points which are represented by homogeneous coordinates.
Let
```math
\\overline{\\mathbf{m}} = \\frac{1}{N} \\sum_{n = 1}^{N} \\mathbf{m}_n
\\quad  \\text{and} \\quad
\\sigma = \\left( \\frac{1}{d \\times n} \\sum_{n = 1}^{N}  \\left \\| \\mathbf{m}_n -
\\overline{\\mathbf{m}} \\right \\|^{2} \\right)^{1/2}
```
represent the centroid of the points and the root-mean-square distance of the
points to the centroid, respectively.

This function returns a new set of points ``\\left \\{ \\tilde{\\mathbf{m}}_n \\right \\}_{n = 1}^{N} ``
where ``\\tilde{\\mathbf{m}}_n = \\mathbf{T} \\mathbf{m}_n`` for each ``n``, and
```math
\\mathbf{T} =
\\begin{bmatrix}
\\sigma^{-1} & 0              &   0     & \\ldots        & -\\sigma^{-1} \\overline{m}_1 \\\\
           0 & \\sigma^{-1}   &   0     & \\ldots        & -\\sigma^{-1} \\overline{m}_2 \\\\
           0 & 0              & \\ddots &  0             &  \\vdots                      \\\\
           \\vdots & \\vdots  &   0     &  \\sigma^{-1}  & -\\sigma^{-1} \\overline{m}_d \\\\
           0 & 0              &   0     &           0    &                             1
\\end{bmatrix}.
```
These new points have the property that their root-mean-square distance to
the origin of the coordinate system is equal to ``\\sqrt{d}``.

"""
function hartley_normalization(ℳ::Vector{<:AbstractArray})
    𝒪, 𝐓 = hartley_normalization!(copy(ℳ))
end

function hartley_normalization!(ℳ::Vector{<:AbstractArray})
    𝐓 = hartley_transformation(ℳ)
    map!(ℳ , ℳ) do 𝐦
         hom⁻¹(𝐓 * hom(𝐦))
    end
     ℳ, 𝐓
end


text:
using GaussianProcesses, BenchmarkTools
using Random, Printf

Random.seed!(1)

const d = 10       # Input observation dimension
const nt = 100     # Number of training points
const λ = rand(d)
const X = 2π * rand(d, nt) # Training data
const cK = Array{Float64}(undef, nt, nt)

const kernelswithtags = Dict((k => ["simple", "isotropic"] for k in
                                (SEIso(1.0, 1.0), Mat12Iso(1.0, 1.0), Mat32Iso(1.0, 1.0),
                                 Mat52Iso(1.0, 1.0), RQIso(1.0, 1.0, 1.0),
                                 Periodic(1.0, 1.0, 1.0)))...,
                             (k => ["simple", "ARD"] for k in
                                (SEArd(λ, 1.0), Mat12Ard(λ, 1.0), Mat32Ard(λ, 1.0),
                                 Mat52Ard(λ, 1.0), RQArd(λ, 1.0, 1.0)))...,
                             (k => ["simple", "non-stationary"] for k in
                                (LinIso(0.0), LinArd(λ)))...,
                             SEIso(1.0, 1.0) + RQIso(1.0, 1.0, 1.0) =>
                                ["composite", "isotropic + isotropic"],
                             SEIso(1.0, 1.0) + LinIso(1.0) =>
                                ["composite", "isotropic + non-stationary"],
                             Mat32Iso(1.0, 1.0) * Periodic(1.0, 1.0, 1.0) =>
                                ["composite", "isotropic * isotropic"],
                             RQArd(λ, 1.0, 1.0) * LinIso(1.0) =>
                                ["composite", "ARD * non-stationary"])
                       # SEIso(1.0, 1.0) + Noise(1.0) => ["composite", "isotropic + Noise"])

# Define benchmarks
const suite = BenchmarkGroup(["kernels"])

for (kernel, tags) in kernelswithtags
    data = GaussianProcesses.KernelData(kernel, X)
    np = GaussianProcesses.num_params(kernel)
    stack = Array{Float64}(undef, nt, nt, np)
    group = suite[string(typeof(kernel))] = BenchmarkGroup(tags)
    group["data"] = @benchmarkable GaussianProcesses.KernelData($kernel, $X)
    group["cov"] =  @benchmarkable GaussianProcesses.cov!($cK, $kernel, $X, $data)
    group["grad"] = @benchmarkable GaussianProcesses.grad_stack!($stack, $kernel, $X, $data)
end

println("Benchmark suite:")
# TODO: Workaround until https://github.com/JuliaCI/BenchmarkTools.jl/issues/96 is fixed
show(stdout, MIME"text/plain"(), suite; verbose = true, limit = max(3, length(suite)))
println()

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `suite` every time the file is included.
const paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals, :samples)
else
    println("\nTuning benchmarks...")
    tune!(suite, verbose=true)
    BenchmarkTools.save(paramspath, params(suite))
    println()
end

println("\nRunning benchmarks...")
# Use @tagged to select more specific tests:
#const test_group = suite[@tagged "non-stationary"]
#const results = run(test_group, verbose = true, seconds = 1)
const results = run(suite, verbose = true, seconds = 1)

# Save results
println("\nSaving results...")
BenchmarkTools.save(joinpath(dirname(@__FILE__), "results.json"), results)

# Print table of results
function printpadright(text, total)
    print(text)
    for _ in (length(text)+1):total
        print(" ")
    end
end

function printresults(results)
    # Compute max number of chars in the first column with kernel names
    maxchar = max(6, maximum(length(k) for k in keys(results)))

    # Print header
    printpadright("kernel", maxchar)
    println(" | data [μs] | cov [μs]  | grad [μs]")
    foreach(x -> print("-"), 1:maxchar)
    println("-+-----------+-----------+----------")

    # Print results for every kernel
    foreach(results) do (key, result)
        printpadright(key, maxchar)
        @printf(" | %9.4f | %9.4f | %9.4f\n", time(result["data"]) / 1000,
                time(result["cov"]) / 1000, time(result["grad"]) / 1000)
    end
end

println("\nResults:")
printresults(results)


text:
using DocSeeker
import Markdown

handle("searchdocs") do data
  @destruct [
    needle = query,
    mod || "Main",
    name_only = nameOnly || false,
    exportedonly = exportedOnly || false,
    allPackages || false
  ] = data
  return searchdocs′(
    needle,
    # kwargs to be passed to `DocSeeker.searchdocs`:
    # TODO: configuration for `maxreturns`
    !allPackages, mod, exportedonly, name_only
  )
end

function searchdocs′(needle, loaded = true, mod = "Main", exportedonly = false, name_only = false)
  searchmod = getmodule(mod)

  items, actualmod = _searchdocs(needle, loaded, searchmod, exportedonly, name_only)
  results = processdocs(items)

  # erase module input if the actual searched module has been changed
  push!(results, :shoulderase => searchmod ≠ actualmod)

  return results
end

function _searchdocs(needle, loaded = true, mod = Main, exportedonly = false, name_only = false)
  identifiers = split(needle, '.')
  head = string(identifiers[1])
  if loaded && head ≠ needle && (nextmod = getfield′(mod, head)) isa Module
    # if `head` is a module, update `needle` and `mod`
    nextneedle = join(identifiers[2:end], '.')
    return _searchdocs(nextneedle, loaded, nextmod, exportedonly, name_only)
  end

  return @errs searchdocs(
    needle;
    loaded = loaded, mod = mod, exportedonly = exportedonly, name_only = name_only
  ), mod
end

function processdocs(items)
  return if items isa EvalError
    errstr = sprint(showerror, items.err)
    err = startswith(errstr, "Please regenerate the") ?
            """
            Please regenerate the documentation cache with the button above.
            Note that this might take minutes. You can track the progress with the progress bar
            in the lower left corner, and should be able to use Juno normally during that time.
            """ : errstr
    Dict(
      :error => true,
      :errmsg => err
    )
  else
    Dict(
      :error => false,
      :items => [renderitem(i[2]) for i in items],
      :scores => [i[1] for i in items]
    )
  end
end

function renderitem(x)
  r = Dict(f => getfield(x, f) for f in fieldnames(DocSeeker.DocObj))
  r[:html] = view(renderMD(x.html))

  mod = getmodule(x.mod)
  name = Symbol(x.name)
  r[:typ], r[:icon], r[:nativetype] = if (name !== :ans || mod === Base) && iskeyword(name)
    "keyword", "k", x.typ
  else
    val = getfield′(mod, name)
    # @NOTE: DocSeeker can show docs for non-loaded packages via `createdocsdb()`
    nativetype = isundefined(val) ? "Undefined or not loaded yet" : x.typ
    wstype(mod, name, val), wsicon(mod, name, val), nativetype
  end
  r
end

handle("moduleinfo") do data
  @destruct [mod] = data
  moduleinfo(mod)
end

function moduleinfo(mod)
  d, items = getmoduleinfo(mod)
  items = [renderitem(i) for i in items]
  Dict(:doc => view(d), :items => items)
end

getmoduleinfo(mod) = ispackage(mod) ? packageinfo(mod) : modinfo(mod)
ispackage(mod) = Base.find_package(mod) ≠ nothing

function packageinfo(mod)
  path = DocSeeker.readmepath(mod)
  readme = ispath(path) ? String(read(path)) : ""
  description = Markdown.parse(Symbol(mod) ∈ stdlib_names ? "## Standard library package `$(mod)`" : readme)

  return  Hiccup.div(
            renderMD(description),
            renderMD("\n---\n## Defined symbols in `$(mod)`:")
          ), modulesymbols(mod)
end

function modinfo(mod)
  header = "## "
  header *= if mod in ("Core", "Base", "Main")
    "Toplevel module:"
  elseif first(split(mod, '.')) == "Core"
    "Core sub module:"
  elseif first(split(mod, '.')) == "Base"
    "Base sub module:"
  else
    "Module:"
  end
  header *= " `$(mod)`"
  header *= "\n---\n## Defined symbols:"

  return renderMD(header), modulesymbols(mod)
end

function modulesymbols(mod)
  syms = filter(x -> x.mod == mod, DocSeeker.alldocs())
  @inbounds sort!(syms, by = x -> x.name)[1:min(100,length(syms))]
end

using Logging: with_logger
using .Progress: JunoProgressLogger

function regeneratedocs()
  with_logger(JunoProgressLogger()) do
    @errs DocSeeker.createdocsdb()
  end
end
handle(regeneratedocs, "regeneratedocs")


text:
#
# Bernoulli Sampler
#
# parameters:
# 	(1) 0 <= p <= 1
#

type BernoulliSampler <: DiscreteDistributionSampler
	paramranges::Vector{Tuple{Float64,Float64}}
	distribution::Bernoulli
	function BernoulliSampler(params=Float64[])
		s = new(Tuple{Float64,Float64}[(0.0,1.0)])
		setparams!(s, isempty(params) ? Float64[0.5] : params)
		s
	end
end

function setparams!(s::BernoulliSampler, params)
	checkparamranges(s, params)
	p = params[1]
	s.distribution = Bernoulli(p)
end

getparams(s::BernoulliSampler) = [s.distribution.p]


text:
export charλ,
       charA,
       charB
"""
charλ(q,ν,k)

char value λ_(ν+k) for Mathieu's equation

y'' + (λ_(ν+k) - 2 q cos( 2z )) y = 0

where

q ∈ ℝ       - parameter
ν ∈ [-1,1]  - fractional part of the non-integer order
k ∈ ℤ⁺      - range of integer parts of the order

"""
function charλ(q::Real, nu_::Real; k::UnitRange=1:1) # reduced = true
    #nu = reduced ? rem(nu_+1,2)-1 : nu_;
    nu = rem(nu_+1,2)-1;

    # Set matrix size using formula from Shirts paper (1993), Eqs. (2.1)-(2.2).
    nu0 = nu + maximum(k)
    C = (8.46 + 0.444*nu0)/(1 + 0.085*nu0)
    D =  (0.24 + 0.0214*nu0)/(1 + 0.059*nu0)
    N = ceil(Int, (nu0 + 2 + C*abs(q)^D)/2) # matrix size is 2N+1

    (two, q2, nu2) = float.(promote(2, q, nu))
    d0 = (two .* (-N:N) .- nu2).^2
    d1 = q2 .* ones(eltype(q2), 2 * N)
    A = SymTridiagonal(d0, d1)
    a = eigvals(A, k)
    return a
end

"""
charA(q; k=0:4)

char value A_k for Mathieu's equation

y'' + (A_k - 2 q cos( 2z )) y = 0

where

q ∈ ℝ  - parameter
k ∈ ℤ⁺ - eigenvalue index

"""
function charA(q::Real; k::UnitRange=1:1)
    all(x -> x >= 0, k) || throw(DomainError(k, "Indices must be non-negative integers."))

    # Boolean indices of even and odd n values
    ie = map(iseven, k)
    io = map(!, ie)

    a = Array{Float64}(undef ,length(k))
    k1 = k .+ 1
    a[ie] = charλ(abs(q), 0.0; k = k1)[ie]
    if q>=0
        a[io] = charλ(q, one(q); k = k1)[io]
    else
        if 0 in k # maybe not the cleanest way to do it
            a[io] = charλ(abs(q), one(q); k = k[2]:last(k))[io[2:end]]
        else
            a[io] = charλ(abs(q), one(q); k=k)[io]
        end
    end
    return a
end

"""
charB(q,k)

char value B_k for Mathieu's equation

y'' + (B_k - 2 q cos( 2z )) y = 0

where

q ∈ ℝ  - parameter
k ∈ ℤ  - eigenvalueindex

"""
function charB(q::Real; k::UnitRange=1:1)
    all(x -> x > 0, k) || throw(DomainError(k, "Indices must be positive integers."))
    # Boolean indices of even and odd n values
    ie = map(iseven, k)
    io = map(!, ie)

    b = Array{Float64}(undef, length(k))
    b[ie] = charλ(q,0.0,k=k)[ie]
    if q>=0
        b[io] = charλ(q,1.0,k=k)[io]
    else
        b[io] = charλ(abs(q),1.0,k = (k .+ 1))[io]
    end
    return b
end


text:
stock_time_series_functions = [:time_series_daily, :time_series_daily_adjusted, :time_series_weekly, :time_series_weekly_adjusted, :time_series_monthly, :time_series_monthly_adjusted]

stock_time_series_functions_test = vcat(:time_series_intraday, stock_time_series_functions[1:MAX_TESTS])

@testset "Stock Time Series" begin
    for f in stock_time_series_functions_test
        @eval begin
            testname = string($f)
            @testset "$testname" begin
                @testset "JSON" begin
                    data = $f("MSFT")
                    @test typeof(data) === Dict{String,Any}
                    @test length(data) === 2
                end

                @testset "CSV" begin
                    data = $f("MSFT", datatype="csv")
                    @test typeof(data) === Tuple{Array{Any, 2}, Array{AbstractString, 2}}
                    @test length(data) === 2
                end

            end
        end
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end
end


text:
# Structured variational inference

using NumericExtensions

# update variational parameter `h` given expectation `γ`
function _update_svp!(fhmm::FHMM,
                      Y::AbstractMatrix,
                      h::Array{Float64, 3}, # shape (K, T, M)
                      γ::Array{Float64, 3}, # shape (K, T, M)
                      C⁻¹::Matrix{Float64},  # shape (D, D)
                      Δ::Matrix{Float64}, # shape (K, M)
                      Ỹ::Array{Float64, 3} # shape (K, M, T)
                      )
    const D, T = size(Y)
    # alias
    W = fhmm.W

    for m=1:fhmm.M
        Δ[:,m] = diag(W[:,:,m]'*C⁻¹*W[:,:,m])
    end

    # Eq. (12b)
    ΣWγ = zeros(Float64, D)
    for m=1:fhmm.M
        for t=1:T
            fill!(ΣWγ, 0.0)
            for l=1:fhmm.M
                if m != l
                    ΣWγ += W[:,:,l] * γ[:,t,l]
                end
            end
            Ỹ[:,m,t] = Y[:,t] - ΣWγ
        end
    end
    @assert !any(isnan(ΣWγ))

    # Eq. (12a)
    for m=1:fhmm.M
        WᵗC⁻¹ = W[:,:,m]'C⁻¹
        for t=1:T
            h[:,t,m] = WᵗC⁻¹*Ỹ[:,m,t] - 0.5Δ[:,m]
            # h[:,t,m] .-= maximum(h[:,t,m])
            # h[:,t,m] = softmax(h[:,t,m])
        end
    end
    h = exp(h)

    @assert !any(isnan(h))

    return h
end

# update expectations using forward-backward algorithm
function _updateE!(fhmm::FHMM,
                   m::Int,
                   Y::AbstractMatrix,    # shape: (D, T)
                   α::Array{Float64, 3},   # shape: (K, T, M)
                   β::Array{Float64, 3},   # shape: (K, T, M)
                   γ::Array{Float64, 3},   # shape: (K, T, M)
                   ξ::Array{Float64, 4}, # shape: (K, K, T-1, M)
                   B::Array{Float64, 3})   # shape: (K, T, M)
    const D, T = size(Y)
    const K = fhmm.K

    # scaling paramter
    c = Array(Float64, T)

    # forward recursion
    α[:,1,m] = fhmm.π[:,m] .* B[:,1,m]
    α[:,1,m] /= (c[1] = sum(α[:,1,m]) + ϵ)
    for t=2:T
        @inbounds α[:,t,m] = (fhmm.P[:,:,m]' * α[:,t-1,m]) .* B[:,t,m]
        @inbounds α[:,t,m] /= (c[t] = sum(α[:,t,m]) + ϵ)
    end
    @assert !any(isnan(α[:,:,m]))

    # backword recursion
    β[:,T,m] = 1.0
    for t=T-1:-1:1
        β[:,t,m] = fhmm.P[:,:,m] * β[:,t+1,m] .* B[:,t+1,m] ./ c[t+1]
    end
    @assert !any(isnan(β[:,:,m]))

    γ[:,:,m] = α[:,:,m] .* β[:,:,m]

    for t=1:T-1
        ξ[:,:,t,m] = fhmm.P[:,:,m] .* α[:,t,m] .* β[:,t+1,m]' .*
            B[:,t+1,m]' ./ c[t+1]
    end

    return γ, ξ
end

# update variational paramter until convergence and then update expectations.
function updateE!(fhmm::FHMM,
                  Y::AbstractMatrix,
                  α::Array{Float64, 3}, # inplace
                  β::Array{Float64, 3}, #
                  γ::Array{Float64, 3}, # shape (K, T, M)
                  ξ::Array{Float64, 4}, # inplace
                  h::Array{Float64, 3};  # shape (K, T, M)
                  maxiter::Int=10,
                  tol::Float64=1.0e-4,
                  verbose::Bool=false)
    const D, T = size(Y)
    const K = fhmm.K
    const M = fhmm.M

    C⁻¹::Matrix{Float64} = fhmm.C^-1
    Δ = Array(Float64, K, M)
    Ỹ = Array(Float64, D, M, T)

    hᵗ = copy(h)
    for iter=1:maxiter
        # update variational parameters
        hᵗ = _update_svp!(fhmm, Y, hᵗ, γ, C⁻¹, Δ, Ỹ)

        # check if converged
        diff = vecnorm(hᵗ-h)
        if diff < tol
            if verbose
                println("$(diff): converged at #$(iter).")
            end
            break
        end

        # update h with new one
        h = copy(hᵗ)
    end

    # update expectations using forward-backward algorithm for
    # each factorized HMM
    for m=1:M
        γ, ξ = _updateE!(fhmm, m, Y, α, β, γ, ξ, h)
    end

    return γ, ξ, h
end

# M-step
function updateM!(fhmm::FHMM,
                  Y::AbstractMatrix,
                  γ::Array{Float64, 3}, # shape (K, T, M)
                  ξ::Array{Float64, 4}; # shape (K, K, T-1, M)
                  fzero::Bool=false)
    const D, T = size(Y)
    const M = fhmm.M
    const K = fhmm.K

    # Eq. (A.4) update initial state prob.
    for m=1:M
        fhmm.π[:,m] = γ[:,1,m] / sum(γ[:,1,m] + ϵ)
    end

    # Eq. (A.5) update transition matrices
    for m=1:fhmm.M
        fhmm.P[:,:,m] = sum(ξ[:,:,:,m], 3) ./ (sum(γ[:,1:end-1,m], 2) + ϵ)
        # fhmm.P[:,:,m] ./= sum(fhmm.P[:,:,m], 1)
    end

    # Eq. (A.3) update observation means
    γʳ = Array(Float64, K*M, T)
    for t=1:T
       γʳ[:,t] = reshape(γ[:,t,:], K*M)
    end
    fhmm.W[:,:,:] = reshape((Y*γʳ') * pinv(γʳ*γʳ'), D, K, M)

    # sparseness contraint for minimum norm col of W
    if fzero
        for m=1:fhmm.M
            minidx = indmin([norm(fhmm.W[:,k,m]) for k=1:fhmm.K])
            fhmm.W[:,minidx,m] .*= 0.99
        end
    end

    # Eq. (A,7) update C
    s = zeros(D, D)
    for t=1:T
        for m=1:fhmm.M
            s += fhmm.W[:,:,m]*γ[:,t,m]*Y[:,t]'
        end
    end
    fhmm.C[:,:] = 1/T * (Y*Y' - s)

    nothing
end

# TODO derive correct lower bound
function bound_sv(fhmm::FHMM,
                  Y::AbstractMatrix,
                  γ::Array{Float64, 3},
                  ξ::Array{Float64, 4})
    const D, T = size(Y)
    likelihood::Float64 = 0.0
    logγ::Array{Float64, 3} = log(γ + ϵ)
    logP::Array{Float64, 3} = log(fhmm.P + ϵ)

    for t=1:T
        for m=1:fhmm.M
            likelihood += (γ[:,t,m]'logγ[:,t,m])[1]
        end
    end

    likelihood -= T*D*log(pi)
    d = det(fhmm.C)
    if d > 0
        likelihood -= 0.5T*log(d)
    end

    C⁻¹::Matrix{Float64} = fhmm.C^-1
    for t=1:T
        ws = sum([fhmm.W[:,:,m]*γ[:,t,m] for m=1:fhmm.M])
        likelihood -= (0.5*(Y[:,t] - ws)'*C⁻¹*(Y[:,t] - ws))[1]
    end

    for m=1:fhmm.M
        likelihood += (γ[:,1,m]'log(fhmm.π[:,m] + ϵ))[1]
    end

    #=
    for t=2:T
        for m=1:fhmm.M
            likelihood += (θ[:,m,t]'logP[:,:,m]*θ[:,m,t-1])[1]
        end
    end
    =#

    @assert !isnan(likelihood)

    return likelihood
end

# Result of structured variational (SV) inference
immutable SVInferenceResult
    likelihoods::Vector{Float64}
    h::Array{Float64, 3}
    α::Array{Float64, 3}
    β::Array{Float64, 3}
    γ::Array{Float64, 3}
    ξ::Array{Float64, 4}
end

function fit_sv!(fhmm::FHMM,
                 Y::AbstractMatrix; # observation matrix, shape (D, T)
                 maxiter::Int=100,
                 tol::Float64=1.0e-5,
                 fzero::Bool=false,
                 verbose::Bool=false)
    const D, T = size(Y)
    const M = fhmm.M
    const K = fhmm.K

    # means of observation probability
    α = Array(Float64, K, T, M)
    β = Array(Float64, K, T, M)
    γ = ones(Float64, K, T, M)/K
    ξ = ones(Float64, K, K, T-1, M)/K
    h = ones(Float64, K, T, M)/K

    likelihood::Vector{Float64} = zeros(1)

    # Roop of EM algorithm
    for iter=1:maxiter
        # compute bound
        score = bound_sv(fhmm, Y, γ, ξ)

        # update expectations
        γ, ξ, h = updateE!(fhmm, Y, α, β, γ, ξ, h, verbose=verbose)

        # update parameters of FHMM
        updateM!(fhmm, Y, γ, ξ, fzero=fzero)

        improvement = (score - likelihood[end]) / abs(likelihood[end])

        if verbose
            println("#$(iter): bound $(likelihood[end])
                    improvement: $(improvement)")
        end

        # check if converged
        #=
        if improvement < 1.0e-7
            if verbose
                println("#$(iter) converged")
            end
            break
        end
        =#

        push!(likelihood, score)
    end

    # remove initial zero
    shift!(likelihood)

    return SVInferenceResult(likelihood, h, α, β, γ, ξ)
end


text:
function transposition(n)

  reps = parse(Int64, ENV["LAMP_REPS"])

  m = Int(floor(n / 10))
  A = randn(m, m)
  B = randn(m, m)
  C = randn(m, m)

  Benchmarker.add_data(csv, "tr_nn_implicit", Benchmarker.measure(reps, no_transpose, A, B))
  Benchmarker.add_data(csv, "tr_nn_explicit", Benchmarker.measure(reps, g_no_transpose, A, B))
  Benchmarker.add_data(csv, "tr_tn_implicit", Benchmarker.measure(reps, first_transpose, A, B))
  Benchmarker.add_data(csv, "tr_tn_explicit", Benchmarker.measure(reps, g_first_transpose, A, B))
  Benchmarker.add_data(csv, "tr_nt_implicit", Benchmarker.measure(reps, second_transpose, A, B))
  Benchmarker.add_data(csv, "tr_nt_explicit", Benchmarker.measure(reps, g_second_transpose, A, B))
  Benchmarker.add_data(csv, "tr_tt_implicit", Benchmarker.measure(reps, both_transpose, A, B))
  Benchmarker.add_data(csv, "tr_tt_explicit", Benchmarker.measure(reps, g_both_transpose, A, B))

  @test(isapprox(no_transpose(A, B)[2], g_no_transpose(A, B)[2]))
  @test(isapprox(first_transpose(A, B)[2], g_first_transpose(A, B)[2]))
  @test(isapprox(second_transpose(A, B)[2], g_second_transpose(A, B)[2]))
  @test(isapprox(both_transpose(A, B)[2], g_both_transpose(A, B)[2]))
end

function no_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,1), size(B,2))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = A * B

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function first_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,2), size(B,2))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = transpose(A)*B

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function second_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,1), size(B,1))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = A * transpose(B)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function both_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,2), size(B,1))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = transpose(A)*transpose(B)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function g_no_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,1), size(B,2))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  LinearAlgebra.BLAS.gemm!('N', 'N', 1.0, A, B, 0.0, result)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function g_first_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,2), size(B,2))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  LinearAlgebra.BLAS.gemm!('T', 'N', 1.0, A, B, 0.0, result)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function g_second_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,1), size(B,1))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  LinearAlgebra.BLAS.gemm!('N', 'T', 1.0, A, B, 0.0, result)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function g_both_transpose(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,2), size(B,1))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  LinearAlgebra.BLAS.gemm!('T', 'T', 1.0, A, B, 0.0, result)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end


text:
# abstract ParameterUpdater
# abstract ParameterUpdaterState

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

"Allows for global storage of a ParameterUpdater, so that you don't need to pass it around."
type CurrentUpdater
  updater::ParameterUpdater
end

const _current_updater = CurrentUpdater(AdaMaxUpdater())

"Get the current global gradient updater used for gradient updates."
current_updater() = _current_updater.updater

"Set the current global gradient updater used for gradient updates."
current_updater!(updater::ParameterUpdater) = (_current_updater.updater = updater)

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

"Allows for global storage of a ModelLoss, so that you don't need to pass it around."
type CurrentModelLoss
  mloss::ModelLoss
end

const _current_mloss = CurrentModelLoss(L2DistLoss())

"Get the current global gradient mloss used for gradient updates."
current_mloss() = _current_mloss.mloss

"Set the current global gradient mloss used for gradient updates."
current_mloss!(mloss::ModelLoss) = (_current_mloss.mloss = mloss)

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

"Allows for global storage of a ParameterLoss, so that you don't need to pass it around."
type CurrentParameterLoss
  ploss::ParameterLoss
end

const _current_ploss = CurrentParameterLoss(NoParameterLoss())

"Get the current global gradient ploss used for gradient updates."
current_ploss() = _current_ploss.ploss

"Set the current global gradient ploss used for gradient updates."
current_ploss!(ploss::ParameterLoss) = (_current_ploss.ploss = ploss)

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

"Enacts a strategy to adjust the learning rate"
abstract LearningRateModel

immutable FixedLearningRate <: LearningRateModel end
OnlineStats.fit!(lrmodel::FixedLearningRate, err::Float64) = nothing

"Adapts learning rate based on relative variance of the changes in the test error"
immutable AdaptiveLearningRate <: LearningRateModel
  updater::ParameterUpdater
  errordiff::Diff
  diffvar::Variance
  adjustmentPct::Float64
  cutoffRatio::Float64
end

function AdaptiveLearningRate(updater::ParameterUpdater,
                              adjustmentPct = 1e-2,
                              cutoffRatio = 1e-1;
                              wgt = ExponentialWeight(20))
  AdaptiveLearningRate(updater, Diff(), Variance(wgt), adjustmentPct, cutoffRatio)
end

# if the error is decreasing at a large rate relative to the variance, increase the learning rate (speed it up)
function OnlineStats.fit!(lrmodel::AdaptiveLearningRate, err::Float64)
  fit!(lrmodel.errordiff, err)
  fit!(lrmodel.diffvar, diff(lrmodel.errordiff))
  m = mean(lrmodel.diffvar)
  s = std(lrmodel.diffvar)
  if s > 0.0
    pct = lrmodel.adjustmentPct * (m / s < -lrmodel.cutoffRatio ? 1.0 : -1.0)
    lrmodel.updater.η *= (1.0 + pct)
  end
end


text:
module OrderStatistics

using Iterators
using Distributions
using IndependentRandomSequences

import Base:length,rand,mean,max,min,maximum,minimum,range,sort
import Distributions:_rand!,pdf,_logpdf,logpdf,cdf,quantile,insupport
import Combinatorics:combinations

export AbstractOrderStatistic,AbstractScalarOrderStatistic,AbstractJointOrderStatistic,
      ScalarOrderStatistic,JointOrderStatistic,
      parentsequence,order,orders,
      orderstatistic,orderstatistics,jointorderstatistic,jointorderstatistics,
      spacing,spacings,jointspacing,jointspacing

include("abstract.jl")
include("helpers.jl")
include("ScalarOrderStatistic.jl")
include("JointOrderStatistic.jl")
include("generic.jl")

dnames = ["uniform","exponential"]
for dname in dnames
    include(joinpath("distributions", "$(dname).jl"))
end

end # module


text:
export DeterministicDistributionModel, get_states, get_actions

import ReinforcementLearningEnvironments: observation_space, action_space

"""
    DeterministicDistributionModel(table::Array{Vector{NamedTuple{(:nextstate, :reward, :prob),Tuple{Int,Float64,Float64}}}, 2})

Store all the transformations in the `table` field.
"""
struct DeterministicDistributionModel <: AbstractDistributionBasedModel
    table::Array{
        Vector{NamedTuple{(:nextstate, :reward, :prob),Tuple{Int,Float64,Float64}}},
        2,
    }
end

observation_space(m::DeterministicDistributionModel) = DiscreteSpace(size(m.table, 1))
action_space(m::DeterministicDistributionModel) = DiscreteSpace(size(m.table, 2))

(m::DeterministicDistributionModel)(s::Int, a::Int) = m.table[s, a]


text:
# Tools for initial digits stuff

char2digit(c::Char) = c - '0'

function first_digit(x::Real, alert::Bool = false)
    if x == 0
        if alert
            warn("Received 0, returning 0")
        end
        return 0
    end

    x = abs(x) # make sure it's positive

    while x < 1
        x *= 10
    end

    c = first(string(x))
    return char2digit(c)
end

function first_counts(x::Array{T,1}) where {T<:Real}
    digs = map(first_digit, x)
    n = length(digs)
    counts = zeros(Int, 9)
    for d in digs
        if d != 0
            counts[d] += 1
        end
    end
    return counts
end

function first_hists(x::Array{T,1}) where {T<:Real}
    counts = first_counts(x)
    hist = zeros(9)
    S = sum(counts)
    for k = 1:9
        hist[k] = counts[k] / S
    end
    return counts, hist
end

function report(x::Array{T,1}) where {T<:Real}
    counts, hist = first_hists(x)
    for k = 1:9
        println(
            k,
            "\t",
            counts[k],
            "\t",
            round(hist[k] * 100, digits = 1),
            "%\t",
            log(10, k + 1) - log(10, k),
        )
    end
end

function report_array(x::Array{T,1}) where {T<:Real}
    counts, hist = first_hists(x)
    A = Array(Any, (9, 3))
    for k = 1:9
        A[k, 1] = k
        A[k, 2] = counts[k]
        A[k, 3] = round(hist[k] * 100, digits = 2)
    end
    return A
end

function experiment()
    n = 1000
    x = rand(n)
    println("Initial distribution")
    report(x)
    println("\n\n")
    sleep(1)

    count = 0
    while count <= 1000
        m = minimum(x)
        count += 1
        if m < 1
            x *= 100.0
        end
        x = x .* rand(n)
        if count % 25 == 0
            println("step = ", count)
            report(x)
            println("-------------------------------")
            sleep(0.1)
        end
    end
end

function digit_split(n::Int)
    map(char2digit, collect(string(n)))
end

# create a 9^d multiplication table
function mult_table_old(d::Int)
    vals = Int[]
    for n = 10^(d-1):10^d-1
        digs = digit_split(n)
        v = prod(digs)
        # println(digs, " --> ", v)
        if v != 0
            push!(vals, v)
        end
    end
    return vals
end

function mult_table(d::Int)
    tic()
    counts = zeros(Int, 9)
    for n = 10^(d-1):10^d-1
        digs = digit_split(n)
        v = prod(digs)
        if v != 0
            d = first_digit(v)
            counts[d] += 1
        end
    end
    toc()
    return counts
end

mult_report_old(d::Int) = report(mult_table(d))

function mult_report(d::Int)
    counts = mult_table(d)
    N = sum(counts)

    A = Array(Any, (9, 3))
    for k = 1:9
        A[k, 1] = k
        A[k, 2] = counts[k]
        A[k, 3] = round(100 * counts[k] / N, 2)
    end
    return A
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using ParallelAccelerator
using DocOpt
if VERSION >= v"0.6.0-pre"
using SpecialFunctions
end

#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)
#CompilerTools.LambdaHandling.set_debug_level(3)
#CompilerTools.LivenessAnalysis.set_debug_level(3)

@acc begin

@inline function cndf2(in)
    out = 0.5 .+ 0.5 .* erf.(0.707106781 .* in)
    return out
end

function blackscholes(sptprice, strike, rate, volatility, time)
    logterm = log10.(sptprice ./ strike)
    powterm = .5 .* volatility .* volatility
    den = volatility .* sqrt.(time)
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den
    d2 = d1 .- den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike .* exp.(- rate .* time)
    c1 = futureValue .* NofXd2
    call = sptprice .* NofXd1 .- c1
    put  = call .- futureValue .+ sptprice
end

end

function run(iterations)
    sptprice   = Float64[ 42.0 for i = 1:iterations ]
    initStrike = Float64[ 40.0 + (i / iterations) for i = 1:iterations ]
    rate       = Float64[ 0.5 for i = 1:iterations ]
    volatility = Float64[ 0.2 for i = 1:iterations ]
    time       = Float64[ 0.5 for i = 1:iterations ]

    tic()
    put = blackscholes(sptprice, initStrike, rate, volatility, time)
    t = toq()
    println("checksum: ", sum(put))
    return t
end

function main()
    doc = """black-scholes.jl

Black-Scholes option pricing model.

Usage:
  black-scholes.jl -h | --help
  black-scholes.jl [--iterations=<iterations>]

Options:
  -h --help                  Show this screen.
  --iterations=<iterations>  Specify a number of iterations [default: 10000000].
"""
    arguments = docopt(doc)

    iterations = parse(Int, arguments["--iterations"])

    srand(0)

    println("iterations = ",iterations)

    tic()
    blackscholes(Float64[], Float64[], Float64[], Float64[], Float64[])
    println("SELFPRIMED ", toq())

    time = run(iterations)
    println("rate = ", iterations / time, " opts/sec")
    println("SELFTIMED ", time)

end

main()


text:
###############################################################################
#
# Tolerance for root solution
#
###############################################################################
"""
Abstract tolerance type

$(TYPEDEF)

"""
abstract type AbstractTolerance{FT} end

"""
Tolerance for target function residual

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
struct ResidualTolerance{FT} <: AbstractTolerance{FT}
    "Tolerance for residual"
    tol::FT
    "limit of iterations"
    n_limit::Int
end

"""
Tolerance for solution

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
struct SolutionTolerance{FT} <: AbstractTolerance{FT}
    "Tolerance for solution"
    tol::FT
    "limit of iterations"
    n_limit::Int
end

"""
Tolerance for 2D and above solution

$(TYPEDEF)

# Fields

$(TYPEDFIELDS)

"""
struct SolutionToleranceND{FT} <: AbstractTolerance{FT}
    "Tolerance for solution"
    tol::Array{FT,1}
    "limit of iterations"
    n_limit::Int
end

###############################################################################
#
# Judge if root is found
#
###############################################################################
"""
Determine whether to stopping finding the solution depending on the tolerance
    type.

$(METHODLIST)

"""
function if_break end

"""
When the tolerance is for target function residual, if the modeled residual is
    lower than the given tolerance, or if the iteration exceeds the maximum
    limit, a `true` will be returned.

    if_break(tol::ResidualTolerance{FT},
             x1::FT,
             x2::FT,
             y::FT,
             n::Int
    ) where {FT<:AbstractFloat}

Determine whether to break, given
- `tol` [`ResidualTolerance`](@ref) type tolerance struct
- `x1` Lower bound of x
- `x2` Upper bound of x
- `y` Residual of y
- `n` Current iteration
"""
if_break(tol::ResidualTolerance{FT},
         x1::FT,
         x2::FT,
         y::FT,
         n::Int
) where {FT<:AbstractFloat} =
    (abs(y) < tol.tol) || (n > tol.n_limit)

"""
When the tolerance is for solution, if the solution range is lower than the
    given tolerance, or if the iteration exceeds the maximum  limit, a `true`
    will be returned.

    if_break(tol::SolutionTolerance{FT},
             x1::FT,
             x2::FT,
             y::FT,
             n::Int
    ) where {FT<:AbstractFloat}

Determine whether to break, given
- `tol` [`SolutionTolerance`](@ref) type tolerance struct
- `x1` Lower bound of x
- `x2` Upper bound of x
- `y` Residual of y
- `n` Current iteration
"""
if_break(tol::SolutionTolerance{FT},
         x1::FT,
         x2::FT,
         y::FT,
         n::Int
) where {FT<:AbstractFloat} =
    (abs(x1-x2) < tol.tol) || (n > tol.n_limit)