text:
"""
Provides a basic interface for working with MDPs/POMDPs
"""
POMDPs

#################################################################
####################### Problem Model ###########################
#################################################################

"""
Base type for a POMDP model, defined by the user
"""
POMDP

"""
Base type for state, action and observation spaces
"""
AbstractSpace

"""
Base type for states
"""
State

"""
Base type for actions
"""
Action

"""
Base type for observation
"""
Observation

"""
    states(pomdp::POMDP)

Returns the complete state space of a POMDP.
"""
states

"""
    actions(pomdp::POMDP)

Returns the entire action space of a POMDP.
"""
actions(pomdp::POMDP)

"""
    actions(pomdp::POMDP, state::State, aspace::AbstractSpace)

Modifies aspace to the action space accessible from the given state and returns it.
"""
actions(pomdp::POMDP, state::State, aspace::AbstractSpace)

"""
    actions(pomdp::POMDP, belief::Belief, aspace::AbstractSpace)

Modifies aspace to the action space accessible from the states with nonzero belief and returns it.
"""
actions(pomdp::POMDP, belief::Belief, aspace::AbstractSpace)

"""
    observations(pomdp::POMDP)

Returns the entire observation space.
"""
observations(pomdp::POMDP)

"""
    observations(pomdp::POMDP, state::State, ospace::AbstractSpace)

Modifies ospace to the observation space accessible from the given state and returns it.
"""
observations(pomdp::POMDP, state::State, ospace::AbstractSpace)

"""
    reward(pomdp::POMDP, state::State, action::Action, statep::State)

Returns the immediate reward for the s-a-s' triple
"""
reward

"""
    transition(pomdp::POMDP, state::State, action::Action)

Returns the transition distribution from the current state-action pair
"""
transition(pomdp::POMDP, state::State, action::Action)

"""
    transition(pomdp::POMDP, state::State, action::Action, distribution::AbstractDistribution)

Modifies distribution to the transition distribution from the current state-action pair and returns it
"""
transition(pomdp::POMDP, state::State, action::Action, distribution::AbstractDistribution)

"""
    observation(pomdp::POMDP, state::State, action::Action, statep::State)

Returns the observation distribution for the s-a-s' tuple (state, action, and next state)
"""
observation(pomdp::POMDP, state::State, action::Action, statep::State)

"""
    observation(pomdp::POMDP, state::State, action::Action, statep::State, distribution::AbstractDistribution)

Modifies distribution to the observation distribution for the s-a-s' tuple (state, action, and next state) and returns it
"""
observation(pomdp::POMDP, state::State, action::Action, statep::State, distribution::AbstractDistribution)

"""
    discount(pomdp::POMDP)

Returns the discount factor
"""
discount

"""
    isterminal(pomdp::POMDP, s::State)

Checks if state s is terminal
"""
isterminal

#################################################################
####################### Distributions ###########################
#################################################################

"""
Base type for a probability distribution
"""
AbstractDistribution

"""
    rand(rng::AbstractRNG, d::AbstractDistribution, sample)

Fill sample with a random element from distribution d. The sample can be a state, action or observation.
"""
rand(rng::AbstractRNG, d::AbstractDistribution, sample)

"""
    rand(rng::AbstractRNG, s::AbstractSpace, sample)

Fill sample with a random element from space s. The sample can be a state, action or observation.
"""
rand(rng::AbstractRNG, s::AbstractSpace, sample)

"""
    pdf(d::AbstractDistribution, x)

Value of probability distribution function at x
"""
pdf

#################################################################
##################### Solvers and Policies ######################
#################################################################

"""
Base type for an MDP/POMDP solver
"""
Solver

"""
Base type for a policy (a map from every possible belief, or more abstract policy state, to an optimal or suboptimal action)
"""
Policy

"""
    solve(solver::Solver, pomdp::POMDP)

Solves the POMDP using method associated with solver, and returns a policy.
"""
solve(solver::Solver, pomdp::POMDP)

"""
    solve(solver::Solver, pomdp::POMDP, policy::Policy)

Solves the POMDP and modifies policy to be the solution of pomdp and returns it
"""
solve(solver::Solver, pomdp::POMDP, policy::Policy)

"""
    action(policy::Policy, belief::Belief)

Returns an action for the current belief, given the policy

"""
action(policy::Policy, belief::Belief)

"""
    action(policy::Policy, belief::Belief, a::Action)

Fills and returns action a for the current belief, given the policy
"""
action(policy::Policy, belief::Belief, a::Action)

"""
    action(policy::Policy, s::State)

Returns an action for the current state, given the policy

"""
action(policy::Policy, s::State)

"""
    action(policy::Policy, s::State, a::Action)

Fills and returns action a for the current state, given the policy
"""
action(policy::Policy, s::State, a::Action)

"""
    value(policy::Policy, s::State)

Returns the utility value from a given state
"""
value(policy::Policy, s::State)

"""
    value(policy::Policy, b::Belief)

Returns the utility value from a given belief
"""
value(policy::Policy, b::Belief)

#################################################################
############################ Belief #############################
#################################################################

"""
Base type for an object representing some knowledge about the state (often a probability distribution)
"""
Belief

"""
Base type for an object that defines how a belief should be updated
"""
BeliefUpdater

"""
    update(updater::BeliefUpdater, belief_old::Belief, action::Action, obs::Observation)

Returns a new instance of an updated belief given the old belief (belief_old) and the latest action and observation
"""
update(updater::BeliefUpdater, belief_old::Belief, action::Action, obs::Observation)

"""
    update(updater::BeliefUpdater, belief_old::Belief, action::Action, obs::Observation, belief_new::Belief)

Modifies belief_new to the belief given the old belief (belief_old) and the latest action and observation and returns
the updated belief
"""
update(updater::BeliefUpdater, belief_old::Belief, action::Action, obs::Observation, belief_new::Belief)

"""
    initial_belief(pomdp::POMDP)

Returns an initial belief for the pomdp
"""
initial_belief

"""
    convert_belief(updater::BeliefUpdater, b::Belief)

Returns a belief that can be updated using updater that has a similar distribution to b
"""
convert_belief

"""
    updater(p::Policy)
Returns a default BeliefUpdater appropriate for the passed in policy
"""
updater

#################################################################
############################ Simulation #########################
#################################################################

"""
Base type for an object defining how a simulation should be carried out
"""
Simulator

"""
    simulate(simulator::Simulator, pomdp::POMDP, policy::Policy, updater::BeliefUpdater, initial_belief::Belief)

Runs a simulation using the specified policy and returns the accumulated reward
"""
simulate

#################################################################
######################### Convenience ###########################
#################################################################

"""
    index(pomdp::POMDP, state::State)

Returns the index of the given state for a discrete POMDP
"""
index(pomdp::POMDP, s::State)

"""
    iterator(space::AbstractSpace)

Returns an iterator over a space
"""
iterator(space::AbstractSpace)

#################################################################
############################ Creators ###########################
#################################################################

"""
    create_state(pomdp::POMDP)

Creates a single state object (for preallocation purposes)
"""
create_state(pomdp::POMDP)

"""
    create_observation(pomdp::POMDP)

Creates a single observation object (for preallocation purposes)
"""
create_observation(pomdp::POMDP)

"""
    create_transition_distribution(pomdp::POMDP)

Returns a transition distribution
"""
create_transition_distribution(pomdp::POMDP)

"""
    create_observation_distribution(pomdp::POMDP)

Returns an observation distribution
"""
create_observation_distribution(pomdp::POMDP)

"""
    create_policy(solver::Solver, pomdp::POMDP)

Creates a policy object (for preallocation purposes)
"""
create_policy(solver::Solver, pomdp::POMDP)

"""
    create_action(pomdp::POMDP)

Creates an action object (for preallocation purposes)
"""
create_action(pomdp::POMDP)

"""
    create_belief(pomdp::POMDP)

Creates a belief either to be used by updater or pomdp
"""
create_belief(pomdp::POMDP)

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
using MINPACK
using Test, Printf, LinearAlgebra

#=
NOTE: the contents of this file were taken from NLsolve.jl.

The license for that package appears below.

The NLsolve.jl package is licensed under the MIT License:
Copyright (c) 2017 NLsolve.jl contributors

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

=#

# If the results should be printed to a
const PRINT_FILE = false

struct DifferentiableMultivariateFunction{Tf,Tg}
    f!::Tf
    g!::Tg
end

@testset "minpack" begin

function rosenbrock()
    function f!(fvec::Vector, x::Vector)
        fvec[1] = 1 - x[1]
        fvec[2] = 10(x[2]-x[1]^2)
    end
    function g!(fjac::Matrix, x::Vector)
        fjac[1,1] = -1
        fjac[1,2] = 0
        fjac[2,1] = -20x[1]
        fjac[2,2] = 10
    end
    (DifferentiableMultivariateFunction(f!, g!), [-1.2; 1.], "Rosenbrock")
end

function powell_singular()
    function f!(fvec::Vector, x::Vector)
        fvec[1] = x[1] + 10x[2]
        fvec[2] = sqrt(5)*(x[3] - x[4])
        fvec[3] = (x[2] - 2x[3])^2
        fvec[4] = sqrt(10)*(x[1] - x[4])^2
    end
    function g!(fjac::Matrix, x::Vector)
        fill!(fjac, 0)
        fjac[1,1] = 1
        fjac[1,2] = 10
        fjac[2,3] = sqrt(5)
        fjac[2,4] = -fjac[2,3]
        fjac[3,2] = 2(x[2] - 2x[3])
        fjac[3,3] = -2fjac[3,2]
        fjac[4,1] = 2sqrt(10)*(x[1] - x[4])
        fjac[4,4] = -fjac[4,1]
    end
    (DifferentiableMultivariateFunction(f!, g!), [3.; -1.; 0.; 1.], "Powell singular")
end

function powell_badly_scaled()
    c1 = 1e4
    c2 = 1.0001
    function f!(fvec::Vector, x::Vector)
        fvec[1] = c1*x[1]*x[2] - 1
        fvec[2] = exp(-x[1]) + exp(-x[2]) - c2
    end
    function g!(fjac::Matrix, x::Vector)
        fjac[1,1] = c1*x[2]
        fjac[1,2] = c1*x[1]
        fjac[2,1] = -exp(-x[1])
        fjac[2,2] = -exp(-x[2])
    end
    (DifferentiableMultivariateFunction(f!, g!), [0.; 1.], "Powell badly scaled")
end

function wood()
    c3 = 2e2
    c4 = 2.02e1
    c5 = 1.98e1
    c6 = 1.8e2

    function f!(fvec::Vector, x::Vector)
        temp1 = x[2] - x[1]^2
        temp2 = x[4] - x[3]^2
        fvec[1] = -c3*x[1]*temp1 - (1 - x[1])
        fvec[2] = c3*temp1 + c4*(x[2] - 1) + c5*(x[4] - 1)
        fvec[3] = -c6*x[3]*temp2 - (1 - x[3])
        fvec[4] = c6*temp2 + c4*(x[4] - 1) + c5*(x[2] - 1)
    end

    function g!(fjac::Matrix, x::Vector)
        fill!(fjac, 0)
        temp1 = x[2] - 3x[1]^2
        temp2 = x[4] - 3x[3]^2
        fjac[1,1] = -c3*temp1 + 1
        fjac[1,2] = -c3*x[1]
        fjac[2,1] = -2*c3*x[1]
        fjac[2,2] = c3 + c4
        fjac[2,4] = c5
        fjac[3,3] = -c6*temp2 + 1
        fjac[3,4] = -c6*x[3]
        fjac[4,2] = c5
        fjac[4,3] = -2*c6*x[3]
        fjac[4,4] = c6 + c4
    end
    (DifferentiableMultivariateFunction(f!, g!), [-3.; -1; -3; -1], "Wood")
end

function helical_valley()
    tpi = 8*atan(1)
    c7 = 2.5e-1
    c8 = 5e-1

    function f!(fvec::Vector, x::Vector)
        if x[1] > 0
            temp1 = atan(x[2]/x[1])/tpi
        elseif x[1] < 0
            temp1 = atan(x[2]/x[1])/tpi + c8
        else
            temp1 = c7*sign(x[2])
        end
        temp2 = sqrt(x[1]^2+x[2]^2)
        fvec[1] = 10(x[3] - 10*temp1)
        fvec[2] = 10(temp2 - 1)
        fvec[3] = x[3]
    end

    function g!(fjac::Matrix, x::Vector)
        temp = x[1]^2 + x[2]^2
        temp1 = tpi*temp
        temp2 = sqrt(temp)
        fjac[1,1] = 100*x[2]/temp1
        fjac[1,2] = -100*x[1]/temp1
        fjac[1,3] = 10
        fjac[2,1] = 10*x[1]/temp2
        fjac[2,2] = 10*x[2]/temp2
        fjac[2,3] = 0
        fjac[3,1] = 0
        fjac[3,2] = 0
        fjac[3,3] = 1
    end
    (DifferentiableMultivariateFunction(f!, g!), [-1.; 0; 0], "Helical Valley")
end

function watson(n::Integer)
    c9 = 2.9e1

    function f!(fvec::Vector, x::Vector)
        fill!(fvec, 0)
        for i = 1:29
            ti = i/c9
            sum1 = 0.0
            temp = 1.0
            for j = 2:n
                sum1 += (j-1)*temp*x[j]
                temp *= ti
            end
            sum2 = 0.0
            temp = 1.0
            for j = 1:n
                sum2 += temp*x[j]
                temp *= ti
            end
            temp1 = sum1 - sum2^2 - 1
            temp2 = 2*ti*sum2
            temp = 1/ti
            for k = 1:n
                fvec[k] += temp*(k-1-temp2)*temp1
                temp *= ti
            end
        end
        temp = x[2] - x[1]^2 - 1
        fvec[1] += x[1]*(1-2temp)
        fvec[2] += temp
    end

    function g!(fjac::Matrix, x::Vector)
        fill!(fjac, 0)
        for i = 1:29
            ti = i/c9
            sum1 = 0.0
            temp = 1.0
            for j = 2:n
                sum1 += (j-1)*temp*x[j]
                temp *= ti
            end
            sum2 = 0.0
            temp = 1.0

            for j = 1:n
                sum2 += temp*x[j]
                temp *= ti
            end
            temp1 = 2(sum1 - sum2^2 - 1)
            temp2 = 2*sum2
            temp = ti^2
            tk = 1.0

            for k = 1:n
                tj = tk
                for j = k:n
                    fjac[k,j] += tj*((Float64(k-1)/ti - temp2)*(Float64(j-1)/ti - temp2) - temp1)
                    tj *= ti
                end
                tk *= temp
            end
        end
        fjac[1,1] += 6x[1]^2 - 2x[2] + 3
        fjac[1,2] -= 2x[1]
        fjac[2,2] += 1
        for k = 1:n
            for j = k:n
                fjac[j,k] = fjac[k,j]
            end
        end
    end
    (DifferentiableMultivariateFunction(f!, g!), zeros(n), "Watson$n")
end

function chebyquad(n::Integer)
    tk = 1/n

    function f!(fvec::Vector, x::Vector)
        fill!(fvec, 0)
        for j = 1:n
            temp1 = 1.0
            temp2 = 2x[j]-1
            temp = 2temp2
            for i = 1:n
                fvec[i] += temp2
                ti = temp*temp2 - temp1
                temp1 = temp2
                temp2 = ti
            end
        end
        iev = -1.0
        for k = 1:n
            fvec[k] *= tk
            if iev > 0
                fvec[k] += 1/(k^2-1)
            end
            iev = -iev
        end
    end

    function g!(fjac::Matrix, x::Vector)
        for j = 1:n
            temp1 = 1.
            temp2 = 2x[j] - 1
            temp = 2*temp2
            temp3 = 0.0
            temp4 = 2.0
            for k = 1:n
                fjac[k,j] = tk*temp4
                ti = 4*temp2 + temp*temp4 - temp3
                temp3 = temp4
                temp4 = ti
                ti = temp*temp2 - temp1
                temp1 = temp2
                temp2 = ti
            end
        end
    end
    (DifferentiableMultivariateFunction(f!, g!), collect(1:n)/(n+1), "Chebyquad")
end

function brown_almost_linear(n::Integer)
    function f!(fvec::Vector, x::Vector)
        sum1 = sum(x) - (n+1)
        for k = 1:(n-1)
            fvec[k] = x[k] + sum1
        end
        fvec[n] = prod(x) - 1
    end

    function g!(fjac::Matrix, x::Vector)
        fill!(fjac, 1)
        fjac[diagind(fjac)] .= 2
        prd = prod(x)
        for j = 1:n
            if x[j] == 0.0
                fjac[n,j] = 1.
                for k = 1:n
                    if k != j
                        fjac[n,j] *= x[k]
                    end
                end
            else
                fjac[n,j] = prd/x[j]
            end
        end
    end
    (DifferentiableMultivariateFunction(f!, g!), 0.5*ones(n), "Brown almost-linear$n")
end

function discrete_boundary_value(n::Integer)
    h = 1/(n+1)

    function f!(fvec::Vector, x::Vector)
        for k = 1:n
            temp = (x[k] + k*h + 1)^3
            if k != 1
                temp1 = x[k-1]
            else
                temp1 = 0.0
            end
            if k != n
                temp2 = x[k+1]
            else
                temp2 = 0.0
            end
            fvec[k] = 2x[k] - temp1 - temp2 + temp*h^2/2
        end
    end

    function g!(fjac::Matrix, x::Vector)
        for k = 1:n
            temp = 3*(x[k]+k*h+1)^2
            for j = 1:n
                fjac[k,j] = 0
            end
            fjac[k,k] = 2 + temp*h^2/2
            if k != 1
                fjac[k,k-1] = -1
            end
            if k != n
                fjac[k,k+1] = -1
            end
        end
    end
    initial_x = collect(1:n)*h
    initial_x = initial_x .* (initial_x .- 1)

    (DifferentiableMultivariateFunction(f!, g!), initial_x, "Discrete boundary value")
end

function discrete_integral_equation(n::Integer)
    h = 1/(n+1)

    function f!(fvec::Vector, x::Vector)
        for k = 1:n
            tk = k*h
            sum1 = 0.0
            for j = 1:k
                tj = j*h
                sum1 += tj*(x[j] + tj + 1)^3
            end
            sum2 = 0.0
            kp1 = k+1
            if n >= kp1
                for j = kp1:n
                    tj = j*h
                    sum2 += (1-tj)*(x[j] + tj + 1)^3
                end
            end
            fvec[k] = x[k] + h*((1-tk)*sum1 + tk*sum2)/2
        end
    end

    function g!(fjac::Matrix, x::Vector)
        for k = 1:n
            tk = k*h
            for j = 1:n
                tj = j*h
                fjac[k,j] = h*min(tj*(1-tk), tk*(1-tj))*3(x[j] + tj + 1)^2/2
            end
            fjac[k,k] += 1
        end
    end

    initial_x = collect(1:n)*h
    initial_x = initial_x .* (initial_x .- 1)

    (DifferentiableMultivariateFunction(f!, g!), initial_x, "Discrete integral equation")
end

function trigonometric(n::Integer)
    function f!(fvec::Vector, x::Vector)
        for j = 1:n
            fvec[j] = cos(x[j])
        end
        sum1 = sum(fvec)
        for k = 1:n
            fvec[k] = n+k-sin(x[k]) - sum1 - k*fvec[k]
        end
    end

    function g!(fjac::Matrix, x::Vector)
        for j = 1:n
            temp = sin(x[j])
            for k = 1:n
                fjac[k,j] = temp
            end
            fjac[j,j] = (j+1)*temp - cos(x[j])
        end
    end
    (DifferentiableMultivariateFunction(f!, g!), ones(n)/n, "Trigonometric")
end

function variably_dimensioned(n::Integer)
    function f!(fvec::Vector, x::Vector)
        sum1 = 0.0
        for j = 1:n
            sum1 += j*(x[j]-1)
        end
        temp = sum1*(1+2sum1^2)
        for k = 1:n
            fvec[k] = x[k] - 1 + k*temp
        end
    end

    function g!(fjac::Matrix, x::Vector)
        sum1 = 0.0
        for j = 1:n
            sum1 += j*(x[j]-1)
        end
        temp = 1 + 6sum1^2
        for k = 1:n
            for j = k:n
                fjac[k,j] = k*j*temp
                fjac[j,k] = fjac[k,j]
            end
            fjac[k,k] += 1
        end
    end
    (DifferentiableMultivariateFunction(f!, g!), 1 .- collect(1:n)/n, "Variably dimensioned")
end

function broyden_tridiagonal(n::Integer)
    function f!(fvec::Vector, x::Vector)
        for k = 1:n
            temp = (3-2x[k])*x[k]
            if k != 1
                temp1 = x[k-1]
            else
                temp1 = 0.0
            end
            if k != n
                temp2 = x[k+1]
            else
                temp2 = 0.0
            end
            fvec[k] = temp - temp1 - 2temp2 + 1
        end
    end

    function g!(fjac::Matrix, x::Vector)
        fill!(fjac, 0)
        for k = 1:n
            fjac[k,k] = 3-4x[k]
            if k != 1
                fjac[k,k-1] = -1
            end
            if k != n
                fjac[k,k+1] = -2
            end
        end
    end
    (DifferentiableMultivariateFunction(f!, g!), -ones(n), "Broyden tridiagonal")
end

function broyden_banded(n::Integer)
    ml = 5
    mu = 1

    function f!(fvec::Vector, x::Vector)
        for k = 1:n
            k1 = max(1, k-ml)
            k2 = min(k+mu, n)
            temp = 0.0
            for j = k1:k2
                if j != k
                    temp += x[j]*(1+x[j])
                end
            end
            fvec[k] = x[k]*(2+5x[k]^2) + 1 - temp
        end
    end

    function g!(fjac::Matrix, x::Vector)
        fill!(fjac, 0)
        for k = 1:n
            k1 = max(1, k-ml)
            k2 = min(k+mu, n)
            for j = k1:k2
                if j != k
                    fjac[k,j] = -(1+2x[j])
                end
            end
            fjac[k,k] = 2+15x[k]^2
        end
    end
    (DifferentiableMultivariateFunction(f!, g!), -ones(n), "Broyden banded")
end

alltests = [rosenbrock();
            powell_singular();
            powell_badly_scaled();
            wood();
            helical_valley();
            watson(6);
            watson(9);
            chebyquad(5);
            chebyquad(6);
            chebyquad(7);
            #chebyquad(8);
            chebyquad(9);
            brown_almost_linear(10);
            brown_almost_linear(30);
            brown_almost_linear(40);
            discrete_boundary_value(10);
            discrete_integral_equation(1);
            discrete_integral_equation(10);
            trigonometric(10);
            variably_dimensioned(10);
            broyden_tridiagonal(10);
            broyden_banded(10);
           ]

if PRINT_FILE; f_out = open("minpack_results.dat", "w"); end

@printf("%-45s   %5s   %5s   %5s   %14s     %10s\n", "Function", "Dim", "NFEV",
        "NJEV", "Final inf-norm", "total time")
println("-"^86)

if PRINT_FILE
    @printf(f_out, "%-45s   %5s   %5s   %5s   %14s\n", "Function", "Dim", "NFEV",
            "NJEV", "Final inf-norm")
    println(f_out, "-"^86)
end

for (df, initial, name) in alltests
    for method in (:hybr, :lm)
        if name == "Trigonometric" && method == :hybr
            continue
        end
        if (Sys.WORD_SIZE == 32 &&
            name in ["Brown almost-linear30", "Brown almost-linear40"] &&
            method == :hybr )
            continue
        end
        tot_time = @elapsed r = fsolve(df.f!, df.g!, initial, method=method, tracing=true)
        @printf("%-45s   %5d   %5d   %5d   %14e   %10e\n", name*"-"*string(method), length(initial),
                r.trace.f_calls, r.trace.g_calls, norm(r.f, Inf), tot_time)
        # @show name, method, r.converged
        @test r.converged
        # without gradient
        if name == "Watson9" && method == :lm
            continue
        end
        tot_time = @elapsed r_no_g = fsolve(df.f!, initial, method=method, tracing=true)
        @printf("%-45s   %5d   %5d   %5d   %14e   %10e\n", name*"-"*string(method)*"-no_g",
                length(initial), r_no_g.trace.f_calls, r_no_g.trace.g_calls, norm(r_no_g.f, Inf), tot_time)
        if PRINT_FILE
            @printf(f_out, "%-45s   %5d   %5d   %5d   %14e\n", name*"-"*string(method)*"-no_g", length(initial),
                    r_no_g.trace.f_calls, r_no_g.trace.g_calls, norm(r_no_g.f, Inf))
            @printf(f_out, "%-45s   %5d   %5d   %5d   %14e\n", name*"-"*string(method), length(initial),
                    r.trace.f_calls, r.trace.g_calls, norm(r.f, Inf))
        end
        # @show name, method, r_no_g.converged
        @test r_no_g.converged
    end
end

if PRINT_FILE; close(f_out); end

end

ext: jl


text:
function powerlevel(x, y)
    serial = 7511
    rackid = x + 10
    power = rackid * y
    power += serial
    power *= rackid
    power = convert(Int, floor(power / 100)) % 10
    return power - 5
end

function foo()
    max_power = -100000000
    max_x = 0
    max_y = 0
    max_size = 0

    data = zeros(Int, 300, 300)

    for x in 1:300
        for y in 1:300
            data[x, y] = powerlevel(x, y)
        end
    end

    for x in 1:300
        for y in 1:300
            for size in 0:(300-max(x,y))
                power = 0
                for i in 0:size
                    for j in 0:size
                        power += data[x + i, y + j]
                    end
                end
                if power > max_power
                    max_power = power
                    max_x = x
                    max_y = y
                    max_size = size + 1
                end
            end
        end
    end

    return max_x, max_y, max_size
end

z = foo()
println(z[1])
println(z[2])
println(z[3])


text:
#!/usr/bin/env julia

using Test

include("parse_id3.jl")
using .ID3: Tags, read_tags_string, parse_tags, print_tags, load_genres, genre

@testset "read_tags_string" begin
  io = IOBuffer()
  ts = string("TAG", repeat("X", 128-3))
  write(io, ts)
  seekstart(io)
  @test read_tags_string(io) == ts

  @testset "invalid" begin
    io = IOBuffer()
    ts = string("TAG", "nope")
    write(io, ts)
    seekstart(io)
    @test read_tags_string(io) == ""
  end
end

@testset "parse_tags" begin
  @testset "invalid strings" begin
    @test parse_tags("") == nothing
    @test parse_tags("i'm not a tags string") == nothing
    @test parse_tags("TAG i'm still not a tags string") == nothing
  end
  @testset "song" begin
    only_song = string("TAGSome song", repeat("\0", 116))
    t = parse_tags(only_song)
    @test t != nothing
    @test t.song == "Some song"

    song1 = repeat("1", 30)
    @test parse_tags(string("TAG", song1, repeat("2", 128-33))).song == song1
  end
  @testset "full fields v1" begin
    song = repeat("S", 30)
    artist = repeat("A", 30)
    album = repeat("a", 30)
    year = "2019"
    comment = repeat("c", 30)
    genre_index = UInt8(42)
    s = string("TAG", song, artist, album, year, comment, Char(genre_index))
    @test length(s) == 128

    t = parse_tags(s)
    @test t != nothing
    @test t.song == song
    @test t.artist == artist
    @test t.album == album
    @test t.year == year
    @test t.comment == comment
    @test t.genre_index == genre_index
    @test t.track_number == nothing
  end

  @testset "partial fields v1" begin
    song = "hey"
    artist = "i"
    album = "just"
    year = "1999"
    comment = "met you"
    genre_index = UInt8(10)
    s = string("TAG",
               rpad(song, 30, '\0'),
               rpad(artist, 30, '\0'),
               rpad(album, 30, '\0'),
               year,
               rpad(comment, 30, '\0'),
               Char(genre_index))
    @test length(s) == 128

    t = parse_tags(s)
    @test t != nothing
    @test t.song == song
    @test t.artist == artist
    @test t.album == album
    @test t.year == year
    @test t.comment == comment
    @test t.genre_index == genre_index
    @test t.track_number == nothing
  end

  @testset "with track number (v1.1)" begin
    @testset "full comment" begin
      comment = repeat("C", 28)
      s = string("TAG", repeat("X", 30+30+30), "1984", comment,
                 '\0', Char(42), Char(18))
      @test length(s) == 128

      t = parse_tags(s)
      @test t != nothing
      @test t.comment == comment
      @test t.track_number == 42
      @test t.genre_index == 18
    end

    @testset "partial comment" begin
      comment = "yo i'm a comment"
      s = string("TAG", repeat("X", 30+30+30), "1984", comment,
                 repeat('\0', 28-length(comment)),
                 '\0', Char(42), Char(18))
      @test length(s) == 128

      t = parse_tags(s)
      @test t != nothing
      @test t.comment == comment
      @test t.track_number == 42
      @test t.genre_index == 18
    end
  end
end

@testset "load_genres" begin
  genres = load_genres()
  @test genres != nothing
  @test length(genres) > 100
end

@testset "genre" begin
  function tag_with_genre(n)
    Tags("", "", "", "", "", n, nothing)
  end
  @test genre(tag_with_genre(1)) == "Classic Rock"
  @test genre(tag_with_genre(0)) == "Blues"
  @test genre(tag_with_genre(200)) == nothing
end

max_stars_repo_name: bfontaine/Katas

max_issues_repo_name: bfontaine/Katas

max_forks_repo_name: bfontaine/Katas


text:
#module Plotfuncs
    using .Plots
    using LinearAlgebra
    using StatsBase
    #using ..TightBinding
    export plot_lattice_2d,calc_band_plot,plot_DOS,calc_band_plot_finite
    export plot_fermisurface_2D
    export get_DOS

    function plot_lattice_2d(lattice)
        dim = lattice.dim
        lw = 0.5
        ls = :dash
        colors = ["red","blue","orange","brown","yellow"]
        println("Plot the 2D lattice structure")
        if dim != 2
            println("Error! dim should be 2 to plot the lattice. dim is $dim")
            return
        end
        pls = plot(size=(500,500), legend=false)
        r0 = [0,0]
        r1 = lattice.vectors[1]
        r2 = lattice.vectors[2]

        na = 2
        nb = 2

        for ia = 0:na
            for ib=0:nb
                rorigin = lattice.vectors[1]*ia +  lattice.vectors[2]*ib
                plot!([r0[1]+rorigin[1], r1[1]+rorigin[1]],
                [r0[2]+rorigin[2], r1[2]+rorigin[2]], color="red", lw=lw, ls=ls)
                plot!([r0[1]+rorigin[1], r2[1]+rorigin[1]],
                [r0[2]+rorigin[2], r2[2]+rorigin[2]], color="red", lw=lw, ls=ls)
                plot!([r1[1]+rorigin[1], r1[1]+r2[1]+rorigin[1]],
                [r1[2]+rorigin[2], r1[2]+r2[2]+rorigin[2]], color="red", lw=lw, ls=ls)
                plot!([r2[1]+rorigin[1], r1[1]+r2[1]+rorigin[1]],
                [r2[2]+rorigin[2], r1[2]+r2[2]+rorigin[2]], color="red", lw=lw, ls=ls)

                for i=1:lattice.numatoms
                    x,y = get_position(lattice,lattice.positions[i])
                    plot!([x+rorigin[1]],[y+rorigin[2]],marker=:circle,color=colors[i])
                end

                for i=1:lattice.numhopps
                    ijpositions = lattice.hoppings[i].ijpositions
                    x,y = get_position(lattice,ijpositions)
                    ijatoms = lattice.hoppings[i].ijatoms
                    x0,y0 = get_position(lattice,lattice.positions[ijatoms[1]])
                    plot!([x0+rorigin[1], x0+x+rorigin[1]],
                    [y0+rorigin[2], y0+y+rorigin[2]], color="black", lw=1)
                end
            end
        end
        plot!(aspect_ratio=:equal)
        return pls
    end

    """
        plot_fermisurface_2D(lattice::Lattice;Eshift = 0.0,nk = 20)

    Show the Fermi surface for 2D system. \\
    This plots the contour with E=0.0+Eshift \\
    Eshift: the energy shift from the chemical potential lattice.μ \\
    nk: the number of meshes in 2D momentum space. The total mesh is nk x nk.
    """
    function plot_fermisurface_2D(lattice;Eshift = 0.0,nk = 50)
        ham = hamiltonian_k(lattice) #Construct the Hamiltonian
        dim = lattice.dim
        n = lattice.numatoms

        if dim != 2
            println("This function only supports 2D case. lattice.dim should be 2. But dim = $dim .")
            return
        end

        k1s = range(-0.51, stop = 0.51, length = nk)
        k2s = range(-0.51, stop = 0.51, length = nk)
        kx_vec = zeros(Float64,nk)
        ky_vec = zeros(Float64,nk)
        energies = zeros(Float64,n,nk,nk)
        colors = ["red","blue","orange","brown","yellow"]
        kxmin,kymin = get_position_kspace(lattice,[-0.5,-0.5])
        kxmax,kymax = get_position_kspace(lattice,[0.5,0.5])

        for i1=1:nk
            for i2=1:nk
                kx,ky = get_position_kspace(lattice,[k1s[i1],k2s[i2]])
                kx_vec[i1] = kx
                ky_vec[i2] = ky
                energy = dispersion(dim,n,ham,[kx,ky])
                for j=1:n
                    energies[j,i1,i2] = energy[j]
                end

            end
        end
#        println("energies ",energies[1,:,:])
        #println(kx_vec)
        i = 1
        pls = contour(kx_vec,ky_vec,energies[i,:,:],levels=[Eshift],
            xlabel="kx",ylabel="ky",aspect_ratio=1,
            xlims=(kxmin,kxmax),ylims=(kymin,kymax)
            )
        for i=2:n
            pls = contour!(kx_vec,ky_vec,energies[i,:,:],levels=[Eshift],
                xlabel="kx",ylabel="ky",aspect_ratio=1,
                xlims=(kxmin,kxmax),ylims=(kymin,kymax)
                )
        end
        return pls

    end

    function calc_band_plot(klines,lattice)
        ham = hamiltonian_k(lattice)
        numlines = klines.numlines
        dim = lattice.dim
        klength = 0.0
        vec_k = []
        numatoms = lattice.numatoms
        energies = zeros(Float64,numatoms,0)
        #plot(x,xticks=([4,5],["G","F"]))
        xticks_values = []
        xticks_labels = []

        for i=1:numlines
            push!(xticks_values,klength)
            push!(xticks_labels,klines.kpoints[i].name_start)
            kmin = klines.kpoints[i].kmin
            kmax = klines.kpoints[i].kmax
            nk = klines.kpoints[i].nk

            kmin_real = kmin[:]
            kmax_real = kmax[:]

            kdistance = sqrt(sum((kmin_real .- kmax_real).^2))
            vec_k_temp,energies_i = calc_band(kmin_real,kmax_real,nk,lattice,ham)
            energies_i = sort(energies_i,dims=1)
            vec_k_i = range(klength,length=nk,stop=klength+kdistance)
            vec_k = vcat(vec_k,vec_k_i)
            energies = hcat(energies,energies_i)
            klength += kdistance
            push!(xticks_values,klength)
            push!(xticks_labels,klines.kpoints[i].name_end)
        end

        pls = plot(vec_k[:],[energies[i,:] for i=1:numatoms],
                        xticks = (xticks_values,xticks_labels))

        return pls
    end

    function calc_band_plot_finite(klines,lattice,direction;periodic = true,nsites = 20)
        ham = hamiltonian_k_1d(lattice,direction,periodic=periodic,nsites=nsites)
        numlines = klines.numlines
        dim = lattice.dim-1
        klength = 0.0
        vec_k = Float64[]
        numatoms = lattice.numatoms
        N = numatoms*nsites
        energies = zeros(Float64,0,N)
        #plot(x,xticks=([4,5],["G","F"]))
        xticks_values = []
        xticks_labels = []

        for i=1:numlines

            push!(xticks_values,klength)
            push!(xticks_labels,klines.kpoints[i].name_start)
            kmin = klines.kpoints[i].kmin
            kmax = klines.kpoints[i].kmax
            nk = klines.kpoints[i].nk

            kmin_real = kmin[:]
            kmax_real = kmax[:]

            kdistance = sqrt(sum((kmin_real .- kmax_real).^2))

            vec_k_temp = zeros(Float64,dim,nk)

            for idim=1:dim
                k = range(kmin_real[idim],length=nk,stop=kmax_real[idim])
                vec_k_temp[idim,:] = k[:]
            end
            energies_i = zeros(Float64,nk,N)
            for ik=1:nk
                energies_i[ik,:] = eigen(ham(vec_k_temp[:,ik])).values[:]
            end

            vec_k_i = range(klength,length=nk,stop=klength+kdistance)
            vec_k = vcat(vec_k,vec_k_i)
            energies = vcat(energies,energies_i)
            klength += kdistance
            push!(xticks_values,klength)
            push!(xticks_labels,klines.kpoints[i].name_end)
        end

        pls = plot(vec_k[:],energies,
                        xticks = (xticks_values,xticks_labels),legend=false)

        return pls
    end

    function plot_DOS(lattice,nk;nbins=100)
        dim = lattice.dim
        hist = get_DOS(lattice,nk;nbins=nbins)
        pls = bar(hist.edges[1][1:end],hist.weights/nk^dim)
        return pls

        #=
        dim = lattice.dim
        ham = hamiltonian_k(lattice)
        n = lattice.numatoms

        if dim == 1
            energies = zeros(Float64,n,nk)
            dk = 1/(nk-1)
            for ik = 1:nk
                k = (ik-1)*dk
            end
        elseif dim == 2
            energies = zeros(Float64,n,nk,nk)
            dk1 = 1/(nk-1)
            dk2 = 1/(nk-1)
            for ik1 = 1:nk
                k1 = (ik1-1)*dk1
                for ik2 = 1:nk
                    k2 = (ik2-1)*dk2
                    kx,ky = get_position_kspace(lattice,[k1,k2])
                    energy = dispersion(dim,n,ham,[kx,ky])
                    energies[:,ik1,ik2] = energy[:]

                end
            end
        elseif dim == 3
            energies = zeros(Float64,nk,nk,nk)
            dk1 = 1/(nk-1)
            dk2 = 1/(nk-1)
            dk3 = 1/(nk-1)

            for ik1 = 1:nk
                k1 = (ik1-1)*dk1
                for ik2 = 1:nk
                    k2 = (ik2-1)*dk2
                    for ik3=1:nk
                        k3= (ik3-1)*dk3
                        kx,ky,kz = get_position_kspace(lattice,[k1,k2,k3])
                        energy = dispersion(dim,n,ham,[kx,ky,kz])
                        energies[:,ik1,ik2,ik3] = energy[:]
                    end

                end
            end
        end
        energies = vec(energies)
#            println(energies)
        pls = histogram(energies,bins=nbins,normalize=true)

        return pls
        =#
    end

#end


text:
using DDD
using Test
using DDD: makeInstanceDict, inclusiveComparison

cd(@__DIR__)

@testset "Geometry" begin
    arr = Int[3; 4; 6]
    @test isapprox(internalAngle(arr[1]), π / 3)
    @test isapprox(internalAngle(arr[2]), π / 2)
    @test isapprox(internalAngle(arr[3]), 2π / 3)
    @test isapprox(externalAngle(arr[1]), 2π / 3)
    @test isapprox(externalAngle(arr[2]), π / 2)
    @test isapprox(externalAngle(arr[3]), π / 3)
    xyz = [1.0, 0.0, 0.0]
    θ = pi / 2
    uvw = [0.0, 5.0, 0.0]
    abc = [0.0, 0.0, 0.0]
    p = rot3D(xyz, uvw, abc, θ)
    @test isapprox(p, [0.0, 0.0, -1.0])
    uvw = [0.0, 0.0, 20.0]
    abc = [0.0, 0.0, 0.0]
    p = rot3D(xyz, uvw, abc, θ)
    @test isapprox(p, [0.0, 1.0, 0.0])
    uvw = [1.0, 0.0, 0.0]
    abc = [0.0, 0.0, 0.0]
    p = rot3D(xyz, uvw, abc, θ)
    @test isapprox(p, xyz)
    xyz = [-23.0, 29.0, -31.0]
    uvw = [11.0, -13.0, 17.0]
    abc = [-2.0, 5.0, 7.0]
    θ = 37 / 180 * pi
    p = rot3D(xyz, uvw, abc, θ)
    @test isapprox(p, [-21.1690, 31.0685, -30.6029]; atol = 1e-4)
    @test compStruct(1, 1.2) == false

    planenorm = Float64[0, 0, 1]
    planepnt = Float64[0, 0, 5]
    raydir = Float64[0, -1, -2]
    raypnt = Float64[0, 0, 10]

    ψ = linePlaneIntersect(planenorm, planepnt, raydir, raypnt)
    @test isapprox(ψ, [0, -2.5, 5.0])

    planenorm = Float64[0, 2, 1]
    planepnt = Float64[0, 0, 5]
    raydir = Float64[0, -1, -2]
    raypnt = Float64[0, 0, 10]
    ψ = linePlaneIntersect(planenorm, planepnt, raydir, raypnt)
    @test isapprox(ψ, [0.0, -1.25, 7.5])

    planenorm = Float64[0, 0, 1]
    planepnt = Float64[0, 0, 5]
    raydir = Float64[0, 1, 2]
    raypnt = Float64[0, 0, 10]

    ψ = linePlaneIntersect(planenorm, planepnt, raydir, raypnt)
    @test isapprox(ψ, [0, -2.5, 5.0])

    planenorm = Float64[0, 0, 1]
    planepnt = Float64[0, 0, 5]
    raydir = Float64[0, 1, -2]
    raypnt = Float64[0, 0, 10]

    ψ = linePlaneIntersect(planenorm, planepnt, raydir, raypnt)
    @test isapprox(ψ, [0, 2.5, 5.0])

    planenorm = Float64[0, 0, 1]
    planepnt = Float64[0, 0, 5]
    raydir = Float64[0, 1, 0]
    raypnt = Float64[0, 0, 5]
    ψ = linePlaneIntersect(planenorm, planepnt, raydir, raypnt)
    @test isinf(ψ)

    planenorm = Float64[0, 0, 1]
    planepnt = Float64[0, 0, 5]
    raydir = Float64[0, 1, 0]
    raypnt = Float64[0, 0, 6]
    ψ = linePlaneIntersect(planenorm, planepnt, raydir, raypnt)
    @test isnothing(ψ)

    x0, x1 = zeros(3), ones(3)
    y0, y1 = zeros(3), zeros(3)
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 0, 0)

    x0, x1 = zeros(3), ones(3)
    y0, y1 = ones(3), ones(3)
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 1, 0)

    y0, y1 = zeros(3), ones(3)
    x0, x1 = zeros(3), zeros(3)
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 0, 0)

    y0, y1 = zeros(3), ones(3)
    x0, x1 = ones(3), ones(3)
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 0, 1)

    x0, x1 = zeros(3), ones(3)
    y0, y1 = 0.5 * ones(3), 0.5 * ones(3)
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 0.5, 0)

    y0, y1 = zeros(3), ones(3)
    x0, x1 = 0.5 * ones(3), 0.5 * ones(3)
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 0, 0.5)

    x0, x1 = zeros(3), ones(3)
    y0, y1 = zeros(3), ones(3) .+ eps(Float64)^2
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 0, 0)

    x0, x1 = zeros(3), ones(3)
    y0, y1 = [1, 0, 0], [2, 1, 1]
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (1, 0, 0, 0)

    x1, x0 = zeros(3), ones(3)
    y0, y1 = [1, 0, 0], [2, 1, 1]
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (6, 0, 1, 1)

    x0, x1 = zeros(3), ones(3)
    y1, y0 = [1, 0, 0], [2, 1, 1]
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (2, 0, 1, 1)

    x1, x0 = zeros(3), ones(3)
    y0, y1 = [1, 1, 0], [1.5, 0.5, 0.5]
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0.5, 0, 0.25, 0.5)

    x0, x1 = zeros(3), zeros(3)
    y0, y1 = zeros(3), zeros(3)
    vx0, vx1 = zeros(3), zeros(3)
    vy0, vy1 = zeros(3), zeros(3)
    distSq, dDistSqDt, L1, L2 = minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    @test (distSq, dDistSqDt, L1, L2) == (0, 0, 0, 0)
end

@testset "Auxiliary" begin
    dict = Dict(
        "intFixDln" => nodeTypeDln(2),
        "noneDln" => nodeTypeDln(0),
        "intMobDln" => nodeTypeDln(1),
        "srfFixDln" => nodeTypeDln(4),
        "extDln" => nodeTypeDln(5),
        "srfMobDln" => nodeTypeDln(3),
        "tmpDln" => nodeTypeDln(6),
    )

    @test makeInstanceDict(nodeTypeDln) == dict
    data = rand(5)
    @test inclusiveComparison(data[rand(1:5)], data...)
    @test !inclusiveComparison(data, data[rand(1:5)] * 6)
end

@testset "Quadrature" begin
    n = 17
    a = -13
    b = 23
    x, w = gausslegendre(n, a, b)

    bma = (b - a) * 0.5
    bpa = (b + a) * 0.5
    x1, w1 = gausslegendre(n)
    x1 = bma * x1 .+ bpa
    w1 = bma * w1

    @test x ≈ x1
    @test w ≈ w1
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
immutable TestBackend
    flag::Int64;
    function TestBackend(n::Int64)
        @assert n>=1 && n<=5;
        return new(n);
    end
end

const testopencl  =TestBackend(1);
const testopenmp  =TestBackend(2);
const testpthreads=TestBackend(3);
const testcuda    =TestBackend(4);
const testserial  =TestBackend(5);

function (==)(rhs::TestBackend,lhs::TestBackend)
    return (rhs.flag==lhs.flag);
end

function get_info_string(backend::TestBackend)

    pinnedCores = [0,0,1,1];

    Serial_Info   ="mode = Serial";
    OpenMP_Info   = "mode = OpenMP  , schedule = compact, chunk = 10";
    OpenCL_Info   = "mode = OpenCL  , platformID = 0, deviceID = 0";
    CUDA_Info     = "mode = CUDA    , deviceID = 0";
    Pthreads_Info = "mode = Pthreads, threadCount = 4, schedule = compact, pinnedCores = [0, 0, 1, 1]";
    #COI_Info      = "mode = COI     , deviceID = 0";

    linfo::String = "";
    if backend == testopencl
        if OCCA.USE_OPENCL
            linfo=OpenCL_Info;
        else
            print("OpenCL support not compiled in OCCA.\n");
            linfo=Serial_Info;
        end
    end
    if backend == testopenmp
        if OCCA.USE_OPENMP
            linfo=OpenMP_Info;
        else
            print("OpenMP support not compiled in OCCA.\n");
            linfo=Serial_Info;
        end
    end
    if backend == testcuda
        if OCCA.USE_CUDA
            linfo = CUDA_Info;
        else
            print("CUDA support not compiled in OCCA.\n");
            linfo = Serial_Info;
        end
    end
    if backend == testpthreads
        if OCCA.USE_PTHREADS
            linfo = Pthreads_Info;
        else
            print("Pthreads support not compiled in OCCA.\n");
            linfo = Serial_Info;
        end
    end

    #OpenMP will default to serial implementation if it is not otherwise built into OCCA.
    if backend == testserial
        linfo = Serial_Info;
    end

    return linfo;
end


text:
using Base.Test, JSON
using ComparativeAutograder

return3good_code = "return3() = 3"

submission_name, submission_fp  = mktemp()
solution_name, solution_fp  = mktemp()
write.([submission_fp, solution_fp], return3good_code)
close.([submission_fp, solution_fp])
testsuite = TestSuite(
    "return3",
    # two tests, each with no args and no kwargs
    [TestCase(()), TestCase(())]
)

test_stdout, test_stderr = @capture_stdstreams begin
    runtestsuite(testsuite, submission_name, submission_name)
end
println(test_stderr)
rm.([submission_name, solution_name])
result = JSON.parse(test_stdout)


text:
module DynamicPolynomials

using Reexport
@reexport using MultivariatePolynomials
const MP = MultivariatePolynomials

include("var.jl")
include("mono.jl")
const DMonomialLike{C} = Union{Monomial{C}, PolyVar{C}}
include("term.jl")
include("monovec.jl")
include("poly.jl")
const TermPoly{C, T} = Union{Term{C, T}, Polynomial{C, T}}
const PolyType{C} = Union{Polynomial{C}, Term{C}, Monomial{C}, PolyVar{C}}
MP.variable_union_type(::Union{PolyType{C}, Type{<:PolyType{C}}}) where {C} = PolyVar{C}
MP.constantmonomial(::Type{<:PolyType{C}}) where {C} = Monomial{C}()
MP.constantmonomial(p::PolyType) = Monomial(_vars(p), zeros(Int, nvariables(p)))
MP.monomialtype(::Type{<:PolyType{C}}) where C = Monomial{C}
MP.monomialtype(::PolyType{C}) where C = Monomial{C}
#function MP.constantmonomial(::Type{Monomial{C}}, vars=PolyVar{C}[]) where {C}
#    return Monomial{C}(vars, zeros(Int, length(vars)))
#end
MP.termtype(::Union{TermPoly{C, T}, Type{<:TermPoly{C, T}}}) where {C, T} = Term{C, T}
MP.termtype(::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{T}) where {C, T} = Term{C, T}
MP.polynomial(p::PolyType) = Polynomial(p)
function MP.polynomial(p::PolyType{C}, ::Type{T}) where {C, T}
    return convert(Polynomial{C, T}, p)
end
MP.polynomialtype(::Type{Term{C, T}}) where {T, C} = Polynomial{C, T}
MP.polynomialtype(::Type{T}, ::Type{<:DMonomialLike{C}}) where {T, C} = Polynomial{C, T}
MP.polynomialtype(::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{T}) where {C, T} = Polynomial{C, T}
_vars(p::AbstractArray{<:PolyType}) = mergevars(_vars.(p))[1]
MP.variables(p::Union{PolyType, MonomialVector, AbstractArray{<:PolyType}}) = _vars(p) # tuple(_vars(p))
MP.nvariables(p::Union{PolyType, MonomialVector, AbstractArray{<:PolyType}}) = length(_vars(p))
MP.similarvariable(p::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{Val{V}}) where {C, V} = PolyVar{C}(string(V))
MP.similarvariable(p::Union{PolyType{C}, Type{<:PolyType{C}}}, V::Symbol) where {C} = PolyVar{C}(string(V))
include("promote.jl")

include("operators.jl")
include("comp.jl")

include("diff.jl")
include("subs.jl")

include("div.jl")

end # module


text:
@testset "Metadata" begin
    #HalfMBB - 2x2

    # 7   8   9
    #  x--x--x      x--x--x
    #  |  |  |      |C3|C4|
    #  x--x--x      x--x--x
    # 4| 5|  |6     |C1|C2|
    #  x--x--x      x--x--x
    # 1   2   3

    #Node coords:
    #JuAFEM.Node{2,Float64}([0.0, 0.0])
    #JuAFEM.Node{2,Float64}([1.0, 0.0])
    #JuAFEM.Node{2,Float64}([2.0, 0.0])
    #JuAFEM.Node{2,Float64}([0.0, 1.0])
    #JuAFEM.Node{2,Float64}([1.0, 1.0])
    #JuAFEM.Node{2,Float64}([2.0, 1.0])
    #JuAFEM.Node{2,Float64}([0.0, 2.0])
    #JuAFEM.Node{2,Float64}([1.0, 2.0])
    #JuAFEM.Node{2,Float64}([2.0, 2.0])

    #Cells
    #JuAFEM.Cell{2,4,4}((1, 2, 5, 4))
    #JuAFEM.Cell{2,4,4}((2, 3, 6, 5))
    #JuAFEM.Cell{2,4,4}((4, 5, 8, 7))
    #JuAFEM.Cell{2,4,4}((5, 6, 9, 8))

    # 7   8   9
    #  x--x--x      x--x--x
    #  |  |  |      |C3|C4|
    #  x--x--x      x--x--x
    # 4| 5|  |6     |C1|C2|
    #  x--x--x      x--x--x
    # 1   2   3

    #Node dofs
    ##1 #2 #3  #4 #5 #6  #7  #8  #9
    #1  3   9  7  5  11  15  13  17
    #2  4  10  8  6  12  16  14  18

    #Cell dofs
    # #1  #2  #3  #4
    # 1   3   7   5
    # 2   4   8   6
    # 3   9   5  11
    # 4  10   6  12
    # 5  11  13  17
    # 6  12  14  18
    # 7   5  15  13
    # 8   6  16  14

    problem = HalfMBB(Val{:Linear}, (2,2), (1.,1.), 1., 0.3, 1.);

    coords = [(0.0, 0.0),
              (1.0, 0.0),
              (2.0, 0.0),
              (0.0, 1.0),
              (1.0, 1.0),
              (2.0, 1.0),
              (0.0, 2.0),
              (1.0, 2.0),
              (2.0, 2.0)]

    for (i, node) in enumerate(problem.ch.dh.grid.nodes)
        @test node.x.data == coords[i]
    end

    cells = [JuAFEM.Cell{2,4,4}((1, 2, 5, 4)),
        JuAFEM.Cell{2,4,4}((2, 3, 6, 5)),
        JuAFEM.Cell{2,4,4}((4, 5, 8, 7)),
        JuAFEM.Cell{2,4,4}((5, 6, 9, 8))]
    @test problem.ch.dh.grid.cells == cells

    node_dofs = [1 3  9 7 5 11 15 13 17;
                 2 4 10 8 6 12 16 14 18]
    @test problem.metadata.node_dofs == node_dofs

    cell_dofs = [1  3  7 5;
                 2  4  8 6;
                 3  9  5 11;
                 4 10  6 12;
                 5 11 13 17;
                 6 12 14 18;
                 7  5 15 13;
                 8  6 16 14]
    @test problem.metadata.cell_dofs == cell_dofs

    dof_cells = problem.metadata.dof_cells
    cell_dofs = problem.metadata.cell_dofs
    for i in 1:JuAFEM.ndofs(problem.ch.dh)
        d_cells = dof_cells[i]
        for c in d_cells
            (cellid, localdof) = c
            @test i == cell_dofs[localdof, cellid]
        end
    end

    node_first_cells = [(1, 1),
                        (1, 2),
                        (2, 2),
                        (1, 4),
                        (1, 3),
                        (2, 3),
                        (3, 4),
                        (3, 3),
                        (4, 3)]
    @test problem.metadata.node_first_cells == node_first_cells

end


text:
# Routines to generate a hierarchical rectangular grid in gmsh

function _validate_gmsh_generate_rectangular_grid_input(bb::NTuple{4, Float64},
                                                        x::Vector{Vector{Float64}},
                                                        y::Vector{Vector{Float64}})
    # bb has valid values
    xmin, xmax, ymin, ymax = bb
    dx = xmax - xmin
    dy = ymax - ymin
    if dx <= 0 || dy <= 0
        @error "Invalid bounding box. Must be of the form [xmin, xmax, ymin, ymax]"
    end
    # x and y have same number of levels
    if length(x) != length(y)
        @error "x and y must have the same length"
    end
    # monotonic x, y
    for xvector in x
        for i = 1:length(xvector)-1
            if xvector[i] > xvector[i+1]
                @error "x is not monotonically increasing."
            end
        end
    end
    for yvector in y
        for i = 1:length(yvector)-1
            if yvector[i] > yvector[i+1]
                @error "y is not monotonically increasing."
            end
        end
    end
    # x, y in bb
    for xvector in x
        for i = 1:length(xvector)
            if !(xmin <= xvector[i] <= xmax)
                errval = xvector[i]
                @error "x value ($errval) is not in the bounding box."
            end
        end
    end
    for yvector in y
        for i = 1:length(yvector)
            if !(ymin <= yvector[i] <= ymax)
                errval = yvector[i]
                @error "y value ($errval) is not in the bounding box."
            end
        end
    end

    # Add bb limits to x, y and propagate all low level (coarse) divisions to higher levels (fine)
    x_full = deepcopy(x)
    y_full = deepcopy(y)
    append!(x_full[1], [xmin, xmax])
    append!(y_full[1], [ymin, ymax])
    for i = 2:length(x_full)
        append!(x_full[i], x_full[i-1])
        append!(y_full[i], y_full[i-1])
    end
    for i = 1:length(x_full)
        sort!(unique!(x_full[i]))
        sort!(unique!(y_full[i]))
    end

    # Check that all modular geometry sizes are the same (nlevels - 1)
    nlevels = length(x_full)
    if 1 < nlevels
        dx = x_full[nlevels-1][2] - x_full[nlevels-1][1]
        dy = y_full[nlevels-1][2] - y_full[nlevels-1][1]
        for i = 2:length(x_full[nlevels-1])-1
            if !(x_full[nlevels-1][i+1] - x_full[nlevels-1][i] ≈ dx)
                @error string("Grid level $nlevels must have equal x-divisions so that all ",
                             "modular geometry have the same size!")
            end
        end
        for i = 2:length(y_full[nlevels-1])-1
            if !(y_full[nlevels-1][i+1] - y_full[nlevels-1][i] ≈ dy)
                @error string("Grid level $nlevels must have equal y-divisions so that all ",
                             "modular geometry have the same size!")
            end
        end
    end
    return x_full, y_full
end

function gmsh_generate_rectangular_grid(bb::NTuple{4, Float64},
                                        x::Vector{Vector{Float64}},
                                        y::Vector{Vector{Float64}};
                                        material::String = "MATERIAL_WATER")
    @debug "Generating rectangular grid in gmsh"
    # Validate input
    x_full, y_full = _validate_gmsh_generate_rectangular_grid_input(bb, x, y)

    # Create the grid
    # We only need to make the smallest rectangles and group them into larger ones
    grid_tags_coords = Tuple{Int32,Float64,Float64}[]
    nlevels = length(x_full)
    x_small = x_full[nlevels]
    y_small = y_full[nlevels]
    for (yi, yv) in enumerate(y_small[1:length(y_small)-1])
        for (xi, xv) in enumerate(x_small[1:length(x_small)-1])
            tag = gmsh.model.occ.add_rectangle(xv, yv, 0, x_small[xi+1] - xv, y_small[yi+1] - yv)
            push!(grid_tags_coords, (tag, xv, yv))
        end
    end
    @debug "Synchronizing model"
    gmsh.model.occ.synchronize()

    # Label the rectangles with the appropriate grid level and location
    # Create a dictionary holding all the physical group names and tags corresponding to
    # each group name.
    grid_levels_tags = Dict{String,Array{Int32,1}}()
    max_grid_digits = max(length(string(length(x_small)-1)),
                          length(string(length(y_small)-1)))
    # Create each grid name
    for lvl in 1:nlevels
        for j in 1:length(y_full[lvl])-1
            for i in 1:length(x_full[lvl])-1
                grid_str = string("GRID_L", lvl, "_", lpad(i, max_grid_digits, "0"), "_",
                                                      lpad(j, max_grid_digits, "0"))
                grid_levels_tags[grid_str] = Int32[]
            end
        end
    end
    # For each rectangle, find which grid level/index it belongs to.
    for (tag, x0, y0) in grid_tags_coords
        for lvl in 1:nlevels
            i = searchsortedlast(x_full[lvl], x0)
            j = searchsortedlast(y_full[lvl], y0)
            grid_str = string("GRID_L", lvl, "_", lpad(i, max_grid_digits, "0"), "_",
                                                  lpad(j, max_grid_digits, "0"))
            push!(grid_levels_tags[grid_str], tag)
        end
    end
    @debug "Setting rectangular grid physical groups"
    for name in keys(grid_levels_tags)
        output_tag = gmsh.model.add_physical_group(2, grid_levels_tags[name])
        gmsh.model.set_physical_name(2, output_tag, name)
    end
    tags = [ tag for (tag, x0, y0) in grid_tags_coords ]
    # If there is already a physical group with this material name, then we need to erase it
    # and make a new physical group with the same name, that contains the previous entities
    # as well as the grid entities
    groups = gmsh.model.get_physical_groups()
    for grp in groups
        name = gmsh.model.get_physical_name(grp[1], grp[2])
        if material == name
            old_tags = gmsh.model.get_entities_for_physical_group(grp[1], grp[2])
            gmsh.model.remove_physical_groups([grp])
            append!(tags, old_tags)
            break
        end
    end
    output_tag = gmsh.model.add_physical_group(2, tags)
    gmsh.model.set_physical_name(2, output_tag, material)

    # Return tags
    return tags
end

function gmsh_generate_rectangular_grid(bb::NTuple{4, Float64},
                                        nx::Vector{Int64},
                                        ny::Vector{Int64};
                                        material::String = "MATERIAL_WATER")
    if any(nx .<= 0) || any(ny .<= 0)
       @error "Can only subdivide into positive, non-zero intervals!"
    end
    x_mult = 1
    xvec = Vector{Float64}[]
    for xv in nx
        push!(xvec, collect(LinRange{Float64}(bb[1], bb[2], xv*x_mult + 1)))
        x_mult *= xv
    end
    y_mult = 1
    yvec = Vector{Float64}[]
    for yv in ny
        push!(yvec, collect(LinRange{Float64}(bb[3], bb[4], yv*y_mult + 1)))
        y_mult *= yv
    end
    tags = gmsh_generate_rectangular_grid(bb, xvec, yvec; material = material)
    return tags
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    supportfun(geometry, direction)

Support function of `geometry` for given `direction`.

## References

* Gilbert, E., Johnson, D., Keerthi, S. 1988. [A fast
  Procedure for Computing the Distance Between Complex
  Objects in Three-Dimensional Space]
  (https://ieeexplore.ieee.org/document/2083)
"""
function supportfun(g::Geometry, d::Vec)
  v = vertices(g)
  c = centroid(g)
  i = argmax([(vᵢ - c) ⋅ d for vᵢ in v])
  v[i]
end

function supportfun(b::Ball, d::Vec)
  r = (radius(b) / norm(d)) * d
  center(b) + r
end


text:
using VLConstraintBasedModelGenerationUtilities
using DataFrames

# setup path to sequence file -
path_to_protein_sequence_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/P-MZ373340.fasta"

# load -
protein_table = build_protein_table(path_to_protein_sequence_file) |> check
seq = protein_table[1,:protein_sequence]
seq_name = protein_table[1,:id]

# build the translation table -
translation_dictionary = Dict{String,DataFrame}()
(number_of_proteins, _) = size(protein_table)
for protein_index = 1:number_of_proteins

    local_seq = protein_table[protein_index,:protein_sequence]
    local_seq_name = protein_table[protein_index,:id]
    translation_table = build_translation_reaction_table(local_seq_name, local_seq) |> check

    # grab -
    translation_dictionary[local_seq_name] = translation_table
end


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

"""
    zigzagenc(x::Signed)::Unsigned

Zigzag-encode x, Google Protocol Buffers compatible.
"""
function zigzagenc end
export zigzagenc

@inline zigzagenc(x::Signed) = unsigned(xor((x << 1), (x >> (8 * sizeof(x) - 1))))

"""
    zigzagdec(x::Unsigned)::Signed

Zigzag-decode x, Google Protocol Buffers compatible.
"""
function zigzagdec end
export zigzagdec

@inline zigzagdec(x::Unsigned) = signed(xor((x >>> 1), (-(x & 1))))


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function write(filename::AbstractString, solution::AbstractDict)::Nothing
    open(filename, "w") do file
        return JSON.print(file, solution, 2)
    end
    return
end


text:
using Documenter, DiffOpt

makedocs(;
    modules=[DiffOpt],
    doctest = false,
    clean = true,
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.MathJax()
    ),
    pages=[
        "Home" => "index.md",
        "Introduction" => "intro.md",
        "Manual" => "manual.md",
        "Usage" => "usage.md",
        "Reference" => "reference.md",
        "Examples" => [
            "Solving an LP" => "solve-LP.md",
            "Solving a QP" => "solve-QP.md",
            "Solving conic with PSD and SOC constraints" => "solve-conic-1.md",
            "Differentiating a simple QP by hand" => "matrix-inversion-manual.md"
        ]
    ],
    strict = true,  # See https://github.com/JuliaOpt/JuMP.jl/issues/1576
    repo="https://github.com/jump-dev/DiffOpt.jl/blob/{commit}{path}#L{line}",
    sitename="DiffOpt.jl",
    authors="JuMP Community"
)

deploydocs(
    repo   = "github.com/jump-dev/DiffOpt.jl.git",
)


text:
using DiffRules, SpecialFunctions, NaNMath

for (M, f, arity) in DiffRules.diffrules()
  arity == 1 || continue
  @eval begin
    @adjoint $M.$f(x::Real) = $M.$f(x),
      Δ -> (Δ * $(DiffRules.diffrule(M, f, :x)),)
  end
end

for (M, f, arity) in DiffRules.diffrules()
  arity == 2 || continue
  da, db = DiffRules.diffrule(M, f, :a, :b)
  @eval begin
    @adjoint $M.$f(a::Real, b::Real) = $M.$f(a, b),
      Δ -> (Δ * $da, Δ * $db)
  end
end

@adjoint Base.convert(T::Type{<:Real}, x::Real) = convert(T, x), Δ -> (nothing, Δ)

for T in Base.uniontypes(Core.BuiltinInts)
    @adjoint (::Type{T})(x::Core.BuiltinInts) = T(x), Δ -> (Δ,)
end

@adjoint Base.:+(xs...) = +(xs...), Δ -> map(_ -> Δ, xs)

@adjoint function sincos(x)
  s, c = sincos(x)
  (s, c), ((s̄, c̄),) -> (s̄*c - c̄*s,)
end

@adjoint a // b = (a // b, c̄ -> (c̄ * 1//b, - c̄ * a // b // b))

@nograd floor, ceil, trunc, round, hash


text:
# ---
# title: 1508. Range Sum of Sorted Subarray Sums
# id: problem1508
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/range-sum-of-sorted-subarray-sums/description/>
# hidden: true
# ---
#
# Given the array `nums` consisting of `n` positive integers. You computed the
# sum of all non-empty continous subarrays from the array and then sort them in
# non-decreasing order, creating a new array of `n * (n + 1) / 2` numbers.
#
# _Return the sum of the numbers from index_`left` _to index_`right` ( **indexed
# from 1** ) _, inclusive, in the  new array. _Since the answer can be a huge
# number return it modulo 10^9 + 7.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4], n = 4, left = 1, right = 5
#     Output: 13
#     Explanation: All subarray sums are 1, 3, 6, 10, 2, 5, 9, 3, 7, 4. After sorting them in non-decreasing order we have the new array [1, 2, 3, 3, 4, 5, 6, 7, 9, 10]. The sum of the numbers from index le = 1 to ri = 5 is 1 + 2 + 3 + 3 + 4 = 13.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4], n = 4, left = 3, right = 4
#     Output: 6
#     Explanation: The given array is the same as example 1. We have the new array [1, 2, 3, 3, 4, 5, 6, 7, 9, 10]. The sum of the numbers from index le = 3 to ri = 4 is 3 + 3 = 6.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2,3,4], n = 4, left = 1, right = 10
#     Output: 50
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^3`
#   * `nums.length == n`
#   * `1 <= nums[i] <= 100`
#   * `1 <= left <= right <= n * (n + 1) / 2`
#
#
## @lc code=start
using LeetCode

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


text:
using Documenter, IntervalTrees

makedocs(
    modules = [IntervalTrees],
    format = :html,
    sitename = "IntervalTrees",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Types" => "man/types.md",
            "Creating IntervalTrees" => "man/insertion.md",
            "Dictionary Operations" => "man/dict.md",
            "Iterating over intervals" => "man/iteration.md",
            "Intersection" => "man/intersection.md"
        ]
    ],
    authors = "Daniel Jones, Kenta Sato, Ben J. Ward, The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/IntervalTrees.jl.git",
    julia = "0.7",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
abstract type AbstractModelContainer end

abstract type VariableType end
abstract type ConstraintType end
abstract type AuxVariableType end
abstract type ParameterType end
abstract type InitialConditionType end
abstract type ExpressionType end


text:
#---------------------------------------------------------------------#
#This function computes the Initial & Exact Solutions.
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function exact_solution(coord,Npoin,time,case,DFloat)

    #Set some constants
    xmin=0
    xmax=+2
    xc=0.5*(xmax+xmin)
    xl=xmax-xmin
    rc=0.125
    sigma=1/16; #best for right convergence rates
    # sigma=1/8
    # sigma=1/4; #fatter Gaussian to compare against 2D code

    #Initialize
    qe=zeros(DFloat,Npoin)

    #timec=time - floor(time)
    timec=time
    #Generate Grid Points
    for I=1:Npoin
        x=coord[I]
        if (case == 1) #Sin
            qe[I]=sin(pi*x) + 0.01;
        end
    end #I
    return (qe)
end


text:
using CoupledConditionalSMC
import SequentialMonteCarlo.SMCModel
using SMCExamples.FiniteFeynmanKac
import SMCExamples.Particles.Int64Particle
using RNGPool

# setRNGs(12345)

d = 2
n = 3
ffk = FiniteFeynmanKac.randomFiniteFK(d, n)

model = FiniteFeynmanKac.makeSMCModel(ffk)
densities = Vector{Float64}(undef, d^n)
for i = 1:length(densities)
  densities[i] = FiniteFeynmanKac.fullDensity(ffk, FiniteFeynmanKac.Int642Path(i, d, n))
end
densities ./= sum(densities)

lM = FiniteFeynmanKac.makelM(ffk)

function _getFreqs(model::SMCModel, lM::F, ccsmcio::CCSMCIO,
  algorithm::Symbol, m::Int64, d::Int64) where F<:Function
  ccsmcio.ref1 .= FiniteFeynmanKac.Int642Path(1, d, ccsmcio.n)
  ccsmcio.ref2 .= FiniteFeynmanKac.Int642Path(d^ccsmcio.n, d, ccsmcio.n)
  counts1 = zeros(Int64, d^ccsmcio.n)
  counts2 = zeros(Int64, d^ccsmcio.n)
  result1 = Vector{Float64}(undef, d^ccsmcio.n)
  result2 = Vector{Float64}(undef, d^ccsmcio.n)
  for i = 1:m
    ccXpf!(model, lM, ccsmcio, algorithm)
    counts1[FiniteFeynmanKac.Path2Int64(ccsmcio.ref1, d)] += 1
    counts2[FiniteFeynmanKac.Path2Int64(ccsmcio.ref2, d)] += 1
  end
  result1 .= counts1 ./ m
  result2 .= counts2 ./ m
  return result1, result2
end

nsamples = 2^22

# ccsmcio = CCSMCIO{model.particle, model.pScratch}(16, model.maxn)
ccsmcio = CCSMCIO{model.particle, model.pScratch}(4, model.maxn)
_getFreqs(model, error, ccsmcio, :AT, 4, d)
@time freqs1, freqs2 = _getFreqs(model, error, ccsmcio, :AT, nsamples, d)

println(sum(abs.(densities-freqs1)))
println(sum(abs.(densities-freqs2)))

_getFreqs(model, lM, ccsmcio, :BS, 4, d)
@time freqs1, freqs2 = _getFreqs(model, lM, ccsmcio, :BS, nsamples, d)

println(sum(abs.(densities-freqs1)))
println(sum(abs.(densities-freqs2)))

_getFreqs(model, lM, ccsmcio, :AS, 4, d)
@time freqs1, freqs2 = _getFreqs(model, lM, ccsmcio, :AS, nsamples, d)

println(sum(abs.(densities-freqs1)))
println(sum(abs.(densities-freqs2)))


text:
"""
    series(curves;
        linewidth=2,
        color=:lighttest,
        solid_color=nothing,
        labels=nothing,
        # scatter arguments, if any is set != nothing, a scatterplot is added
        marker=nothing,
        markersize=nothing,
        markercolor=automatic,
        strokecolor=nothing,
        strokewidth=nothing)

Curves can be:
* `AbstractVector{<: AbstractVector{<: Point2}}`: the native representation of a series as a vector of lines
* `AbstractMatrix`: each row represents y coordinates of the line, while `x` goes from `1:size(curves, 1)`
* `AbstractVector, AbstractMatrix`: the same as the above, but the first argument sets the x values for all lines
* `AbstractVector{<: Tuple{X<: AbstractVector, Y<: AbstractVector}}`: A vector of tuples, where each tuple contains a vector for the x and y coordinates

"""
@recipe(Series, curves) do scene
    Attributes(
        linewidth=2,
        color=:lighttest,
        solid_color=nothing,
        labels=nothing,

        marker=nothing,
        markersize=nothing,
        markercolor=automatic,
        strokecolor=nothing,
        strokewidth=nothing,
    )
end

function categorical_colors(cols::AbstractVector{<: Colorant}, categories::Integer)
    if length(cols) < categories
        error("Not enough colors for number of categories. Categories: $(categories), colors: $(length(cols))")
    end
    return to_colormap(cols)
end

function categorical_colors(cols::AbstractVector, categories::Integer)
    return categorical_colors(to_color.(cols), categories)
end

function categorical_colors(cs::Union{String, Symbol}, categories::Integer)
    cs_string = string(cs)
    if cs_string in all_gradient_names
        cols = PlotUtils.get_colorscheme(Symbol(cs_string)).colors
        categorical_colors(cols, categories)
    else
        error(
            """
            There is no color gradient named $cs.
            See `available_gradients()` for the list of available gradients,
            or look at http://makie.juliaplots.org/dev/generated/colors#Colormap-reference.
            """
        )
    end
end

replace_missing(x) = ismissing(x) ? NaN : x

function convert_arguments(T::Type{<: Series}, y::AbstractMatrix)
    convert_arguments(T, 1:size(y, 2), y)
end

function convert_arguments(::Type{<: Series}, x::AbstractVector, ys::AbstractMatrix)
    return (map(1:size(ys, 1)) do i
        Point2f.(replace_missing.(x), replace_missing.(view(ys, i, :)))
    end,)
end

function convert_arguments(::Type{<: Series}, arg::AbstractVector{<: Tuple{X, Y}}) where {X, Y}
    return (map(arg) do (x, y)
        Point2f.(replace_missing.(x), replace_missing.(y))
    end,)
end

function convert_arguments(T::Type{<: Series}, arg::Tuple{<:AbstractVector, <:AbstractVector})
    return convert_arguments(T, [arg])
end

function convert_arguments(::Type{<: Series}, arg::AbstractVector{<: AbstractVector{<:Point2}})
    return (map(arg) do points
        Point2f.(replace_missing.(first.(points)), replace_missing.(last.(points)))
    end,)
end

function plot!(plot::Series)
    @extract plot (curves, labels, linewidth, color, solid_color)
    sargs = [:marker, :markersize, :strokecolor, :strokewidth]
    scatter = Dict((f => plot[f] for f in sargs if !isnothing(plot[f][])))
    nseries = length(curves[])
    colors = lift(color, solid_color) do color, scolor
        if isnothing(scolor)
            return categorical_colors(color, nseries)
        else
            return scolor
        end
    end

    for i in 1:nseries
        label = @lift isnothing($labels) ? "series $(i)" : $labels[i]
        positions = @lift $curves[i]
        series_color = @lift $colors isa AbstractVector ? $colors[i] : $colors
        if !isempty(scatter)
            mcolor = plot.markercolor
            markercolor = @lift $mcolor == automatic ? $series_color : $mcolor
            scatterlines!(plot, positions;
                linewidth=linewidth, color=series_color, markercolor=series_color,
                label=label[], scatter...)
        else
            lines!(plot, positions; linewidth=linewidth, color=series_color, label=label)
        end
    end
end

function MakieLayout.get_plots(plot::Series)
    return plot.plots
end


text:
function get_raiberthopper(;
    timestep=0.05,
    gravity=[0.0; 0.0; -9.81],
    spring=0.0,
    damper=0.1,
    contact_foot=true,
    contact_body=true,
    T=Float64)
    #TODO: make customizable

    # Parameters
    leg_axis = [0.0; 0.0; 1.0]
    leg_length_nominal = 0.5
    body_radius = 0.1
    foot_radius = 0.05
    body_mass = 4.18 # from MuJoCo model
    foot_mass = 0.52 # from MuJoCo model

    # Links
    origin = Origin{Float64}()
    body = Sphere(body_radius, body_mass)
    foot = Sphere(foot_radius, foot_mass)
    links = [body, foot]

    # Joint Constraints
    joint_origin_body = JointConstraint(Floating(origin, body))
    joint_body_foot = JointConstraint(Prismatic(body, foot, leg_axis;
        parent_vertex=szeros(Float64, 3),
        child_vertex=szeros(Float64, 3),
        spring=spring,
        damper=damper))
    joints = [joint_origin_body, joint_body_foot]

    # Mechanism
    if contact_foot
         # Contact
        contact_normal = [0.0; 0.0; 1.0]
        friction_coefficient = 0.5

        # foot
        foot_contacts = contact_constraint(foot, contact_normal,
            friction_coefficient=friction_coefficient,
            contact_point=[0.0; 0.0; 0.0],
            offset=[0.0; 0.0;
            foot_radius])

        contacts = [foot_contacts]

        # body
        if contact_body
            body_contacts = contact_constraint(body, contact_normal,
                friction_coefficient=friction_coefficient,
                contact_point=[0.0; 0.0; 0.0],
                offset=[0.0; 0.0; body_radius])
            push!(contacts, body_contacts)
        end

        mech = Mechanism(origin, links, joints, contacts,
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    else
        mech = Mechanism(origin, links, joints,
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    end
    return mech
end

function initialize_raiberthopper!(mech::Mechanism{T,Nn,Ne,Nb};
    body_position=[0.0, 0.0, 0.05],
    leg_length=0.5,
    body_linear_velocity=zeros(3),
    body_angular_velocity=zeros(3)) where {T,Nn,Ne,Nb}

    pbody = mech.bodies[1]
    cbody = mech.bodies[2]
    joint2 = mech.joints[2]
    tra2 = joint2.translational

    # origin to body
    set_maximal_configurations!(mech.origin, pbody,
        Δx=[body_position[1:2]; leg_length + body_position[3]])
    set_maximal_velocities!(pbody,
        v=body_linear_velocity,
        ω=body_angular_velocity)

    # body to foot
    set_maximal_configurations!(pbody, cbody,
        Δx=[0.0; 0.0; -leg_length],
        Δq=RotX(0.0))
    set_maximal_velocities!(pbody, cbody,
        parent_vertex=tra2.vertices[1],
        child_vertex=tra2.vertices[2],
        Δv=zeros(3),
        Δω=zeros(3))
end


text:
"""

"""

export @namekw

"""
    @namekw func(expr)

    Alters a function call to include a name keyword based
        on the expression for the first argument.

    @namekw( func(expr) ) == func(expr; name = string( :(expr) ) )
"""
macro namekw(expr)
    ( isa(expr,Expr) && expr.head == :call ) ||
        error("@namekw must be used on a method call.")
    push!( expr.args, Expr( :kw, :name, string(expr.args[2]) ) )
    esc( expr )
end

"This macro is meant to replace the function of napari.utils.naming.magic_name"
function call_with_name_keyword(func, data, kwargs...)
    # We need to replace the function of magic_name
    local name_given = any([a.args[1] == :name for a in kwargs if isa(a,Expr)])
    local x = [esc(kw) for kw in kwargs]
    if name_given
        return :( $(func)( $( esc(data ) ); $(x...) ) )
    else
        local name = string(data)
        return :( $(func)( $( esc(data) ); name=$name, $(x...) ) )
    end
end

# We duplicate this function
"This macro is meant to replace the function of napari.utils.naming.magic_name"
function call_viewer_with_name_keyword(func, viewer, data, kwargs...)
    # We need to replace the function of magic_name
    local name_given = any([a.args[1] == :name for a in kwargs if isa(a,Expr)])
    local x = [esc(kw) for kw in kwargs]
    if name_given
        return :( $(func)( $( esc(viewer) ), $( esc(data) ); $(x...) ) )
    else
        local name = string(data)
        return :( $(func)( $( esc(viewer) ), $( esc(data) ); name=$name, $(x...) ) )
    end
end

#=
macro view_image(image, kwargs...)
    call_with_name_keyword( :( napari.view_image ),
                           image, kwargs...)
end
 =#

# see napari/view_layers.py
# see napari.layers.NAMES
for layer in layers
    local macro_name = Symbol("view_", layer)
    # local view = Expr( :. , :(napari), QuoteNode( macro_name ) )
    local view = QuoteNode( macro_name )
    eval( Expr( :export, Symbol('@', macro_name) ) )
    @eval begin
        """
            @view_{layer}(expr)

            Equivalent to napari.view_*(expr; name = string( :(expr) )).
            layer is one of image, points, labels, shapes, surface, vectors
            or tracks.
        """
        macro $macro_name($layer, kwargs...)
            call_with_name_keyword( $view , $layer, kwargs... )
        end
    end

    local macro_name = Symbol("add_", layer)
    local add = QuoteNode(macro_name)
    eval( Expr( :export, Symbol('@', macro_name) ) )
    @eval begin
        """
            @add_{layer}(viewer, expr)

            Equivalent to viewer.add_*(expr; name = string( :(expr) ))
            viewer is a napari.viewer.Viewer created by napari.Viewer()
            layer is one of image, points, labels, shapes, surface, vectors
            or tracks.
        """
        macro $macro_name(viewer, data, kwargs...)
            call_viewer_with_name_keyword( $add, viewer, data, kwargs... )
        end
    end
end


text:
module BadSeq

using BioSequences
using FASTX
using BioAlignments
using DataStructures

######## NamedSeq ################
struct NamedSeq{T <: BioSequence}
    name::String
    seq::T
end

Base.length(s::NamedSeq) = length(s.seq)
Base.setindex!(x::NamedSeq, i::Integer, v) = (x.seq[i] = v)
Base.getindex(x::NamedSeq, i::Integer) = x.seq[i]

function Base.print(io::IO, s::NamedSeq)
    print(io, '>')
    println(io, s.name)
    print(io, s.seq)
end

function Base.read(reader::FASTA.Reader, ::Type{NamedSeq})
    record = read!(reader, FASTA.Record())
    header = FASTA.identifier(record)
    FASTA.hasdescription(record) && (header *= (" " * FASTA.description(record)))
    seq = FASTA.sequence(record)
    return NamedSeq(header, seq)
end

function read_fasta(io::IO)
    io = io isa FASTA.Reader ? io : FASTA.Reader(io)
    result = NamedSeq[]
    while !eof(io)
        push!(result, read(io, NamedSeq))
    end
    result
end

macro seq_str(str, name)
    :(NamedSeq($(name), LongDNASeq($(str))))
end

"MAFFT align sequences. By default prints a log to /tmp"
function mafft(seqs, stderr=tempname("/tmp"))
    inpath = tempname("/tmp")

    open(inpath, "w") do file
        for seq in seqs
            println(file, seq)
        end
    end
    out = IOBuffer()
    if stderr isa String
        open(stderr, "w") do file
            run(pipeline(pipeline(`mafft $inpath`; stderr=file), out))
        end
    else
        run(pipeline(pipeline(`mafft $inpath`; stderr=stderr), out))
    end
    rm(inpath)
    return read_fasta(seekstart(out))
end

"Trim '-' from both ends of an iterable of Seqs"
function trim(seqs)
    strings = [String(s.seq) for s in seqs]
    leftpad  = maximum([length(s) - length(lstrip(s, ['-'])) for s in strings])
    rightpad = maximum([length(s) - length(rstrip(s, ['-'])) for s in strings])
    return [NamedSeq(s.name, typeof(s.seq)(st[leftpad+1:end-rightpad])) for (s, st) in zip(seqs, strings)]
end
#########################

######### DNA - AA alignment ############

function findstarts(seq::LongSequence{DNAAlphabet{2}})
	starts = Int[]
	query = dna"ATG"
	p = 1
	loc = findfirst(query, seq, 1)
	while loc !== nothing
		push!(starts, first(loc))
		p = last(loc) + 1
		loc = findfirst(query, seq, p)
	end
	return starts
end

function orflength(seq::LongSequence{DNAAlphabet{2}}, start)
	it = each(DNACodon, seq[start:end], 3)
	stop = 0
	for (pos, fw, rv) in it
		stop = pos-1
		if (fw == mer"TAG") | (fw == mer"TGA") | (fw == mer"TAA")
			break
		end
	end
	return stop
end

function likely_orfs_forward(seq::LongSequence{DNAAlphabet{2}}, len=1)
	starts = findstarts(seq)
	likely_orfs = UnitRange{Int}[]
	cutoff = round(Int, 0.8 * len)
	for start in starts
		length(seq) - start < cutoff && break
		orflen = orflength(seq, start)
		orflen >= cutoff && push!(likely_orfs, start:start+orflen-1)
	end
	return likely_orfs
end

function best_alignment_forward(seq::LongSequence{DNAAlphabet{2}}, gene::AminoAcidSeq)
	likely_orfs = likely_orfs_forward(seq, 3*length(gene))
	isempty(likely_orfs) && return (1:0, typemin(Int))

	scoremodel = AffineGapScoreModel(BLOSUM62, gap_open=-5, gap_extend=-1)
	alignments = []

	for orf in likely_orfs
		orfgene = translate(seq[orf])
		push!(alignments, pairalign(GlobalAlignment(), orfgene, gene, scoremodel))
	end

	bestindex = argmax(map(score, alignments))
	return (likely_orfs[bestindex], score(alignments[bestindex]))
end

function findbest(gene::AminoAcidSeq, seq::LongSequence{<:NucleicAcidAlphabet})
	seq_ = convert(LongSequence{DNAAlphabet{2}}, seq)
	orf_fw, score_fw = best_alignment_forward(seq_, gene)
	orf_rv, score_rv = best_alignment_forward(reverse_complement(seq_), gene)

    reversed = score_rv > score_fw
    return reversed, ifelse(reversed, orf_rv, orf_fw)
end

###########
"Get all maximum-length orfs separated by a stop codon, if range satisfies f"
function get_orfs(seq::BioSequence{<:NucleicAcidAlphabet})
    orfs = UnitRange{Int}[]
    for frame in 1:3
        lastindex = length(seq) - (length(seq) - frame + 1) % 3
        aas = translate(seq[frame:lastindex])
        start = frame
        local stop = 0
        for i in eachindex(aas)
            index = 3*(i-1) + frame
            if aas[i] == AA_Term
                stop = index - 1
                stop > start && push!(orfs, start:stop)
                start = index + 3
            end
        end
        # If no stop was detected at end.
        start < lastindex && push!(orfs, start:lastindex)
    end
    return orfs
end

function get_met_orfs(seq::BioSequence{<:NucleicAcidAlphabet})
    orfs = UnitRange{Int}[]
    start = typemax(Int)
    stop = 0
    for frame in 1:3
        lastindex = length(seq) - (length(seq) - frame + 1) % 3
        aas = translate(seq[frame:lastindex])
        for i in eachindex(aas)
            index = 3*(i-1) + frame
            if aas[i] == AA_Term
                stop = index - 1
                stop > start && push!(orfs, start:stop)
                start = typemax(Int)
            elseif (start == typemax(Int)) & (aas[i] == AA_M)
                start = index
            end
        end
        stop > start && push!(orfs, start:stop)
    end
    return orfs
end

longest(xs) = maximum(((length(i), i) for i in xs))[2]

"Get (is_rc, largest ORF) in the sequence. `met` is whether it must start with Met"
function largest_orf(sequence::BioSequence{<:NucleicAcidAlphabet}, met=false)
    seq = copy(sequence)
    f = ifelse(met, get_met_orfs, get_orfs)
    fw = longest(f(seq))
    rv = longest(f(reverse_complement!(seq)))
    reversed = length(rv) > length(fw)
    return reversed, ifelse(reversed, rv, fw)
end

function largest_aa(sequence::BioSequence{<:NucleicAcidAlphabet}, met=false)
    seq = copy(sequence)
    isrv, orfrange = largest_orf(sequence, met)
    return translate((isrv ? reverse_complement(sequence) : sequence)[orfrange])
end

end # module

max_stars_repo_name: jakobnissen/BadBio


text:
"""
    ConstantDecrementRule

Fixed temperature decrement rule. The temperature is decrease by a constant
factor after each Markov chain.

Parameters
==========
factor: the constant factor multiplying the temperature.
"""
@with_kw struct ConstantDecrementRule{T}
    factor::T = 0.9
end

function (decrement_rule::ConstantDecrementRule)(state::AnnealingState)
    return decrement_rule.factor * state.temperature
end

"""
    AVLDecrementRule

Aarts and van Laarhoven temperature decrement rule.

Parameters
==========
distance_parameter

Reference
=========
Varanelli, eq. 2.42.
"""
@with_kw struct AVLDecrementRule{T}
    distance_parameter::T = 0.085
end

function (decrement_rule::AVLDecrementRule)(state::AnnealingState)
    T = state.temperature
    σ = std(state.energies)
    δ = decrement_rule.distance_parameter

    return T/(1 + T*log(1 + δ)/3σ)
end


text:
using Test
import Pluto: Configuration, Notebook, ServerSession, ClientSession, update_run!, Cell, WorkspaceManager, SessionActions
import Pluto.Configuration: Options, EvaluationOptions
import Distributed
using Pluto.WorkspaceManager: poll
import Pkg

@testset "Reload from file" begin

    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = false
    🍭.options.server.auto_reload_from_file = true

    timeout_between_tests = 🍭.options.server.auto_reload_from_file_cooldown * 1.5

    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    notebook = SessionActions.new(🍭; run_async=false)
    fakeclient.connected_notebook = notebook

    ###
    sleep(timeout_between_tests)

    nb1 = Notebook([
        Cell("x = 123"),
        Cell("x + x"),
        Cell("rand()"),
    ])
    file1 = sprint(Pluto.save_notebook, nb1)

    write(notebook.path, file1)

    @test poll(30) do
        length(notebook.cells) == 3
    end
    @test poll(5) do
        notebook.cells[1].output.body == "123"
    end
    @test poll(5) do
        all(c -> !c.running, notebook.cells)
    end

    @test notebook.cells[2].output.body == "246"
    @test notebook.cells[3].errored == false

    original_rand_output = notebook.cells[3].output.body

    ###
    sleep(timeout_between_tests)

    nb2 = Notebook(reverse(nb1.cells))

    file2 = sprint(Pluto.save_notebook, nb2)
    write(notebook.path, file2)

    @test poll(10) do
        notebook.cells[3].output.body == "123"
    end

    # notebook order reversed, but cell should not re-run
    @test original_rand_output == notebook.cells[1].output.body

    ###
    sleep(timeout_between_tests)

    file3 = replace(file1, "123" => "6")
    write(notebook.path, file3)

    @test poll(10) do
        notebook.cells[1].output.body == "6"
    end
    @test poll(5) do
        all(c -> !c.running, notebook.cells)
    end
    @test notebook.cells[2].output.body == "12"

    # notebook order reversed again, but cell should not re-run
    @test original_rand_output == notebook.cells[3].output.body

    ###
    sleep(timeout_between_tests)

    file4 = read(joinpath(@__DIR__, "packages", "simple_stdlib_import.jl"), String)
    write(notebook.path, file4)

    @test poll(10) do
        notebook.cells[2].output.body == "false"
    end
    @test length(notebook.cells) == 2
    @test notebook.nbpkg_restart_required_msg !== nothing
end


text:
using CmdStan, Distributions, Random

Random.seed!(103)
n_subj = 100
n_obs = 100
μs = rand(Normal(0, 1), n_subj)
data = map(μ->rand(Normal(μ, 1), n_obs), μs)
data = reduce(hcat, data)'

ProjDir = @__DIR__
cd(ProjDir)

CmdStanGaussian = "
data {
  int<lower=0> n_subj;
  int<lower=0> n_obs;
  matrix[n_subj,n_obs] y;
}
parameters {
  real mumu;
  real<lower=0> musigma;
  real<lower=0> sigma;
  vector[n_subj] mus;
}
model {
  mumu ~ normal(0, 1);
  musigma ~ gamma(1, 1);
  mus ~ normal(mumu, musigma);
  sigma ~ gamma(1, 1);
  for(i in 1:n_subj){
    y[i,:] ~ normal(mus[i], sigma);
  }
}
"

CmdStanConfig = Stanmodel(
  name = "CmdStanGaussian", model=CmdStanGaussian, nchains=4, output_format=:mcmcchains,
  Sample(
    num_samples=1000, num_warmup=1000, adapt=CmdStan.Adapt(delta=0.65),
    save_warmup=false
  )
)

observations = Dict("n_obs"=>n_obs, "n_subj"=>n_subj, "y"=>data)

@elapsed rc, samples, cnames = stan(CmdStanConfig, observations, ProjDir, CmdStanDir=CMDSTAN_HOME)


text:
@testset "396.rotate-function.jl" begin
    @test max_rotate_function([4, 3, 2, 6]) == 26
    @test max_rotate_function([40, 3, 44, 39, 49, 6, 29, 14, 8, 50, 3, 18, 12, 33, 40, 17, 22, 8, 25, 38, 2, 15, 25, 27, 30, 9, 36, 40, 3, 40]) == 11852
end


text:
@testset "Step tests" begin
    cfg = Config()
    ccon = rand() * 0.49 + 0.51
    c = const_controller(cfg; c=ccon)
    m = Model(cfg, c)
    nin = rand(5:10)
    nout = rand(5:10)
    random_init!(m, nin, nout)

    @testset "Update cell state" begin
        Neurodevo.update_cell_states!(m)
        for cell in m.cells
            @test all(cell.state .== ccon)
            @test cell.output == ccon
            @test cell.input == 0.0
        end
    end

    @testset "Update conn state" begin
        Neurodevo.update_conn_states!(m)
        for cell in m.cells
            @test cell.input <= 1.0
            @test cell.input >= -1.0
            for conn in cell.conns
                @test all(conn.state .== ccon)
            end
        end
    end

    @testset "Update params" begin
        Neurodevo.update_params!(m)
        for cell in m.cells
            @test all(cell.params .== ccon)
            for conn in cell.conns
                @test all(conn.params .== ccon)
            end
        end
    end

    @testset "Add cells" begin
        ncells = length(m.cells)
        Neurodevo.add_cells!(m)
        @test length(m.cells) == 2 * ncells
        for i in (ncells+1):length(m.cells)
            @test all(m.cells[i].params .== ccon)
        end
    end

    @testset "Remove cells" begin
        Neurodevo.remove_cells!(m)
        @test length(m.cells) == nin + nout
        for i in eachindex(m.cells)
            @test m.cells[i].interface
        end
    end

    @testset "Connect cells" begin
        Neurodevo.connect_cells!(m)
        nconns = 0
        for cell in m.cells
            nconns += length(cell.conns)
        end
        @test nconns <= m.cfg["conns_max"]
    end

    @testset "Disconnect cells" begin
        Neurodevo.disconnect_cells!(m)
        for cell in m.cells
            @test length(cell.conns) == 0
        end
    end
end