text:
using SafeTestsets

@safetestset "Mine Count" begin
    using Test, Random, Minesweeper
    Random.seed!(6564)
    n_mines = 3
    game = Game(dims=(5,5), n_mines=n_mines)
    @test n_mines == count(x->x.has_mine, game.cells)
end

@safetestset "Neighbors" begin
    using Test, Minesweeper

    game = Game(dims=(5,5), n_mines=3)
    neighbors = get_neighbor_indices(game, 1, 1)
    @test length(neighbors) == 3

    neighbors = get_neighbor_indices(game, 2, 2)
    @test length(neighbors) == 8
    v = [-1,0,1]
    for Δx in v, Δy in v
        if (Δx == 0) & (Δy == 0)
            @test CartesianIndex(2 + Δx, 2 + Δy) ∉ neighbors
        else
            @test CartesianIndex(2 + Δx, 2 + Δy) ∈ neighbors
        end
    end
end

@safetestset "reveal zeros" begin
    using Test, Random, Minesweeper

    Random.seed!(6564)
    game = Game(dims=(5,5), n_mines=3)
    reveal_zeros!(game, 3, 2)
    @test count(x->x.revealed, game.cells) == 0

    Random.seed!(6564)
    game = Game(dims=(5,5), n_mines=3)
    reveal_zeros!(game, 1, 4)
    @test count(x->x.revealed, game.cells) == 6
    revealed_indices = ((1,3),(1,4),(1,5),(2,3),(2,4),(2,5))
    for i in revealed_indices
        @test game.cells[i...].revealed
    end
    revealed = filter(x->x.revealed, game.cells)
    @test !any(x->x.has_mine, revealed)
end

# @safetestset "gui tests" begin
#     using Test, Random, Minesweeper
#     Random.seed!(59808)
#     import_gui()
#     game = Game(dims=(10,10), n_mines=10)
#     gui = generate_gui(game)
#     cell = game.cells[1]
#     flag_cell!(cell, game, gui)
#     @test cell.flagged
#     @test game.mines_flagged == 1

#     idx = findfirst(x->x.mine_count ==1, game.cells)
#     select_cell!(game, idx)
#     @test game.cells[idx].revealed

#     idx = findfirst(x->x.has_mine, game.cells)
#     select_cell!(game, idx)
#     @test game.cells[idx].revealed
#     @test game.mine_detonated
#     @test game_over(game)

#     game = Game(dims=(10,10), n_mines=10)

#     for cell in game.cells
#         if cell.has_mine
#             cell.flagged = true
#         else
#             cell.revealed = true
#         end
#     end
#     @test game_over(game)
# end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "csrsm" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        @testset "csr" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSR(sparse(A))
            info = CUSPARSE.sm_analysis('N','U',d_A,'O')
            d_Y = CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = to_host(d_y)
            y = A\X
            @test y ≈ h_y
            d_X = CudaArray(rand(elty,n,n))
            @test_throws DimensionMismatch CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.sm_analysis('T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end
        @testset "csc" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSC(sparse(A))
            info = CUSPARSE.sm_analysis('N','U',d_A,'O')
            d_Y = CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = to_host(d_y)
            y = A\X
            @test y ≈ h_y
            d_X = CudaArray(rand(elty,n,n))
            @test_throws DimensionMismatch CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.sm_analysis('T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end
    end
end

ext: jl


text:
module MLBasedESC

using Base: front, tail
using LinearAlgebra
using Revise
using Printf

import Flux
import Flux.Zygote
using Flux: Chain, Dense, elu

import DiffEqFlux
using DiffEqFlux: FastChain, FastDense

using DynamicPolynomials

using Distributions
import ReverseDiff

using OrdinaryDiffEq, SciMLBase

include("fluxhelpers.jl")
include("idapbc.jl")
include("soslayers.jl")
include("constants.jl")
include("systems.jl")
include("neuralpbc.jl")
include("optimize.jl")

end


text:
if basename(pwd()) == "aoc"
    cd("2017/23")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module ExampleExperiment

# Learning Updates

struct QLearning
    γ::Float64
end

struct TD
    γ::Float64
end

struct NStepQLearning
    γ::Float64
    n::Int
end

struct NStepTDLearning
    γ::Float64
    n::Int
end

struct QLambda
    γ::Float64
    λ::Float64
end

struct TDLambda
    γ::Float64
    λ::Float64
end

macro create_lu_build_trait(lu_type, trait)
    fn = :lu_build_trait
    quote
        function $(esc(fn))(::Type{$lu_type})
            $trait()
        end
    end
end

lu_build_trait(lu) = @error "$(lu) not supported! Please implement `construct_rnn_layer`"

# Fun with constructors:
struct LUBuildVanilla end

@create_lu_build_trait QLearning LUBuildVanilla
@create_lu_build_trait TD LUBuildVanilla

struct LUBuildNStep end
# If there are a lot use a for loop
for lu in [NStepQLearning, NStepTDLearning]
    @eval begin
        @create_lu_build_trait $(lu) LUBuildNStep
    end
end

struct LUBuildTrace end

@create_lu_build_trait QLambda LUBuildTrace
@create_lu_build_trait TDLambda LUBuildTrace

"""
    construct_learning_update(config)

Construct the learning update. Types of constructors:
- LUBuildVanilla
- LUBuildNStep
- LUBuildTrace
"""
function construct_learning_update(config)
    lu = if isdefined(ExampleExperiment, Symbol(config["learning_update"]))
        getproperty(ExampleExperiment, Symbol(config["learning_update"]))
    else
        @error "Oh No! Learning Update not defined"
    end

    construct_learning_update(lu, config)
end

construct_learning_update(lu, config) =
    construct_learning_update(lu_build_trait(lu), lu, config)

function construct_learning_update(::LUBuildVanilla, lu, config)
    lu(config["gamma"])
end

function construct_learning_update(::LUBuildNStep, lu, config)
    lu(config["gamma"], config["nstep"])
end

function construct_learning_update(::LUBuildTrace, lu, config)
    lu(config["gamma"], config["lambda"])
end

# Macro Tools are for doing better configs.
include("macro_tools.jl")

end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/../" "${BASH_SOURCE[0]}" "$@" -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

const feature_types = (
    two_vertical = (1, 2),
    two_horizontal = (2, 1),
    three_horizontal = (3, 1),
    three_vertical = (1, 3),
    four = (2, 2)
    )

abstract type HaarFeatureAbstractType end

"""
    mutable struct HaarLikeObject{I <: Integer, F <: AbstractFloat}

        Struct representing a Haar-like feature.

    feature_type::Tuple{I, I}
    position::Tuple{I, I}
    top_left::Tuple{I, I}
    bottom_right::Tuple{I, I}
    width::I
    height::I
    threshold::I
    polarity::I
    weight::F
"""
mutable struct HaarLikeObject{I <: Integer, F <: AbstractFloat} <: HaarFeatureAbstractType
    #parametric struct to store the ints and floats efficiently
    feature_type::Tuple{I, I}
    position::Tuple{I, I}
    top_left::Tuple{I, I}
    bottom_right::Tuple{I, I}
    width::I
    height::I
    threshold::I
    polarity::I
    weight::F
end # end structure

"""
    HaarLikeObject(
        feature_type::Tuple{Integer, Integer},
        position::Tuple{Integer, Integer},
        width::Integer,
        height::Integer,
        threshold::Integer,
        polarity::Integer
    ) -> HaarLikeObject
"""
function HaarLikeObject(
    feature_type::Tuple{Integer, Integer},
    position::Tuple{Integer, Integer},
    width::Integer,
    height::Integer,
    threshold::Integer,
    polarity::Integer
)

    # make sure that everything is of the same size
    p₁, p₂ = position
    f₁, f₂ = feature_type
    p₁, p₂, f₁, f₂, width, height, threshold, polarity =
        promote(p₁, p₂, f₁, f₂, width, height, threshold, polarity)
    position = (p₁, p₂)
    feature_type = (f₁, f₂)
    top_left = position
    bottom_right = (first(position) + width, last(position) + height)
    weight = float(one(p₁)) #to make a float of the same size

    HaarLikeObject(feature_type, position, top_left, bottom_right, width, height, threshold, polarity, weight)
end

"""
    get_score(feature::HaarLikeObject, int_img::Array) -> Tuple{Number, Number}

Get score for given integral image array.  This is the feature cascade.

# Arguments

- `feature::HaarLikeObject`: given Haar-like feature (parameterised replacement of Python's `self`)
- `int_img::AbstractArray`: Integral image array

# Returns

- `score::Number`: Score for given feature
"""
function get_score(feature::HaarLikeObject{I, F}, int_img::Array) where {I, F}
    score = zero(I)
    faceness = zero(I)
    _2f = F(2)
    _3f = F(3)
    _half = F(0.5)
    _one_third = F(1.0 / 3.0)

    if feature.feature_type == feature_types.two_vertical
        _first = sum_region(int_img, feature.top_left, (first(feature.top_left) + feature.width, round(I, last(feature.top_left) + feature.height / 2)))
        second = sum_region(int_img, (first(feature.top_left), round(I, last(feature.top_left) + feature.height / 2)), feature.bottom_right)
        score = _first - second
        faceness = I(1)
    elseif feature.feature_type == feature_types.two_horizontal
        _first = sum_region(int_img, feature.top_left, (round(I, first(feature.top_left) + feature.width / 2), last(feature.top_left) + feature.height))
        second = sum_region(int_img, (round(I, first(feature.top_left) + feature.width / 2), last(feature.top_left)), feature.bottom_right)
        score = _first - second
        faceness = I(2)
    elseif feature.feature_type == feature_types.three_horizontal
        _first = sum_region(int_img, feature.top_left, (round(I, first(feature.top_left) + feature.width / 3), last(feature.top_left) + feature.height))
        second = sum_region(int_img, (round(I, first(feature.top_left) + feature.width / 3), last(feature.top_left)), (round(I, first(feature.top_left) + 2 * feature.width / 3), last(feature.top_left) + feature.height))
        third = sum_region(int_img, (round(I, first(feature.top_left) + 2 * feature.width / 3), last(feature.top_left)), feature.bottom_right)
        score = _first - second + third
        faceness = I(3)
    elseif feature.feature_type == feature_types.three_vertical
        _first = sum_region(int_img, feature.top_left, (first(feature.bottom_right), round(I, last(feature.top_left) + feature.height / 3)))
        second = sum_region(int_img, (first(feature.top_left), round(I, last(feature.top_left) + feature.height / 3)), (first(feature.bottom_right), round(I, last(feature.top_left) + 2 * feature.height / 3)))
        third = sum_region(int_img, (first(feature.top_left), round(I, last(feature.top_left) + 2 * feature.height / 3)), feature.bottom_right)
        score = _first - second + third
        faceness  = I(4)
    elseif feature.feature_type == feature_types.four
        # top left area
        _first = sum_region(int_img, feature.top_left, (round(I, first(feature.top_left) + feature.width / 2), round(I, last(feature.top_left) + feature.height / 2)))
        # top right area
        second = sum_region(int_img, (round(I, first(feature.top_left) + feature.width / 2), last(feature.top_left)), (first(feature.bottom_right), round(I, last(feature.top_left) + feature.height / 2)))
        # bottom left area
        third = sum_region(int_img, (first(feature.top_left), round(I, last(feature.top_left) + feature.height / 2)), (round(I, first(feature.top_left) + feature.width / 2), last(feature.bottom_right)))
        # bottom right area
        fourth = sum_region(int_img, (round(I, first(feature.top_left) + feature.width / 2), round(I, last(feature.top_left) + feature.height / 2)), feature.bottom_right)
        score = _first - second - third + fourth
        faceness = I(5)
    end

    return score, faceness
end

"""
    get_vote(feature::HaarLikeObject, int_img::AbstractArray) -> Integer

Get vote of this feature for given integral image.

# Arguments

- `feature::HaarLikeObject`: given Haar-like feature (parameterised replacement of Python's `self`)
- `int_img::AbstractArray`: Integral image array [type: Abstract Array]

# Returns

 - `vote::Integer`:
    1       ⟺ this feature votes positively
    -1      otherwise
"""
function get_vote(feature::HaarLikeObject, int_img::AbstractArray)
    score = first(get_score(feature, int_img)) # we only care about score here, not faceness
    # return (feature.weight * score) < (feature.polarity * feature.threshold) ? one(Int8) : -one(Int8)
    # return feature.weight * (score < feature.polarity * feature.threshold ? one(Int8) : -one(Int8))
    return score < feature.polarity * feature.threshold ? feature.weight : -feature.weight
    # self.weight * (1 if score < self.polarity * self.threshold else -1)
end


text:
using GraphPlot, Compose #, ImageMagick, LightGraphs # RoME, DistributedFactorGraphs, Gadfly
# include(joinpath(pkgdir(RoME), "test", "runtests.jl")) # February 2021: some tests fail here
# include(joinpath(pkgdir(Gadfly), "test", "runtests.jl")) # February 2021: some tests fail here
include(joinpath(pkgdir(GraphPlot), "test", "runtests.jl"))
# include(joinpath(pkgdir(DistributedFactorGraphs), "test", "runtests.jl")) # February 2021: CI already errors?
# include(joinpath(pkgdir(Compose), "test", "runtests.jl")) # February 2021: some tests fail here


text:
"""
    compact_mnm_hyper(α::AbstractArray{Float32,1})

[Dirichlet Distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution)
"""
struct compact_mnm_hyper <: distribution_hyper_params
    α::AbstractArray{Float32,1}
end

mutable struct compact_mnm_sufficient_statistics <: sufficient_statistics
    N::Float32
    points_sum::AbstractArray{Float32,1}
end

function calc_posterior(prior:: compact_mnm_hyper, suff_statistics::Vector{Tuple{sufficient_statistics,Number}})
    if length(suff_statistics) == 0
        return prior
    end
    points_num = sum([post_kernel(x[2],global_time)*x[1].N for x in suff_statistics])
    if points_num == 0
        return prior
    end

    suff_sum = sum([post_kernel(x[2],global_time)*x[1].points_sum for x in suff_statistics],dims=1)[1]
    return compact_mnm_hyper(prior.α + suff_sum)
end

function sample_distribution(hyperparams::compact_mnm_hyper)
    return compact_mnm_dist(log.(rand(Dirichlet(Float64.(hyperparams.α)))))
end

function create_sufficient_statistics(hyper::compact_mnm_hyper,posterior::compact_mnm_hyper,points::AbstractArray{Float32,2}, pts_to_group = 0)
    if length(points) == 0
        return compact_mnm_sufficient_statistics(size(points,2),zeros(Int64,length(hyper.α)))
    end
    pt_count = counts(Int.(points),length(hyper.α))
    return compact_mnm_sufficient_statistics(size(points,2),pt_count)
end

function log_marginal_likelihood(hyper::compact_mnm_hyper, posterior_hyper::compact_mnm_hyper, suff_stats::Vector{Tuple{sufficient_statistics,Number}})
    D = length(suff_stats[1][1].points_sum)
    logpi = log(pi)
    val = logabsgamma(sum(hyper.α))[1] -logabsgamma(sum(posterior_hyper.α))[1] + sum((x-> logabsgamma(x)[1]).(posterior_hyper.α) - (x-> logabsgamma(x)[1]).(hyper.α))
    return val
end

function aggregate_suff_stats(suff_l::compact_mnm_sufficient_statistics, suff_r::compact_mnm_sufficient_statistics)
    return compact_mnm_sufficient_statistics(suff_l.N+suff_r.N, suff_l.points_sum + suff_r.points_sum)
end

function posterior_predictive!(r::AbstractArray,x::AbstractMatrix,posterior_hyper::compact_mnm_hyper)
    if length(posterior_hyper.α) == 0
        return
    end
    post_sum = float(sum(posterior_hyper.α))
    @inbounds for i in eachindex(r)
        r[i] = 0.0
        @inbounds for j=1:size(x,1)
            if x[j,i] > 0
                r[i] += log(posterior_hyper.α[Int64.(x[j,i])]/post_sum)
            end
        end
    end
end


text:
__precompile__()

module GtkBuilderAid
using Gtk
using Cairo

export @GtkBuilderAid, @GtkFunctionTable, start_application

include("connect_signals.jl")
include("aidbuild.jl")
include("misc.jl")

end # module


text:
# # Image classification
#
# *See [`ImageClassification`](#) for complete documentation of its arguments.*
#
# Let's explore what you can do with the [`LearningMethod`](#) interface implemented. We're using
# [DLDatasets.jl](https://github.com/lorenzoh/DLDatasets.jl) to access *ImageNette*, a small
# image classification dataset. Install the package using
#
#   `]add https://github.com/lorenzoh/DLDatasets.jl`
# {cell=main style="display:none" result=false}

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
ENV["CI"] = "false"

# {cell=main output=false}

using DLPipelines
using DLDatasets
using LearnBase: getobs

# With the packages imported, loading the dataset is straightforward:
# {cell=main}

dataset = DLDatasets.loaddataset(ImageNette, "v2_160px")

# As you can see, every observation consists of an image and a category label:
# {cell=main}
image, category = getobs(dataset, 4000)
image
# {cell=main, result = false style="display:none;"}
@show category

# We can also retrieve the 10 different categories from the dataset's metadata:
# {cell=main}

categories = DLDatasets.metadata(ImageNette).labels

# ## Image classification `LearningMethod`
#
# With the dataset ready, we can create an instance of [`ImageClassification`](#).
#
# Note that you could use any image classification dataset, as long as `getobs(ds, idx)`
# returns an image and a category.
# {cell=main}

method = ImageClassification(categories, sz = (128, 128))

# We can now use this `method` with a [`Context`](#) to transform the data.
# The image is encoded as a normalized 3D-array:
# {cell=main}

x = encodeinput(method, Training(), image)
summary(x)

# And the category is one-hot encoded:
# {cell=main}

y = encodetarget(method, Training(), category)

# You can also use [`MethodDataset`](#) to create a wrapper around your existing dataset
# that directly returns encoded observations:
# {cell=main}

methoddataset = MethodDataset(dataset, method, Training())
x, y = getobs(methoddataset, 1)
summary.((x, y))


text:
#///////////////////////////////////////
#// File Name: crowd_nav_controller_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/04/08
#// Description: Test code for src/crowd_nav_controller.jl
#///////////////////////////////////////

using LinearAlgebra
using PyCall

@testset "Crowd Nav Controller Test" begin
    Cep = Matrix(0.5I, 2, 2);                                                           # quadratic position cost matrix
    Cu = Matrix(0.2I, 2, 2);                                                            # quadratic control cost matrix
    β_pos = 0.1;                                                                        # relative weight between instant and terminal pos cost
    α_col = 100.0;                                                                      # scale parameter for exponential collision cost
    β_col = 0.1;                                                                        # relative weight between instant and terminal pos cost
    λ_col = 0.2;                                                                        # bandwidth parameter for exponential collision cost
    σ_risk = 0.0;                                                                       # risk-sensitiveness parameter

    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    sim_param = SimulationParameter(0.01, 0.4, 12, 10, cost_param);

    tcalc = 0.1;
    dtr = 0.4;
    target_speed = 0.7;
    model_dir = normpath(joinpath(@__DIR__,
                                  "../CrowdNav/crowd_nav/data/output_om_sarl_radius_0.4"))
    cnt_param = CrowdNavControlParameter(model_dir,
                                         "env.config", "policy.config", "sarl",
                                         [0.0, 4.0], tcalc, dtr, target_speed);

    crowd_nav_controller = CrowdNavController(sim_param, cnt_param);
    @test !crowd_nav_controller.rl_robot.visible
    @test crowd_nav_controller.rl_robot.policy.time_step == 0.4;
    @test crowd_nav_controller.rl_robot.v_pref == 0.7;
    @test crowd_nav_controller.rl_robot.time_step == 0.4;
    @test crowd_nav_controller.rl_robot.gx == 0.0;
    @test crowd_nav_controller.rl_robot.gy == 4.0;

    e_init = RobotState([0.0, -4.0, 0.0, 0.0]);
    u_vel = get_action!(crowd_nav_controller, e_init,
                        Dict{String, Vector{Float64}}());
    @test u_vel ≈ [0.0, 8.0]./norm([0.0, 8.0])*crowd_nav_controller.rl_robot.v_pref

    e_init = RobotState([0.0, 0.0, 0.0, 0.0]);
    ado_state_dict = Dict{String, Vector{Float64}}();
    ado_state_dict["PEDESTRIAN/1"] = [5.0, 4.0, -1.0, -3.0];
    ado_state_dict["PEDESTRIAN/2"] = [-3.0, 0.0, 1.0, 1.0];

    u_vel = get_action!(crowd_nav_controller, e_init, ado_state_dict)

    ~, u_acc = crowdnav_control_update!(crowd_nav_controller, e_init, ado_state_dict)
    #@test all(u_vel.*0.4 .≈ get_velocity(e_init).*0.4 .+ 0.5*u_acc.*(0.4)^2)
    @test all(u_vel .≈ get_velocity(e_init) .+ u_acc.*0.4)
end


text:
using Distributions, Random

#Simulates with the cool off policy
function simCoolOff(N,σ,γ,c,tMax; seed = 1974)
    Random.seed!(seed)
    state = [0 for _ in 1:N] #actual state of arms
    allowedToSwitchTo = [true for _ in 1:N] #state of belief state of arms (true => passed the threshold)
    currentArm = 1
    t = 0.0
    reward = 0.0
    λ, μ = 1.0, 1/γ - 1
    rλ() = rand(Exponential(1/λ))
    rμ() = rand(Exponential(1/μ))
    nextArm(i) = (i%N)+1

    #2N event types in this discrete event simulation
    timeJumps = [rλ() for _ in 1:N] #first N events are for state change
    append!(timeJumps,fill(σ,N)) #initilize N events for "ok to switch to arm"
    timeJumps[currentArm + N] = Inf

    #assumes the currentArm is bad (in 0) and searches for a better arm afterwards
    function searchForNextArm()
        na = nextArm(currentArm)
        for _ in 1:N-1 #loop untill either (1) found a good arm (2) no where to switch to
            if allowedToSwitchTo[na] #if allowed to switch - lets do it
                timeJumps[currentArm + N] = t + σ #leaving arm so set the belief update
                allowedToSwitchTo[currentArm] = false
                currentArm = na
                reward -= c
               # @show timeJumps
                if state[currentArm] == 1 #we are happy and leave
                    break
                end
            end
            na = nextArm(na)
        end
        timeJumps[currentArm + N] = Inf
    end

    while t<tMax
       tLast = t
       t,i = findmin(timeJumps)
       reward += state[currentArm]*(t-tLast)
        if i >= N+1 #pass the omega threshold event
            allowedToSwitchTo[i-N] = true
            if state[currentArm] == 0
                searchForNextArm()
            else
                timeJumps[i] = Inf
            end
        else #arm i is changing state event
           if state[i] == 0
              state[i] = 1
              timeJumps[i] = t + rμ()
           else #state was 1
                state[i] = 0
                timeJumps[i] = t + rλ()
                if allowedToSwitchTo[nextArm(currentArm)] && state[currentArm] == 0
                    searchForNextArm()
                end
           end
       end
    end
    reward/tMax
end


text:
# This is based on Higham's own Matlab implementation based on his paper:
#   By N. J. Higham, 13/6/01, updated 30/1/13.
#   Reference:  N. J. Higham, Computing the nearest correlation
#   matrix---A problem from finance. IMA J. Numer. Anal.,
#   22(3):329-343, 2002.
#
# as presented on: http://nickhigham.wordpress.com/2013/02/13/the-nearest-correlation-matrix/
# It is typically not the fastest algorithm but should be ok for small (<500) dimensions...

# NEARCORR    Nearest correlation matrix.
#    (X, iter) = NEARCORR(A,TOL,FLAG,MAXITS,N_POS_EIG,W,PRNT)
#    finds the nearest correlation matrix to the symmetric matrix A.
#    TOL is a convergence tolerance, which defaults to 16*EPS.
#    If using FLAG == 1, TOL must be a 2-vector, with first component
#    the convergence tolerance and second component a tolerance
#    for defining "sufficiently positive" eigenvalues.
#    FLAG = 0: solve using full eigendecomposition (EIG).
#    FLAG = 1: treat as "highly non-positive definite A" and solve
#              using partial eigendecomposition (EIGS).
#    MAXITS is the maximum number of iterations (default 100, but may
#    need to be increased).
#    N_POS_EIG (optional) is the known number of positive eigenvalues of A.
#    W is a vector defining a diagonal weight matrix diag(W).
function nearcorr(A, tol = [16*Base.eps()], flag = 0, maxits = 100,
   n_pos_eig = false, w = ones(size(A, 1),1))

   if A != A'
      throw( ArgumentError("The matrix must be symmetric.") )
   end

   if length(tol) == 1 && flag == 1
      tol = size(A, 1) * Base.eps() * [1.0, 1.0]
   end

   n = size(A, 1)
   if flag == 1 && n_pos_eig == false
      d, V = eig(A)
      n_pos_eig = sum(d >= tol[2] * d[n])
   end

   X = copy(A)
   Y = copy(A)
   iter = 1
   rel_diffXY = rel_diffY = rel_diffX = Inf
   dS = zeros(size(A))
   w = w[:]
   Whalf = sqrt(w*w')

   while max(rel_diffX, rel_diffY, rel_diffXY) > tol[1]

      Xold = X
      R = X - dS
      R_wtd = Whalf .* R

      if flag == 0
         X = proj_spd(R_wtd)
      elseif flag == 1
         X, np = proj_spd_eigs(R_wtd, n_pos_eig, tol[2])
      end

      X =  X ./ Whalf
      dS = X - R
      Yold = Y
      Y = proj_unitdiag(X)
      normy = normfro(Y)
      rel_diffX = normfro(X-Xold)/normfro(X)
      rel_diffY = normfro(X-Xold)/normy
      rel_diffXY = normfro(Y-X)/normy

      iter = iter + 1;
      if iter > maxits
         throw( OverflowError("Stopped after $(maxits) iterations. Try increasing maxits.") )
      end

      X = Y

   end

   return X, iter
end

function proj_spd(A)
   d, V = eig(A)
   d[d .< 0.0] = 0.0 # Set negative eigen values to 0
   B = V * diagm(d) * V'
   (B+B')/2 # Ensure symmetric
end

function proj_spd_eigs(A, n_pos_eig, tol)
   k = n_pos_eig + 10 # 10 is safety factor.
   if k > size(A, 1)
      k = n_pos_eig
   end
   d, V = eigs(A, k; which = "LM")
   j = d .> tol*maximum(d)
   n_pos_eig_found = sum(j)
   B = V * diagm(d[j]) * V'
   B = (B+B')/2 # Ensure symmetric
   return B, n_pos_eig_found
end

function proj_unitdiag(A)
   B = copy(A)
   B[diagind(B)] = 1.0
   B
end


text:
export sage

"""
`sage(A)` takes a one- or two-dimensional array and prints it in a form
that can be pasted into SageMath.
"""
function sage(A::Array{T,2}) where {T}
    r, c = size(A)
    println("Matrix([")
    for i = 1:r
        print("[")
        for j = 1:c
            print(A[i, j])
            if j < c
                print(",")
            else
                print("]")
            end
        end
        if i < r
            println(",")
        else
            println("])")
        end
    end
end

function sage(A::Array{T,1}) where {T}
    n = length(A)
    print("Matrix([")
    for i = 1:n
        print("[$(A[i])]")
        if i < n
            print(",")
        else
            println("])")
        end
    end
end


text:
#This file is currently in development

using ISA
using Test
@testset "ISA.jl" begin

true

end


text:
using ClassifierMetrics
using DataFrames
using CSV
using Test
using LinearAlgebra
using Random
using Plots
Random.seed!(1234)

#utils
add_noise(label::Bool, λ=0.0) = label ? 1 - λ*rand() : λ*rand()

@testset "ClassifierMetrics" begin
  @testset "Binrary Metrics" begin
    include("binarymetrics.jl")
  end
  @testset "testing ROC" begin
    include("roc.jl")
  end
  @testset "testing PR" begin
    include("pr.jl")
  end
  @testset "testing DET" begin
    include("det.jl")
  end
  @testset "testing Plots" begin
    include("plots.jl")
  end
end


text:
# Add Max Trust Constraints

using SparseArrays

include("../constraints/constraints.jl")

"""
    makeMaxThrustConstraint(NSteps::Int64, nDim::Int64, tMax::Float64)

Creates an [`AL_Multiple_simpleCone`](@ref) constraint such that the trust
(or equivalently control) vector is capped by `tMax`.

Usage: `NSteps` states the number of timesteps in the trajectory. `nDim` states
the dimensionality of the problem (e.g. 1D, 2D, or 3D). `tMax` is the max
thrust.
"""
function makeMaxThrustConstraint(NSteps::Int64, nDim::Int64, tMax::Float64)
    # sizeTraj = 3 * NSteps * nDim + 2 * nDim

    # This helper function gets the index of ux for any k
    startUInd(k) = (2 * nDim + 1) + (3 * nDim * (k - 1))

    indicatorList = Int64[startUInd(k) for k in 1:NSteps]

    return makeAL_Multiple_simpleCone(tMax, nDim, indicatorList)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
module SimplePackage

  using Boot

  include_folder(SimplePackage, @__FILE__)

end


text:
"""
Collision avoidance system type. Partially overrides pilot intent if threat is detected.
"""
abstract type CAS end

"""
Simple rule-based collision avoidance system.
"""
Base.@kwdef mutable struct RuleBasedCAS <: CAS
	 hra::Symbol = :none  # horizontal resolution advisory
	 vra::Symbol = :none  # vertical resolution advisory
	 t_predict::R = 5.0   # prediction horizon [s]
	 d::R = 50.0          # lowest tolerable separation [ft]
	 psi_crit::R = 15.0   # determines reversal regions
end

reverse(ra::Symbol) = ra == :cw ? :ccw : ra == :ccw ? :cw : :none

"""
Determine relative heading region of aircraft for input to ruleset. Regions partition the
angle space defined in the `kernel` function, with horizontal velocity measured clockwise
from the +y direction. Intruder aircraft angle is flipped for symmetry.
"""
function region(psi::R, psi_crit::R; own::Bool=true)
	# reflect intruder angle
	xi = rad2deg(own ? psi : -psi)

	# discretize
	bins = [-180, -90 - psi_crit, -90, -90 + psi_crit, 0, 90]
	regs = [:IV, :iv, :i, :I, :II, :III]
	idx = searchsortedlast(bins, xi)
	return regs[idx]
end

function initialize(cas::RuleBasedCAS)
	cas.hra = :none
	cas.vra = :none
end

function update(cas::RuleBasedCAS, acs::Vector{Aircraft}, i::Z)
    # can only resolve pairwise interactions
    for j in filter(j -> j != i, 1:length(acs))
        advise(cas, acs[i], acs[j], i, j)
    end
end

"""
Use rule-based CAS to generate horizontal and vertical resolution advisories for a pair of
aircraft. There is a general turning rule, two reversal rules, and a vertical rule. The
aircraft containing the CAS is referred to as the ownship, while the other aircraft is the
intruder.

General turn rule: if ownship is aimed to the right of the intruder, turn clockwise;
otherwise, turn counterclockwise.

Reversal rules (override general turn rule):
1. If ownship is aimed within `psi_crit` of the intruder and intruder is aimed
more than `psi_crit` from the ownship in the same direction, turn in the opposite direction.
2. If both aircraft are aimed within `psi_crit` of each other in the same direction, use id
number to break the tie, where the aircraft with the higher id has lower precendence and
must switch direction.

Vertical rule: If ownship is above the intruder, ascend; if ownship is below the intruder,
descend. If aircraft are at the same altitude, lower-precendence aircraft must descend.
"""
function advise(cas::RuleBasedCAS, ac1::Aircraft, ac2::Aircraft, i::Z, j::Z)
	# assess threat within prediction window
	threat = assess_threat(cas, ac1, ac2)
	if !threat
		cas.hra = cas.vra = :none
		return
	end

	# get regions and establish precedence
	state = kernel(ac1, ac2; scale=false)
	r1 = region(state[2], cas.psi_crit, own=true)
	r2 = region(state[4], cas.psi_crit, own=false)
    defer = i > j

	# general turning rule
    cas.hra = r1 in (:i, :I, :II) ? :cw : :ccw

	# reversal rules
	rule1_ccw = r1 in (:i,) && r2 in (:I, :II)
    rule1_cw = r1 in (:iv,) && r2 in (:III, :IV)
	rule2 = r1 in (:i, :iv) && r1 == r2 && defer
    if rule1_cw || rule1_ccw || rule2
        cas.hra = reverse(cas.hra)
    end

	# vertical rule
	z1, z2 = ac1.r[3], ac2.r[3]
    if z1 > z2
        cas.vra = :asc
    elseif z1 < z2
        cas.vra = :desc
    else
        cas.vra = defer ? :desc : :none
    end
end

"""
Return true if projected aircraft trajectories pass within critical distance of each other
during prediction window. Analytical calculation naively assumes constant velocity.
"""
function assess_threat(cas::RuleBasedCAS, ac1::Aircraft, ac2::Aircraft)
	Δr = ac1.r - ac2.r
	Δv = ac1.v - ac2.v

	a = dot(Δv, Δv)
	b = 2 * dot(Δr, Δv)
	c = dot(Δr, Δr) - cas.d^2
	D = b^2 - 4 * a * c
	D < 0 && return false

	ts = (-b .+ [1, -1] * sqrt(D)) / (2 * a) # times of closest approach (may be in past)
	threat = any(0 .< ts .<= cas.t_predict)
	return threat
end

vec(cas::CAS) = [cas.hra; cas.vra]


text:
module Expect

import Base: convert, isassigned, getindex
export Expected, unexpected

struct Expected{T,E}
    v::Union{T,E}
    ok::Bool
end
convert(::Type{Expected{T,E}}, v::T) where {T, E} = Expected{T, E}(v, true)
convert(::Type{Expected{T,E}}, v::E) where {T, E} = Expected{T, E}(v, false)
convert(::Type{Expected{T,T}}, v::T) where {T}    = Expected{T, T}(v, true)

# `Unexpected` needed to disambiguate when `T == E` and want return an Unexpected
struct Unexpected{T}
    v::T
end
unexpected(x) = Unexpected(x)
convert(::Type{Expected{T,E}}, v::Unexpected{E}) where {T, E} = Expected{T, E}(v.v, false)

isassigned(e::Expected) = e.ok
@inline (getindex(e::Expected{T})::T) where {T} = isassigned(e) ? e.v : throw(e.v)

end # module


text:
include("../../dd/construction.jl")
include("../cut/cutgenerator.jl")

"""
Perform ECP using DD cuts.

# Input
- `m::JuMP.Model`: The current model to which the cuts will be added.
- `dd_list::Array{DecisionDiagram}`: An array of decision diagrams corresponding to inequalities of the current model.
- `constraint_list::Array{Function}`: An array of functions representing variable part of constraints with one-to-one correspondence to the DD array.
- `constant_list`: Array of constant part of constraints where all moved to the lhs of <= 0 form.
- `time_limit::Float64`: Time limit to run the algorithm (optional).
- `iter_max_ecp`: The maximum number of iterations to execute the algorithm (optional).
- `max_cut_num`: The maximum number of cuts (one cut for each violated inequality) to be added at EACH iteration (optional).
- `lp_to_ip_iter`: The number of iterations to solve as LP before switching to IP (optional).
- `lp_to_ip_tol`: The relative constraint violation threshold in the LP mode before swithcing to IP (optional).
- `obj_improvement_tol`: The improvement tolerance for the objective function (optional).
- `tol_num_ecp`: The number of past objective values (including the current one) wrt which the objective improvement tolerance is computed (optional).

# Output
- `x_output`: Optimal solution.
- `obj_output`: Optimal value.
- `status_output`: Status of the output model.
- `initial_lp_obj`: Objective value of the starting LP.
- `iter_ecp`: Total number of ECP iterations.
- `total_cuts`: Total number of added cuts.
- `ecp_time`: Total time elapsed.
- `relative_obj_improvement`: The relative improvement of the objective function value over the last few (tol_num_ecp) iterations.
"""

function ECP!(m::JuMP.Model, dd_list::Array{DecisionDiagram, 1}, constraint_list::Array{Function, 1}, constant_list::Array{Float64, 1};
    time_limit::Float64 = 10000.0, iter_max_ecp::Int = 1000, max_cut_num::Int = 10,
    lp_to_ip_iter::Int = 500, lp_to_ip_tol::Float64 = 0.00001, obj_improvement_tol::Float64 = 1e-8, tol_num_ecp::Int = 10)::
    Tuple{Array{Float64, 1}, Float64, Symbol, Float64, Int, Int, Float64, Float64}

    time_start::Float64 = time()

    @assert (length(dd_list) == length(constraint_list))
    @assert (length(constant_list) == length(constraint_list))

    n::Int = MathProgBase.numvar(m)
    # stopping criteria for the ECP algorithm
    obj_history_ecp::Array{Float64, 1} = [Inf for i=1:tol_num_ecp]    # the vector that stores the previous (improved) objective values for tolerance check
    relative_obj_improvement::Float64 = NaN    # The relative improvement of the objective function value over the last few (tol_num_ecp) iterations

    # function to define what stopping criterion must be used:
    # 0: time limit
    # 1: iteration number
    # 2: objective improvement: the relative difference between k consequtive improvements
    # 3: error on the optimality gap (TO DO: need a primal heuristic method)
    function stop_rule_ecp(stp::Array{Int64})::Bool   # multiple stopping criteria can be passed as an array
        p = true
        for v::Int in stp
            if v ==  0      # time limit criterion
                clock_end::Float64 = time()
                time_elapsed::Float64 = float(clock_end - clock_start)
                p = p && (time_elapsed <= time_limit)
                #println("v==0: ", p)
            end
            if v == 1       # the number of iteration criterion
                p = p && (iter_ecp <= iter_max_ecp)
                #println("v==1: ", p)
            end
            if v == 2       # uses the objective improvement criterion
                if cut_added   # if a cut has been added to the model, the `obj_history_ecp` vector must be updated and the criterion is checked
                    push!(obj_history_ecp, obj_val)
                    #popfirst!(obj_history_ecp)         #not supported for Julia .6
                    deleteat!(obj_history_ecp, 1)
                    relative_obj_improvement = (obj_history_ecp[1] - obj_val)/(abs(obj_val) + 1e-8)
                    p = p && (relative_obj_improvement > obj_improvement_tol)
                end
            end
        end
        return p    # is false when the stopping criterion is met
    end

    #*******************************************************************************

    # initializing cut gnerating paratmeters
    subgrad_flag::Bool = true       # if subgradient method is used to generate cuts
    CGLP_flag::Bool = false             # if CGLP method is used to generate cuts
    use_CGLP::Bool = false          # determines whether CGLP will be used at all in the ECP
    st_rule::Array{Int, 1} = [0, 1, 2]        # sets the stopping rule for the algorithm
    cut_added::Bool = false       # indicates whether a cut (of any type) has been added to the model
    con_num::Int = length(constraint_list)          # the number of constraint
    subgrad_start_point::Array{Float64, 2} = zeros(con_num,n)  # stores the starting point for each constraint to be used in the subgradient method
    violation_val::Array{Float64, 1} = Array{Float64}(  con_num)     # stores the violation value for each constraint at the current optimal point
    normalized_violation_val::Array{Float64, 1} = Array{Float64}(con_num)     # stores the violation value after being normalized wrt the rhs
    iter_ecp::Int = 1                                # iteration counter
    total_cuts::Int = 0                                 # total number of cuts added

    # set parameters to switch between solving LP and IP at iterations
    lp_flag::Bool = true          # if the LP method is used to solve the current model
    ip_flag::Bool = false         # if the IP method is used to solve the current model
    use_ip::Bool = false          # if the model is supposed to swotch to IP at some point (must be false for continuous programs)

    # initializing time parameters
    clock_start::Float64 = time()
    CplexSolver(CPX_PARAM_TILIM = time_limit)        #sets time limit parameter for the solver
    #ClpSolver(MaximumSeconds = time_limit)

    # get the initial solution of the model
    x = getindex(m, :x)                 # redefines variables in model m
    status::Symbol = solve(m)
    obj_val::Float64 = getobjectivevalue(m)
    initial_lp_obj = obj_val
    x_val = Array{Float64, 1}(n)    # DD variables and constraint functions are based on the flattened out form in MathProgBase, seo we need to store variable values in that format
    for i=1:n
        x_val[i] = getvalue(Variable(m,i))
        #x_val[i] = getvalue(x[i])
    end
    m_output::JuMP.Model = m
    x_output::Array{Float64, 1} = x_val         #array copy, change of x_val changes x_output
    obj_output::Float64 = obj_val
    status_output::Symbol = :Not_Optimal

    #************************************************************
    #performing the ecp method until the stopping criteria is met

    #println("@@@@@@@@@@@@@@@@@@@@", stop_rule_ecp(st_rule))

    while stop_rule_ecp(st_rule)

        # if the previous iteration has not added any cuts to the model:
        if cut_added == false && iter_ecp > 1
            # 1. if the algorithm is allowed to switch to CGLP
            if use_CGLP == true
                # 1.1. if CGLP method has not been activated: activate it
                if CGLP_flag == false
                    CGLP_flag = true
                    subgrad_flag = false
                # 1.2. if the CGLP has been solved to optimality: no cuts can be obtained from the current OA, B&B is needed to improve further
                elseif status_CGLP == :Optimal
                    println("OA (NOT global) optimal solution found!")
                    status_output = :OA_Optimal
                    break
                # 1.3. if the CGLP was not solved to optimality: cannot decide on the optimality of the current OA
                else
                    println("CGLP not optimal!")
                    status_output = :CGLP_Limit
                    break
                end
            # 2. the subgradient algorithm has not been able to detect a cut: cannot decide on the optimiality of the current OA
            else
                println("Subgradient method has not added cuts!")
                status_output = :Subgrad_Limit
                break
            end
        end

        #computing the violation value of the current point at each constraint
        for i::Int=1:con_num
            violation_val[i] = constraint_list[i](x_val) + constant_list[i]
            normalized_violation_val[i] = violation_val[i]/max(1e-9, abs(constant_list[i]))
        end

        sort_ind::Array{Int, 1} = sortperm(normalized_violation_val)     #stores the indices of the ascending sorted elements of violation_val

        # if the current point satisfies all constraints, it is global optimal
        if violation_val[sort_ind[end]] <= 0 + 1e-8
            println("Global optimal solution found!")
            status_output = :Global_Optimal
            break
        end

        # check whether to solve IP or LP based on the violation tolerance and the iteration number
        if (lp_flag == true) && (use_ip == true)
            if (normalized_violation_val[sort_ind[end]] <= lp_to_ip_tol) || (iter_ecp > lp_to_ip_iter)
                lp_flag = false
                ip_flag = true
                println("Switching to solve IP...")
                for i=1:n
                    setcategory(Variable(m,i), :Int)        # makes variable type with "linear index" i in its flattened out form, integer
                end
            end
        end

        cut_num = 0       # number of cuts added so far
        CGLP_cut = false    # CGLP cut has been added
        subgrad_cut = false     # subgradient cut has been added
        cut_added = false      # any cut has been added
        status_CGLP = :Optimal # it is assumed that CGLP will solve to optimality by default

        # analyze each constraint
        for i::Int=con_num:-1:1

            clock_end::Float64 = time()
            time_elapsed::Float64 = float(clock_end - clock_start)
            c_index::Int = sort_ind[i]

            # stop the loop if the inequality is satisfied or the number of allowable cuts is fulfilled or time exceeds the limit
            if (max_cut_num <= cut_num) || (violation_val[c_index] <= 0 + 1e-8) || (time_elapsed > time_limit)
                break
            end

            time_remain::Float64 = max(0.0, time_limit - time_elapsed)     # compute the remaining time

            if subgrad_flag == true       # uses subgradient method to generate cut
                obj1::Float64, cut_coef1::Array{Float64, 1}, cut_rhs1::Float64 = subgradient(dd_list[c_index], x_val, step_rule = 0, starting_point = subgrad_start_point[c_index,:], tilim = time_remain)

                if obj1 > 0 + 1e-7            # if the inequality cuts off the point
                    #println("\nsubgrad alg used\n")
                    cut_num += 1
                    subgrad_cut = true
                    #aff::JuMP.AffExpr = AffExpr(x, cut_coef1, 0)      # stores the affine expression of the inequality
                    aff = sum(Variable(m, i)*cut_coef1[i] for i=1:n)      # operates on variables in their original form (converted from flattened form)
                    @constraint(m, aff <= cut_rhs1)
                    subgrad_start_point[c_index,:] = cut_coef1      # uses the current cut's normal vector as the initial point for the next subgrad iteration
                end
            end

            if (CGLP_flag == true) && (use_CGLP == true)       # uses CGLP to generate cut

                println("SOLVING CGLP NOW")

                obj1, cut_coef1, cut_rhs1, status1 = CGLP(dd_list[c_index], x_val, tilim = time_remain)

                if (status1 == :Optimal) && (obj1 > 0 + 1e-8)
                    cut_num += 1
                    CGLP_cut = true
                    #aff = AffExpr(x, cut_coef1, 0)      #stores the affine expression of the inequality
                    aff = sum(Variable(m, i)*cut_coef1[i] for i=1:n)      # operates on variables in their original form (converted from flattened form)
                    @constraint(m, aff <= cut_rhs1)
                elseif status1 != :Optimal
                    status_CGLP == status1
                    info("CGLP terminated, unbounded or infeasible!")
                end

            end

        end

        cut_added = CGLP_cut || subgrad_cut

        #getting the initial solution of the model if a cut has been added
        if cut_added

            clock_end = time()
            time_elapsed = float(clock_end - clock_start)
            time_remain = max(0.0, time_limit - time_elapsed)     # compute the remaining time
            CplexSolver(CPX_PARAM_TILIM = time_remain)        #sets time limit parameter for the solver
            #ClpSolver(MaximumSeconds = time_remain)

            status = solve(m)
            println("\n\nLinear CUT:\t", MathProgBase.numlinconstr(m), "\n\n")
            obj_val = getobjectivevalue(m)
            println("Objective values: ", obj_val)
            if status == :Optimal               # the model has been solved to optimality
                for i=1:n
                    x_val[i] = getvalue(Variable(m,i))
                    #x_val[i] = getvalue(x[i])
                end
                x_output = x_val       # the change already takes effect because of array copy
                obj_output = obj_val
                status_output = :Not_Optimal
            end
        end

        total_cuts += cut_num
        iter_ecp += 1

    end

    time_end::Float64 = time()
    total_time_elapsed::Float64 = float(time_end - time_start)

    return x_output, obj_output, status_output, initial_lp_obj, iter_ecp, total_cuts, total_time_elapsed, relative_obj_improvement

end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
julia> using Sockets

julia> getaddrinfo("www.kame.net")
ip"203.178.141.194"

julia> getaddrinfo("www.kame.net", IPv6)
ip"2001:200:dff:fff1:216:3eff:feb1:44d7"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Load required package
using FisherySim
using Distributions
using StatsBase
using Test
using Random
using PDMats
using LinearAlgebra
using NeutralLandscapes

## Make reproducible
Random.seed!(1234)

## -----------------------------------------------------------------------------
## Construct a GriddedFisheryDomain
origin = (0.0, 0.0)
antipode = (100.0, 100.0)
n = (50, 50)
Ω = GriddedFisheryDomain(origin, antipode, n)
T = 10 # Number of years

include("test-fisherydomain.jl")

## -----------------------------------------------------------------------------
## Construct covariance kernels and covariance matrices over fishery domains
σ² = 5.0
ρ = 40.0

expkern = ExpCov(σ², ρ)
matkern = Matérn32Cov(σ², ρ)
ar1kern = AR1(0.1, 0.5)

expΣ = cov(expkern, Ω)
matΣ = cov(matkern, Ω)
ar1Σ = cov(ar1kern, T)

include("test-covkernels.jl")

## -----------------------------------------------------------------------------
## Define distributions over the domain
lognorm = DomainDistribution(LogNormal(-0.2^2 / 2, 0.2), Ω)
mvlognorm = DomainDistribution(MvLogNormal, ones(length(Ω)), matΣ, Ω)
mvnorm = DomainDistribution(MvNormal(zeros(length(Ω)), matΣ), Ω)
# matlognorm = DomainDistribution(MatrixLogNormal, 0.2ones(length(Ω), T), matΣ, ar1Σ, Ω)

ln_real = rand(lognorm)
mvln_real = rand(mvlognorm)
mvn_real = rand(mvnorm)
# matlognorm_real = rand(matlognorm)

dd1 = DomainDistribution(MidpointDisplacement(), Ω)
cdd1 = ClassifiedDomainDistribution(MidpointDisplacement(), Ω, 3)
cdd2 = ClassifiedDomainDistribution(MidpointDisplacement(), Ω, [0.1, 0.4])
ddist_vec = [mvnorm, dd1]
mdd1 = MultiDomainDistribution(ddist_vec)
bdd1 = BlendedDomainDistribution(mdd1, [0.3, 0.5])
bdd2 = BlendedDomainDistribution(ddist_vec, Ω, [0.5, 0.9])
bdd3 = BlendedDomainDistribution([mvlognorm, cdd1], Ω, [0.3, 2.0])

include("test-domaindistribution.jl")

## -----------------------------------------------------------------------------
## Habitats
hab1 = Habitat(rand(dd1))
hab2 = Habitat(rand.(ddist_vec))
include("test-habitat.jl")

## -----------------------------------------------------------------------------
## Generate bathymetry
μv = zeros(length(Ω))
# μm = zeros(size(Ω)...)
for idx in eachindex(Ω)
    μv[idx] = -1 / 2500 * first(Ω[idx]) * (first(Ω[idx]) - 100)
end

# hab_model = DomainDistribution(MvLogNormal, μv, matΣ, Ω)
hab_model = DomainDistribution(MvNormal(μv, matΣ), Ω)
hab = rand(hab_model)

hab2 = Habitat(rand(mdd1))

## -----------------------------------------------------------------------------
## Construct a movement model
h1pref(h) = exp(-h ^ 2 / 10)
h2pref(h) = 1
hab2pref = HabitatPreference(h1pref, h2pref)
dist = MovementRate(d -> exp(-d ^ 2 / 100))
move = MovementModel(hab2, hab2pref, dist, Ω)

# Find the equilibrium distribution using `eigvecs` to compare to to approximate
# method.
ews = eigvecs(move.M)
eqdist_eigvecs = real.(reshape(ews[:, end], size(Ω)...))
eqdist_eigvecs ./= sum(eqdist_eigvecs)

eqdist_ap0 = eqdist(move)
eqdist_ap = eqdist(move, 100.0)

include("test-movement.jl")

## -----------------------------------------------------------------------------
## Define Schaefer and PellaTomlinson dynamics models
r = 0.05
K = 100.0

schaef = Schaefer(r, K)
P1 = schaef(eqdist_ap)
Phalf = PopState(P1.P ./ 2)
Phalf1 = schaef(Phalf)

pt = PellaTomlinson(r, K, 3.39)
P1_pt = pt(eqdist_ap)
P1half_pt = pt(Phalf1)

unispt = StochasticProduction(PellaTomlinson(r, K, 3.39),
                              LogNormal(-0.2^2 / 2, 0.2))
multispt = StochasticProduction(
    PellaTomlinson(r, K, 3.39),
    mvlognorm)
P1_uspt = unispt(P1)
P1_mspt = multispt(P1)

include("test-pop_dynamics.jl")

## -----------------------------------------------------------------------------
## Define vessels and fleets
target_rand = RandomTargeting()
survey_stations = vec(LinearIndices(n)[2:2:48, 2:8:48])
target_fixed = FixedTargeting(survey_stations)
target_pref = PreferentialTargeting(10 .* eqdist_ap.P, Ω)

fixed_t = target(target_fixed, Ω)
reset!(target_fixed)
rand_t = target(target_rand, Ω)
reset!(target_rand)
pref_t = target(target_pref, Ω)
reset!(target_pref)

q_const = Catchability(0.2)
q_diff = Catchability(lognorm, q_const)
q_spat = Catchability(mvlognorm, q_const)
# q_sptemp = Catchability(matlognorm, q_const)
q_hab = HabitatCatchability(hab2, 0.1, h1 -> 1.0, h2 -> 1.2 * h2)

## These have to be fairly high to get many non-zero catches.
## These give ~10% zeros. These (especially ϕ?) may be good for testing,
## as they result in completely fishing out some cells, testing that check.
ξ = 1.9
ϕ = 1.9

v1 = Vessel(target_fixed, q_const, ξ, ϕ)
v2 = Vessel(target_rand, q_diff, ξ, ϕ)
v3 = Vessel(target_pref, q_spat, ξ, ϕ)
v4 = Vessel(target_pref, q_hab, ξ, ϕ)

P = PopState(ones(size(Ω)...) / length(Ω))
c1 = fish!(P, v1, Ω)
reset!(v1.target)

fleet = Fleet([v1, v2, v3, v4], [length(v1.target), 100, 1000, 1000])

c2 = fish!(P, fleet, Ω)
total_catch = sum(getfield.(c2, :catch_biomass))

include("test-vessels.jl")

## -----------------------------------------------------------------------------
## Fish a population using above definitions
P1 = eqdist_ap
Pvec, Cvec = simulate(P1, fleet, move, schaef, Ω, T)

Psums = sum.(Pvec)

function getcatch(CV::Vector{C}) where {C<:Catch}
    timevec = getfield.(CV, :time)
    catchvec = getfield.(CV, :catch_biomass)
    catchtot = zeros(length(unique(timevec)))
    for (t, c) in zip(timevec, catchvec)
        catchtot[t] += c
    end
    catchtot
end

Csums = getcatch(Cvec)

p2 = begin
         rempop = Psums[1] - Csums[1]
         rempop + rempop * r * (1 - rempop / K)
     end

include("test-simulate.jl")


text:
using Plots
using FFTW
using LinearAlgebra

include("fdtd.jl")

function A2b()
    # Spatial discretization
    n = 2000 + 1
    Δx = 1.0e-3

    # Time discretization
    Δt = Δx

    #steps to take
    steps = 2000

    #result at first and last step
    snapshots = Int.(round.(collect(range(1, stop = steps, length = 20))))

    #simulate
    result = FDTD1D(steps, snapshots, n, Δx, Δt)

    #save
    SaveFDTDResult("A2b.csv", result)

    #plot result
    ps = []
    for s = 1:length(snapshots)
        p = plot(result.xe, result.efield[s,:], ylims = [-1.5, 1.5])
        plot!(p, result.xh, result.hfield[s,:])
        push!(ps, p)
    end
    plot(ps..., layout = (length(snapshots), 1), show = true, size = [1200, 600])
end

function A2c()
    # Spatial discretization
    n = 2000 + 1
    Δx = 1.0e-3

    # Time discretization
    Δt = 1.01 * Δx
    steps = 130
    snapshots = [130, 140]
    result = FDTD1D(steps, snapshots, n, Δx, Δt)
    SaveFDTDResult("A2c.csv", result)
    p = plot(result.xe, result.efield[1,:], ylims = [-1.5, 1.5])
    plot!(p, result.xh, result.hfield[1,:])
end

function A2d()
    n = 2000 + 1
    Δx = 1.0e-3
    Δt = Δx
    steps = 2000
    snapshots = collect(1:steps)
    result = FDTD1D(steps, snapshots, n, Δx, Δt, 1.0)
    e = result.efield[:, 1000]
    h = result.hfield[:, 1000]
    ef = fft(e .* ((-1.0).^collect(1:length(e))))
    hf = fft(h .* ((-1.0).^collect(1:length(h))))
    s = 0.5 .* real.(ef .* hf)
    plot(abs.(s[950:1050]))
end

function A2f()
    n = 2000 + 1
    Δx = 1.0e-3
    Δt = Δx
    steps = 1500
    snapshots = Int.(round.(collect(range(1, stop = steps, length = 20))))
    result = FDTD1D(steps, snapshots, n, Δx, Δt, 1.0)
    SaveFDTDResult("A2f.csv", result)
end

function main()
    #A2b()
    #A2c()
    A2d()
    #A2f()
end

main()

max_stars_repo_name: JulienKluge/ComputationalPhotonics

max_issues_repo_name: JulienKluge/ComputationalPhotonics

max_forks_repo_name: JulienKluge/ComputationalPhotonics


text:
"""

The exported names are
$(EXPORTS)
"""
module ReferenceFEs

using Test
using DocStringExtensions
using LinearAlgebra
using Combinatorics
using FillArrays

using Gridap.Helpers
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials

using QuadGK: gauss
using FastGaussQuadrature: gaussjacobi
using FastGaussQuadrature: gausslegendre

import Gridap.Arrays: return_cache
import Gridap.Arrays: evaluate!
import Gridap.Arrays: return_type
import Gridap.Fields: evaluate
import Gridap.Fields: lazy_map
import Gridap.Polynomials: MonomialBasis

import Gridap.Polynomials: get_order
import Gridap.Polynomials: get_orders

import Gridap.Io: to_dict
import Gridap.Io: from_dict

import Base: ==

export Polytope
export ExtrusionPolytope
export get_extrusion
export get_faces
export get_dimranges
export get_dimrange
export get_vertex_coordinates
export get_facet_normal
export get_facet_orientations
export get_edge_tangent
export get_vertex_permutations
export get_face_dimranges
export get_face_coordinates
export is_simplex
export is_n_cube
export simplexify
export num_dims
export num_cell_dims
export num_point_dims
export num_vertices
export num_faces
export num_facets
export num_edges
export num_dims
export num_point_dims
export get_facedims
export get_offsets
export get_offset
export get_face_vertices
export get_reffaces
export get_face_type
export get_bounding_box
export get_face_vertex_permutations
export get_order
export get_orders
export test_polytope
export VERTEX
export SEGMENT
export TRI
export QUAD
export TET
export HEX
export WEDGE
export PYRAMID
export HEX_AXIS
export TET_AXIS
export INVALID_PERM

export ContraVariantPiolaMap

export Dof
export get_nodes
export get_face_moments
export get_face_nodes_dofs
export get_nodes
export evaluate!
export evaluate_dof
export return_cache
export return_type
export test_dof
export test_dof_array
# export evaluate_dof_array

export ReferenceFE
export ReferenceFEName
export GenericRefFE
export get_polytope
export get_prebasis
export get_dof_basis
export Conformity
export get_face_own_dofs
export get_face_own_dofs_permutations
export get_face_dofs
export get_own_dofs_permutations
export get_shapefuns
export compute_shapefuns
export test_reference_fe
export num_dofs
export expand_cell_data
export compress_cell_data

export LagrangianRefFE
export GenericLagrangianRefFE
export get_face_own_nodes_permutations
export get_own_nodes_permutations
export get_node_coordinates
export get_dof_to_node
export get_dof_to_comp
export get_node_and_comp_to_dof
export get_vertex_node
export num_nodes
export test_lagrangian_reference_fe

export LagrangianDofBasis
export compute_monomial_basis
export compute_own_nodes
export compute_face_orders
export compute_nodes
export compute_own_nodes_permutations
export compute_lagrangian_reffaces
export is_first_order
export is_Q
export is_P
export is_S

export MomentBasedDofBasis
export get_face_own_nodes
export get_face_nodes

export VERTEX1
export SEG2
export TRI3
export QUAD4
export TET4
export HEX8

export CONT
export DISC

export Conformity
export L2Conformity
export H1Conformity
export GradConformity
export CurlConformity
export DivConformity
export CDConformity

export SerendipityRefFE
export RaviartThomasRefFE
export NedelecRefFE
export BezierRefFE

export Lagrangian
export RaviartThomas
export Nedelec
export Bezier

export lagrangian
export raviart_thomas
export nedelec
export bezier

export Quadrature
export QuadratureName
export GenericQuadrature
export num_points
export get_coordinates
export get_weights
export get_name
export num_dims
export num_point_dims
export test_quadrature
export tensor_product
export duffy
export strang

include("Polytopes.jl")

include("ExtrusionPolytopes.jl")

include("Dofs.jl")

include("LagrangianDofBases.jl")

include("ReferenceFEInterfaces.jl")

include("LagrangianRefFEs.jl")

include("CLagrangianRefFEs.jl")

include("SerendipityRefFEs.jl")

include("PDiscRefFEs.jl")

include("CDLagrangianRefFEs.jl")

include("Quadratures.jl")

include("TensorProductQuadratures.jl")

include("DuffyQuadratures.jl")

include("StrangQuadratures.jl")

include("RaviartThomasRefFEs.jl")

include("NedelecRefFEs.jl")

include("MockDofs.jl")

include("BezierRefFEs.jl")

end # module


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Coefficients for the IGRF v12 to be used in the function `igrf12syn`.
#
#   This values were obtained from the FORTRAN source code [1].
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] https://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# The IGRF constant array in `igrf12.f` is composed by the concatenation of
# several arrays (`g0`, `g1`, ..., `gr`), each one with the coefficients related
# to a year. Hence, the final array to be used in the algorithm (`gh`) must be
# assembled using:
#
#   gh = [g0; g1; g2; ...; gr]
#
# The arrays have the following form:
#
#   -31543.,-2298., 5922., -677., 2905.,-1061.,  924., 1121., 1900
#
# The last element, which is the year, **must** be removed before those
# sub-arrays are concatenated.
#
# Furthermore, if the last element of the array has the following format:
#
#    N*0.0 ,
#
# it means that it has more `N` fields filed with 0.0. Then, in julia, it should
# be replaced by `zeros(N)...`.

const gh_igrf12 = [
      -31543.,-2298., 5922., -677., 2905.,-1061.,  924., 1121., # 1900
        1022.,-1469., -330., 1256.,    3.,  572.,  523.,  876., # 1900
         628.,  195.,  660.,  -69., -361., -210.,  134.,  -75., # 1900
        -184.,  328., -210.,  264.,   53.,    5.,  -33.,  -86., # 1900
        -124.,  -16.,    3.,   63.,   61.,   -9.,  -11.,   83., # 1900
        -217.,    2.,  -58.,  -35.,   59.,   36.,  -90.,  -69., # 1900
          70.,  -55.,  -45.,    0.,  -13.,   34.,  -10.,  -41., # 1900
          -1.,  -21.,   28.,   18.,  -12.,    6.,  -22.,   11., # 1900
           8.,    8.,   -4.,  -14.,   -9.,    7.,    1.,  -13., # 1900
           2.,    5.,   -9.,   16.,    5.,   -5.,    8.,  -18., # 1900
           8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., # 1900
          -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,   -1., # 1900
          -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1900
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1900
           0.,   -2.,    2.,    4.,    2.,    0.,    0.,   -6., # 1900
      -31464.,-2298., 5909., -728., 2928.,-1086., 1041., 1065., # 1905
        1037.,-1494., -357., 1239.,   34.,  635.,  480.,  880., # 1905
         643.,  203.,  653.,  -77., -380., -201.,  146.,  -65., # 1905
        -192.,  328., -193.,  259.,   56.,   -1.,  -32.,  -93., # 1905
        -125.,  -26.,   11.,   62.,   60.,   -7.,  -11.,   86., # 1905
        -221.,    4.,  -57.,  -32.,   57.,   32.,  -92.,  -67., # 1905
          70.,  -54.,  -46.,    0.,  -14.,   33.,  -11.,  -41., # 1905
           0.,  -20.,   28.,   18.,  -12.,    6.,  -22.,   11., # 1905
           8.,    8.,   -4.,  -15.,   -9.,    7.,    1.,  -13., # 1905
           2.,    5.,   -8.,   16.,    5.,   -5.,    8.,  -18., # 1905
           8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., # 1905
          -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,    0., # 1905
          -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1905
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1905
           0.,   -2.,    2.,    4.,    2.,    0.,    0.,   -6., # 1905
      -31354.,-2297., 5898., -769., 2948.,-1128., 1176., 1000., # 1910
        1058.,-1524., -389., 1223.,   62.,  705.,  425.,  884., # 1910
         660.,  211.,  644.,  -90., -400., -189.,  160.,  -55., # 1910
        -201.,  327., -172.,  253.,   57.,   -9.,  -33., -102., # 1910
        -126.,  -38.,   21.,   62.,   58.,   -5.,  -11.,   89., # 1910
        -224.,    5.,  -54.,  -29.,   54.,   28.,  -95.,  -65., # 1910
          71.,  -54.,  -47.,    1.,  -14.,   32.,  -12.,  -40., # 1910
           1.,  -19.,   28.,   18.,  -13.,    6.,  -22.,   11., # 1910
           8.,    8.,   -4.,  -15.,   -9.,    6.,    1.,  -13., # 1910
           2.,    5.,   -8.,   16.,    5.,   -5.,    8.,  -18., # 1910
           8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., # 1910
          -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,    0., # 1910
          -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1910
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1910
           0.,   -2.,    2.,    4.,    2.,    0.,    0.,   -6., # 1910
      -31212.,-2306., 5875., -802., 2956.,-1191., 1309.,  917., # 1915
        1084.,-1559., -421., 1212.,   84.,  778.,  360.,  887., # 1915
         678.,  218.,  631., -109., -416., -173.,  178.,  -51., # 1915
        -211.,  327., -148.,  245.,   58.,  -16.,  -34., -111., # 1915
        -126.,  -51.,   32.,   61.,   57.,   -2.,  -10.,   93., # 1915
        -228.,    8.,  -51.,  -26.,   49.,   23.,  -98.,  -62., # 1915
          72.,  -54.,  -48.,    2.,  -14.,   31.,  -12.,  -38., # 1915
           2.,  -18.,   28.,   19.,  -15.,    6.,  -22.,   11., # 1915
           8.,    8.,   -4.,  -15.,   -9.,    6.,    2.,  -13., # 1915
           3.,    5.,   -8.,   16.,    6.,   -5.,    8.,  -18., # 1915
           8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., # 1915
          -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,    0., # 1915
          -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1915
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1915
           0.,   -2.,    1.,    4.,    2.,    0.,    0.,   -6., # 1915
      -31060.,-2317., 5845., -839., 2959.,-1259., 1407.,  823., # 1920
        1111.,-1600., -445., 1205.,  103.,  839.,  293.,  889., # 1920
         695.,  220.,  616., -134., -424., -153.,  199.,  -57., # 1920
        -221.,  326., -122.,  236.,   58.,  -23.,  -38., -119., # 1920
        -125.,  -62.,   43.,   61.,   55.,    0.,  -10.,   96., # 1920
        -233.,   11.,  -46.,  -22.,   44.,   18., -101.,  -57., # 1920
          73.,  -54.,  -49.,    2.,  -14.,   29.,  -13.,  -37., # 1920
           4.,  -16.,   28.,   19.,  -16.,    6.,  -22.,   11., # 1920
           7.,    8.,   -3.,  -15.,   -9.,    6.,    2.,  -14., # 1920
           4.,    5.,   -7.,   17.,    6.,   -5.,    8.,  -19., # 1920
           8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., # 1920
          -3.,    1.,   -2.,   -2.,    9.,    2.,   10.,    0., # 1920
          -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1920
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1920
           0.,   -2.,    1.,    4.,    3.,    0.,    0.,   -6., # 1920
      -30926.,-2318., 5817., -893., 2969.,-1334., 1471.,  728., # 1925
        1140.,-1645., -462., 1202.,  119.,  881.,  229.,  891., # 1925
         711.,  216.,  601., -163., -426., -130.,  217.,  -70., # 1925
        -230.,  326.,  -96.,  226.,   58.,  -28.,  -44., -125., # 1925
        -122.,  -69.,   51.,   61.,   54.,    3.,   -9.,   99., # 1925
        -238.,   14.,  -40.,  -18.,   39.,   13., -103.,  -52., # 1925
          73.,  -54.,  -50.,    3.,  -14.,   27.,  -14.,  -35., # 1925
           5.,  -14.,   29.,   19.,  -17.,    6.,  -21.,   11., # 1925
           7.,    8.,   -3.,  -15.,   -9.,    6.,    2.,  -14., # 1925
           4.,    5.,   -7.,   17.,    7.,   -5.,    8.,  -19., # 1925
           8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., # 1925
          -3.,    1.,   -2.,   -2.,    9.,    2.,   10.,    0., # 1925
          -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1925
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1925
           0.,   -2.,    1.,    4.,    3.,    0.,    0.,   -6., # 1925
      -30805.,-2316., 5808., -951., 2980.,-1424., 1517.,  644., # 1930
        1172.,-1692., -480., 1205.,  133.,  907.,  166.,  896., # 1930
         727.,  205.,  584., -195., -422., -109.,  234.,  -90., # 1930
        -237.,  327.,  -72.,  218.,   60.,  -32.,  -53., -131., # 1930
        -118.,  -74.,   58.,   60.,   53.,    4.,   -9.,  102., # 1930
        -242.,   19.,  -32.,  -16.,   32.,    8., -104.,  -46., # 1930
          74.,  -54.,  -51.,    4.,  -15.,   25.,  -14.,  -34., # 1930
           6.,  -12.,   29.,   18.,  -18.,    6.,  -20.,   11., # 1930
           7.,    8.,   -3.,  -15.,   -9.,    5.,    2.,  -14., # 1930
           5.,    5.,   -6.,   18.,    8.,   -5.,    8.,  -19., # 1930
           8.,   10.,  -20.,    1.,   14.,  -12.,    5.,   12., # 1930
          -3.,    1.,   -2.,   -2.,    9.,    3.,   10.,    0., # 1930
          -2.,   -2.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1930
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1930
           0.,   -2.,    1.,    4.,    3.,    0.,    0.,   -6., # 1930
      -30715.,-2306., 5812.,-1018., 2984.,-1520., 1550.,  586., # 1935
        1206.,-1740., -494., 1215.,  146.,  918.,  101.,  903., # 1935
         744.,  188.,  565., -226., -415.,  -90.,  249., -114., # 1935
        -241.,  329.,  -51.,  211.,   64.,  -33.,  -64., -136., # 1935
        -115.,  -76.,   64.,   59.,   53.,    4.,   -8.,  104., # 1935
        -246.,   25.,  -25.,  -15.,   25.,    4., -106.,  -40., # 1935
          74.,  -53.,  -52.,    4.,  -17.,   23.,  -14.,  -33., # 1935
           7.,  -11.,   29.,   18.,  -19.,    6.,  -19.,   11., # 1935
           7.,    8.,   -3.,  -15.,   -9.,    5.,    1.,  -15., # 1935
           6.,    5.,   -6.,   18.,    8.,   -5.,    7.,  -19., # 1935
           8.,   10.,  -20.,    1.,   15.,  -12.,    5.,   11., # 1935
          -3.,    1.,   -3.,   -2.,    9.,    3.,   11.,    0., # 1935
          -2.,   -2.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1935
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1935
           0.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6., # 1935
      -30654.,-2292., 5821.,-1106., 2981.,-1614., 1566.,  528., # 1940
        1240.,-1790., -499., 1232.,  163.,  916.,   43.,  914., # 1940
         762.,  169.,  550., -252., -405.,  -72.,  265., -141., # 1940
        -241.,  334.,  -33.,  208.,   71.,  -33.,  -75., -141., # 1940
        -113.,  -76.,   69.,   57.,   54.,    4.,   -7.,  105., # 1940
        -249.,   33.,  -18.,  -15.,   18.,    0., -107.,  -33., # 1940
          74.,  -53.,  -52.,    4.,  -18.,   20.,  -14.,  -31., # 1940
           7.,   -9.,   29.,   17.,  -20.,    5.,  -19.,   11., # 1940
           7.,    8.,   -3.,  -14.,  -10.,    5.,    1.,  -15., # 1940
           6.,    5.,   -5.,   19.,    9.,   -5.,    7.,  -19., # 1940
           8.,   10.,  -21.,    1.,   15.,  -12.,    5.,   11., # 1940
          -3.,    1.,   -3.,   -2.,    9.,    3.,   11.,    1., # 1940
          -2.,   -2.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1940
          -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., # 1940
           0.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6., # 1940
      -30594.,-2285., 5810.,-1244., 2990.,-1702., 1578.,  477., # 1945
        1282.,-1834., -499., 1255.,  186.,  913.,  -11.,  944., # 1945
         776.,  144.,  544., -276., -421.,  -55.,  304., -178., # 1945
        -253.,  346.,  -12.,  194.,   95.,  -20.,  -67., -142., # 1945
        -119.,  -82.,   82.,   59.,   57.,    6.,    6.,  100., # 1945
        -246.,   16.,  -25.,   -9.,   21.,  -16., -104.,  -39., # 1945
          70.,  -40.,  -45.,    0.,  -18.,    0.,    2.,  -29., # 1945
           6.,  -10.,   28.,   15.,  -17.,   29.,  -22.,   13., # 1945
           7.,   12.,   -8.,  -21.,   -5.,  -12.,    9.,   -7., # 1945
           7.,    2.,  -10.,   18.,    7.,    3.,    2.,  -11., # 1945
           5.,  -21.,  -27.,    1.,   17.,  -11.,   29.,    3., # 1945
          -9.,   16.,    4.,   -3.,    9.,   -4.,    6.,   -3., # 1945
           1.,   -4.,    8.,   -3.,   11.,    5.,    1.,    1., # 1945
           2.,  -20.,   -5.,   -1.,   -1.,   -6.,    8.,    6., # 1945
          -1.,   -4.,   -3.,   -2.,    5.,    0.,   -2.,   -2., # 1945
      -30554.,-2250., 5815.,-1341., 2998.,-1810., 1576.,  381., # 1950
        1297.,-1889., -476., 1274.,  206.,  896.,  -46.,  954., # 1950
         792.,  136.,  528., -278., -408.,  -37.,  303., -210., # 1950
        -240.,  349.,    3.,  211.,  103.,  -20.,  -87., -147., # 1950
        -122.,  -76.,   80.,   54.,   57.,   -1.,    4.,   99., # 1950
        -247.,   33.,  -16.,  -12.,   12.,  -12., -105.,  -30., # 1950
          65.,  -55.,  -35.,    2.,  -17.,    1.,    0.,  -40., # 1950
          10.,   -7.,   36.,    5.,  -18.,   19.,  -16.,   22., # 1950
          15.,    5.,   -4.,  -22.,   -1.,    0.,   11.,  -21., # 1950
          15.,   -8.,  -13.,   17.,    5.,   -4.,   -1.,  -17., # 1950
           3.,   -7.,  -24.,   -1.,   19.,  -25.,   12.,   10., # 1950
           2.,    5.,    2.,   -5.,    8.,   -2.,    8.,    3., # 1950
         -11.,    8.,   -7.,   -8.,    4.,   13.,   -1.,   -2., # 1950
          13.,  -10.,   -4.,    2.,    4.,   -3.,   12.,    6., # 1950
           3.,   -3.,    2.,    6.,   10.,   11.,    3.,    8., # 1950
      -30500.,-2215., 5820.,-1440., 3003.,-1898., 1581.,  291., # 1955
        1302.,-1944., -462., 1288.,  216.,  882.,  -83.,  958., # 1955
         796.,  133.,  510., -274., -397.,  -23.,  290., -230., # 1955
        -229.,  360.,   15.,  230.,  110.,  -23.,  -98., -152., # 1955
        -121.,  -69.,   78.,   47.,   57.,   -9.,    3.,   96., # 1955
        -247.,   48.,   -8.,  -16.,    7.,  -12., -107.,  -24., # 1955
          65.,  -56.,  -50.,    2.,  -24.,   10.,   -4.,  -32., # 1955
           8.,  -11.,   28.,    9.,  -20.,   18.,  -18.,   11., # 1955
           9.,   10.,   -6.,  -15.,  -14.,    5.,    6.,  -23., # 1955
          10.,    3.,   -7.,   23.,    6.,   -4.,    9.,  -13., # 1955
           4.,    9.,  -11.,   -4.,   12.,   -5.,    7.,    2., # 1955
           6.,    4.,   -2.,    1.,   10.,    2.,    7.,    2., # 1955
          -6.,    5.,    5.,   -3.,   -5.,   -4.,   -1.,    0., # 1955
           2.,   -8.,   -3.,   -2.,    7.,   -4.,    4.,    1., # 1955
          -2.,   -3.,    6.,    7.,   -2.,   -1.,    0.,   -3., # 1955
      -30421.,-2169., 5791.,-1555., 3002.,-1967., 1590.,  206., # 1960
        1302.,-1992., -414., 1289.,  224.,  878., -130.,  957., # 1960
         800.,  135.,  504., -278., -394.,    3.,  269., -255., # 1960
        -222.,  362.,   16.,  242.,  125.,  -26., -117., -156., # 1960
        -114.,  -63.,   81.,   46.,   58.,  -10.,    1.,   99., # 1960
        -237.,   60.,   -1.,  -20.,   -2.,  -11., -113.,  -17., # 1960
          67.,  -56.,  -55.,    5.,  -28.,   15.,   -6.,  -32., # 1960
           7.,   -7.,   23.,   17.,  -18.,    8.,  -17.,   15., # 1960
           6.,   11.,   -4.,  -14.,  -11.,    7.,    2.,  -18., # 1960
          10.,    4.,   -5.,   23.,   10.,    1.,    8.,  -20., # 1960
           4.,    6.,  -18.,    0.,   12.,   -9.,    2.,    1., # 1960
           0.,    4.,   -3.,   -1.,    9.,   -2.,    8.,    3., # 1960
           0.,   -1.,    5.,    1.,   -3.,    4.,    4.,    1., # 1960
           0.,    0.,   -1.,    2.,    4.,   -5.,    6.,    1., # 1960
           1.,   -1.,   -1.,    6.,    2.,    0.,    0.,   -7., # 1960
      -30334.,-2119., 5776.,-1662., 2997.,-2016., 1594.,  114., # 1965
        1297.,-2038., -404., 1292.,  240.,  856., -165.,  957., # 1965
         804.,  148.,  479., -269., -390.,   13.,  252., -269., # 1965
        -219.,  358.,   19.,  254.,  128.,  -31., -126., -157., # 1965
         -97.,  -62.,   81.,   45.,   61.,  -11.,    8.,  100., # 1965
        -228.,   68.,    4.,  -32.,    1.,   -8., -111.,   -7., # 1965
          75.,  -57.,  -61.,    4.,  -27.,   13.,   -2.,  -26., # 1965
           6.,   -6.,   26.,   13.,  -23.,    1.,  -12.,   13., # 1965
           5.,    7.,   -4.,  -12.,  -14.,    9.,    0.,  -16., # 1965
           8.,    4.,   -1.,   24.,   11.,   -3.,    4.,  -17., # 1965
           8.,   10.,  -22.,    2.,   15.,  -13.,    7.,   10., # 1965
          -4.,   -1.,   -5.,   -1.,   10.,    5.,   10.,    1., # 1965
          -4.,   -2.,    1.,   -2.,   -3.,    2.,    2.,    1., # 1965
          -5.,    2.,   -2.,    6.,    4.,   -4.,    4.,    0., # 1965
           0.,   -2.,    2.,    3.,    2.,    0.,    0.,   -6., # 1965
      -30220.,-2068., 5737.,-1781., 3000.,-2047., 1611.,   25., # 1970
        1287.,-2091., -366., 1278.,  251.,  838., -196.,  952., # 1970
         800.,  167.,  461., -266., -395.,   26.,  234., -279., # 1970
        -216.,  359.,   26.,  262.,  139.,  -42., -139., -160., # 1970
         -91.,  -56.,   83.,   43.,   64.,  -12.,   15.,  100., # 1970
        -212.,   72.,    2.,  -37.,    3.,   -6., -112.,    1., # 1970
          72.,  -57.,  -70.,    1.,  -27.,   14.,   -4.,  -22., # 1970
           8.,   -2.,   23.,   13.,  -23.,   -2.,  -11.,   14., # 1970
           6.,    7.,   -2.,  -15.,  -13.,    6.,   -3.,  -17., # 1970
           5.,    6.,    0.,   21.,   11.,   -6.,    3.,  -16., # 1970
           8.,   10.,  -21.,    2.,   16.,  -12.,    6.,   10., # 1970
          -4.,   -1.,   -5.,    0.,   10.,    3.,   11.,    1., # 1970
          -2.,   -1.,    1.,   -3.,   -3.,    1.,    2.,    1., # 1970
          -5.,    3.,   -1.,    4.,    6.,   -4.,    4.,    0., # 1970
           1.,   -1.,    0.,    3.,    3.,    1.,   -1.,   -4., # 1970
      -30100.,-2013., 5675.,-1902., 3010.,-2067., 1632.,  -68., # 1975
        1276.,-2144., -333., 1260.,  262.,  830., -223.,  946., # 1975
         791.,  191.,  438., -265., -405.,   39.,  216., -288., # 1975
        -218.,  356.,   31.,  264.,  148.,  -59., -152., -159., # 1975
         -83.,  -49.,   88.,   45.,   66.,  -13.,   28.,   99., # 1975
        -198.,   75.,    1.,  -41.,    6.,   -4., -111.,   11., # 1975
          71.,  -56.,  -77.,    1.,  -26.,   16.,   -5.,  -14., # 1975
          10.,    0.,   22.,   12.,  -23.,   -5.,  -12.,   14., # 1975
           6.,    6.,   -1.,  -16.,  -12.,    4.,   -8.,  -19., # 1975
           4.,    6.,    0.,   18.,   10.,  -10.,    1.,  -17., # 1975
           7.,   10.,  -21.,    2.,   16.,  -12.,    7.,   10., # 1975
          -4.,   -1.,   -5.,   -1.,   10.,    4.,   11.,    1., # 1975
          -3.,   -2.,    1.,   -3.,   -3.,    1.,    2.,    1., # 1975
          -5.,    3.,   -2.,    4.,    5.,   -4.,    4.,   -1., # 1975
           1.,   -1.,    0.,    3.,    3.,    1.,   -1.,   -5., # 1975
      -29992.,-1956., 5604.,-1997., 3027.,-2129., 1663., -200., # 1980
        1281.,-2180., -336., 1251.,  271.,  833., -252.,  938., # 1980
         782.,  212.,  398., -257., -419.,   53.,  199., -297., # 1980
        -218.,  357.,   46.,  261.,  150.,  -74., -151., -162., # 1980
         -78.,  -48.,   92.,   48.,   66.,  -15.,   42.,   93., # 1980
        -192.,   71.,    4.,  -43.,   14.,   -2., -108.,   17., # 1980
          72.,  -59.,  -82.,    2.,  -27.,   21.,   -5.,  -12., # 1980
          16.,    1.,   18.,   11.,  -23.,   -2.,  -10.,   18., # 1980
           6.,    7.,    0.,  -18.,  -11.,    4.,   -7.,  -22., # 1980
           4.,    9.,    3.,   16.,    6.,  -13.,   -1.,  -15., # 1980
           5.,   10.,  -21.,    1.,   16.,  -12.,    9.,    9., # 1980
          -5.,   -3.,   -6.,   -1.,    9.,    7.,   10.,    2., # 1980
          -6.,   -5.,    2.,   -4.,   -4.,    1.,    2.,    0., # 1980
          -5.,    3.,   -2.,    6.,    5.,   -4.,    3.,    0., # 1980
           1.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6., # 1980
      -29873.,-1905., 5500.,-2072., 3044.,-2197., 1687., -306., # 1985
        1296.,-2208., -310., 1247.,  284.,  829., -297.,  936., # 1985
         780.,  232.,  361., -249., -424.,   69.,  170., -297., # 1985
        -214.,  355.,   47.,  253.,  150.,  -93., -154., -164., # 1985
         -75.,  -46.,   95.,   53.,   65.,  -16.,   51.,   88., # 1985
        -185.,   69.,    4.,  -48.,   16.,   -1., -102.,   21., # 1985
          74.,  -62.,  -83.,    3.,  -27.,   24.,   -2.,   -6., # 1985
          20.,    4.,   17.,   10.,  -23.,    0.,   -7.,   21., # 1985
           6.,    8.,    0.,  -19.,  -11.,    5.,   -9.,  -23., # 1985
           4.,   11.,    4.,   14.,    4.,  -15.,   -4.,  -11., # 1985
           5.,   10.,  -21.,    1.,   15.,  -12.,    9.,    9., # 1985
          -6.,   -3.,   -6.,   -1.,    9.,    7.,    9.,    1., # 1985
          -7.,   -5.,    2.,   -4.,   -4.,    1.,    3.,    0., # 1985
          -5.,    3.,   -2.,    6.,    5.,   -4.,    3.,    0., # 1985
           1.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6., # 1985
      -29775.,-1848., 5406.,-2131., 3059.,-2279., 1686., -373., # 1990
        1314.,-2239., -284., 1248.,  293.,  802., -352.,  939., # 1990
         780.,  247.,  325., -240., -423.,   84.,  141., -299., # 1990
        -214.,  353.,   46.,  245.,  154., -109., -153., -165., # 1990
         -69.,  -36.,   97.,   61.,   65.,  -16.,   59.,   82., # 1990
        -178.,   69.,    3.,  -52.,   18.,    1.,  -96.,   24., # 1990
          77.,  -64.,  -80.,    2.,  -26.,   26.,    0.,   -1., # 1990
          21.,    5.,   17.,    9.,  -23.,    0.,   -4.,   23., # 1990
           5.,   10.,   -1.,  -19.,  -10.,    6.,  -12.,  -22., # 1990
           3.,   12.,    4.,   12.,    2.,  -16.,   -6.,  -10., # 1990
           4.,    9.,  -20.,    1.,   15.,  -12.,   11.,    9., # 1990
          -7.,   -4.,   -7.,   -2.,    9.,    7.,    8.,    1., # 1990
          -7.,   -6.,    2.,   -3.,   -4.,    2.,    2.,    1., # 1990
          -5.,    3.,   -2.,    6.,    4.,   -4.,    3.,    0., # 1990
           1.,   -2.,    3.,    3.,    3.,   -1.,    0.,   -6., # 1990
      -29692.,-1784., 5306.,-2200., 3070.,-2366., 1681., -413., # 1995
        1335.,-2267., -262., 1249.,  302.,  759., -427.,  940., # 1995
         780.,  262.,  290., -236., -418.,   97.,  122., -306., # 1995
        -214.,  352.,   46.,  235.,  165., -118., -143., -166., # 1995
         -55.,  -17.,  107.,   68.,   67.,  -17.,   68.,   72., # 1995
        -170.,   67.,   -1.,  -58.,   19.,    1.,  -93.,   36., # 1995
          77.,  -72.,  -69.,    1.,  -25.,   28.,    4.,    5., # 1995
          24.,    4.,   17.,    8.,  -24.,   -2.,   -6.,   25., # 1995
           6.,   11.,   -6.,  -21.,   -9.,    8.,  -14.,  -23., # 1995
           9.,   15.,    6.,   11.,   -5.,  -16.,   -7.,   -4., # 1995
           4.,    9.,  -20.,    3.,   15.,  -10.,   12.,    8., # 1995
          -6.,   -8.,   -8.,   -1.,    8.,   10.,    5.,   -2., # 1995
          -8.,   -8.,    3.,   -3.,   -6.,    1.,    2.,    0., # 1995
          -4.,    4.,   -1.,    5.,    4.,   -5.,    2.,   -1., # 1995
           2.,   -2.,    5.,    1.,    1.,   -2.,    0.,   -7., # 1995
                                zeros(75)...,                   # 1995
      -29619.4,-1728.2, 5186.1,-2267.7, 3068.4,-2481.6, 1670.9, # 2000
        -458.0, 1339.6,-2288.0, -227.6, 1252.1,  293.4,  714.5, # 2000
        -491.1,  932.3,  786.8,  272.6,  250.0, -231.9, -403.0, # 2000
         119.8,  111.3, -303.8, -218.8,  351.4,   43.8,  222.3, # 2000
         171.9, -130.4, -133.1, -168.6,  -39.3,  -12.9,  106.3, # 2000
          72.3,   68.2,  -17.4,   74.2,   63.7, -160.9,   65.1, # 2000
          -5.9,  -61.2,   16.9,    0.7,  -90.4,   43.8,   79.0, # 2000
         -74.0,  -64.6,    0.0,  -24.2,   33.3,    6.2,    9.1, # 2000
          24.0,    6.9,   14.8,    7.3,  -25.4,   -1.2,   -5.8, # 2000
          24.4,    6.6,   11.9,   -9.2,  -21.5,   -7.9,    8.5, # 2000
         -16.6,  -21.5,    9.1,   15.5,    7.0,    8.9,   -7.9, # 2000
         -14.9,   -7.0,   -2.1,    5.0,    9.4,  -19.7,    3.0, # 2000
          13.4,   -8.4,   12.5,    6.3,   -6.2,   -8.9,   -8.4, # 2000
          -1.5,    8.4,    9.3,    3.8,   -4.3,   -8.2,   -8.2, # 2000
           4.8,   -2.6,   -6.0,    1.7,    1.7,    0.0,   -3.1, # 2000
           4.0,   -0.5,    4.9,    3.7,   -5.9,    1.0,   -1.2, # 2000
           2.0,   -2.9,    4.2,    0.2,    0.3,   -2.2,   -1.1, # 2000
          -7.4,    2.7,   -1.7,    0.1,   -1.9,    1.3,    1.5, # 2000
          -0.9,   -0.1,   -2.6,    0.1,    0.9,   -0.7,   -0.7, # 2000
           0.7,   -2.8,    1.7,   -0.9,    0.1,   -1.2,    1.2, # 2000
          -1.9,    4.0,   -0.9,   -2.2,   -0.3,   -0.4,    0.2, # 2000
           0.3,    0.9,    2.5,   -0.2,   -2.6,    0.9,    0.7, # 2000
          -0.5,    0.3,    0.3,    0.0,   -0.3,    0.0,   -0.4, # 2000
           0.3,   -0.1,   -0.9,   -0.2,   -0.4,   -0.4,    0.8, # 2000
          -0.2,   -0.9,   -0.9,    0.3,    0.2,    0.1,    1.8, # 2000
          -0.4,   -0.4,    1.3,   -1.0,   -0.4,   -0.1,    0.7, # 2000
           0.7,   -0.4,    0.3,    0.3,    0.6,   -0.1,    0.3, # 2000
           0.4,   -0.2,    0.0,   -0.5,    0.1,   -0.9,         # 2000
     -29554.63,-1669.05, 5077.99,-2337.24, 3047.69,-2594.50,    # 2005
       1657.76, -515.43, 1336.30,-2305.83, -198.86, 1246.39,    # 2005
        269.72,  672.51, -524.72,  920.55,  797.96,  282.07,    # 2005
        210.65, -225.23, -379.86,  145.15,  100.00, -305.36,    # 2005
       -227.00,  354.41,   42.72,  208.95,  180.25, -136.54,    # 2005
       -123.45, -168.05,  -19.57,  -13.55,  103.85,   73.60,    # 2005
         69.56,  -20.33,   76.74,   54.75, -151.34,   63.63,    # 2005
        -14.58,  -63.53,   14.58,    0.24,  -86.36,   50.94,    # 2005
         79.88,  -74.46,  -61.14,   -1.65,  -22.57,   38.73,    # 2005
          6.82,   12.30,   25.35,    9.37,   10.93,    5.42,    # 2005
        -26.32,    1.94,   -4.64,   24.80,    7.62,   11.20,    # 2005
        -11.73,  -20.88,   -6.88,    9.83,  -18.11,  -19.71,    # 2005
         10.17,   16.22,    9.36,    7.61,  -11.25,  -12.76,    # 2005
         -4.87,   -0.06,    5.58,    9.76,  -20.11,    3.58,    # 2005
         12.69,   -6.94,   12.67,    5.01,   -6.72,  -10.76,    # 2005
         -8.16,   -1.25,    8.10,    8.76,    2.92,   -6.66,    # 2005
         -7.73,   -9.22,    6.01,   -2.17,   -6.12,    2.19,    # 2005
          1.42,    0.10,   -2.35,    4.46,   -0.15,    4.76,    # 2005
          3.06,   -6.58,    0.29,   -1.01,    2.06,   -3.47,    # 2005
          3.77,   -0.86,   -0.21,   -2.31,   -2.09,   -7.93,    # 2005
          2.95,   -1.60,    0.26,   -1.88,    1.44,    1.44,    # 2005
         -0.77,   -0.31,   -2.27,    0.29,    0.90,   -0.79,    # 2005
         -0.58,    0.53,   -2.69,    1.80,   -1.08,    0.16,    # 2005
         -1.58,    0.96,   -1.90,    3.99,   -1.39,   -2.15,    # 2005
         -0.29,   -0.55,    0.21,    0.23,    0.89,    2.38,    # 2005
         -0.38,   -2.63,    0.96,    0.61,   -0.30,    0.40,    # 2005
          0.46,    0.01,   -0.35,    0.02,   -0.36,    0.28,    # 2005
          0.08,   -0.87,   -0.49,   -0.34,   -0.08,    0.88,    # 2005
         -0.16,   -0.88,   -0.76,    0.30,    0.33,    0.28,    # 2005
          1.72,   -0.43,   -0.54,    1.18,   -1.07,   -0.37,    # 2005
         -0.04,    0.75,    0.63,   -0.26,    0.21,    0.35,    # 2005
          0.53,   -0.05,    0.38,    0.41,   -0.22,   -0.10,    # 2005
         -0.57,   -0.18,   -0.82,                               # 2005
     -29496.57,-1586.42, 4944.26,-2396.06, 3026.34,-2708.54,    # 2010
       1668.17, -575.73, 1339.85,-2326.54, -160.40, 1232.10,    # 2010
        251.75,  633.73, -537.03,  912.66,  808.97,  286.48,    # 2010
        166.58, -211.03, -356.83,  164.46,   89.40, -309.72,    # 2010
       -230.87,  357.29,   44.58,  200.26,  189.01, -141.05,    # 2010
       -118.06, -163.17,   -0.01,   -8.03,  101.04,   72.78,    # 2010
         68.69,  -20.90,   75.92,   44.18, -141.40,   61.54,    # 2010
        -22.83,  -66.26,   13.10,    3.02,  -78.09,   55.40,    # 2010
         80.44,  -75.00,  -57.80,   -4.55,  -21.20,   45.24,    # 2010
          6.54,   14.00,   24.96,   10.46,    7.03,    1.64,    # 2010
        -27.61,    4.92,   -3.28,   24.41,    8.21,   10.84,    # 2010
        -14.50,  -20.03,   -5.59,   11.83,  -19.34,  -17.41,    # 2010
         11.61,   16.71,   10.85,    6.96,  -14.05,  -10.74,    # 2010
         -3.54,    1.64,    5.50,    9.45,  -20.54,    3.45,    # 2010
         11.51,   -5.27,   12.75,    3.13,   -7.14,  -12.38,    # 2010
         -7.42,   -0.76,    7.97,    8.43,    2.14,   -8.42,    # 2010
         -6.08,  -10.08,    7.01,   -1.94,   -6.24,    2.73,    # 2010
          0.89,   -0.10,   -1.07,    4.71,   -0.16,    4.44,    # 2010
          2.45,   -7.22,   -0.33,   -0.96,    2.13,   -3.95,    # 2010
          3.09,   -1.99,   -1.03,   -1.97,   -2.80,   -8.31,    # 2010
          3.05,   -1.48,    0.13,   -2.03,    1.67,    1.65,    # 2010
         -0.66,   -0.51,   -1.76,    0.54,    0.85,   -0.79,    # 2010
         -0.39,    0.37,   -2.51,    1.79,   -1.27,    0.12,    # 2010
         -2.11,    0.75,   -1.94,    3.75,   -1.86,   -2.12,    # 2010
         -0.21,   -0.87,    0.30,    0.27,    1.04,    2.13,    # 2010
         -0.63,   -2.49,    0.95,    0.49,   -0.11,    0.59,    # 2010
          0.52,    0.00,   -0.39,    0.13,   -0.37,    0.27,    # 2010
          0.21,   -0.86,   -0.77,   -0.23,    0.04,    0.87,    # 2010
         -0.09,   -0.89,   -0.87,    0.31,    0.30,    0.42,    # 2010
          1.66,   -0.45,   -0.59,    1.08,   -1.14,   -0.31,    # 2010
         -0.07,    0.78,    0.54,   -0.18,    0.10,    0.38,    # 2010
          0.49,    0.02,    0.44,    0.42,   -0.25,   -0.26,    # 2010
         -0.53,   -0.26,   -0.79,                               # 2010
     -29442.0,-1501.0, 4797.1,-2445.1, 3012.9,-2845.6, 1676.7,  # 2015
       -641.9, 1350.7,-2352.3, -115.3, 1225.6,  244.9,  582.0,  # 2015
       -538.4,  907.6,  813.7,  283.3,  120.4, -188.7, -334.9,  # 2015
        180.9,   70.4, -329.5, -232.6,  360.1,   47.3,  192.4,  # 2015
        197.0, -140.9, -119.3, -157.5,   16.0,    4.1,  100.2,  # 2015
         70.0,   67.7,  -20.8,   72.7,   33.2, -129.9,   58.9,  # 2015
        -28.9,  -66.7,   13.2,    7.3,  -70.9,   62.6,   81.6,  # 2015
        -76.1,  -54.1,   -6.8,  -19.5,   51.8,    5.7,   15.0,  # 2015
         24.4,    9.4,    3.4,   -2.8,  -27.4,    6.8,   -2.2,  # 2015
         24.2,    8.8,   10.1,  -16.9,  -18.3,   -3.2,   13.3,  # 2015
        -20.6,  -14.6,   13.4,   16.2,   11.7,    5.7,  -15.9,  # 2015
         -9.1,   -2.0,    2.1,    5.4,    8.8,  -21.6,    3.1,  # 2015
         10.8,   -3.3,   11.8,    0.7,   -6.8,  -13.3,   -6.9,  # 2015
         -0.1,    7.8,    8.7,    1.0,   -9.1,   -4.0,  -10.5,  # 2015
          8.4,   -1.9,   -6.3,    3.2,    0.1,   -0.4,    0.5,  # 2015
          4.6,   -0.5,    4.4,    1.8,   -7.9,   -0.7,   -0.6,  # 2015
          2.1,   -4.2,    2.4,   -2.8,   -1.8,   -1.2,   -3.6,  # 2015
         -8.7,    3.1,   -1.5,   -0.1,   -2.3,    2.0,    2.0,  # 2015
         -0.7,   -0.8,   -1.1,    0.6,    0.8,   -0.7,   -0.2,  # 2015
          0.2,   -2.2,    1.7,   -1.4,   -0.2,   -2.5,    0.4,  # 2015
         -2.0,    3.5,   -2.4,   -1.9,   -0.2,   -1.1,    0.4,  # 2015
          0.4,    1.2,    1.9,   -0.8,   -2.2,    0.9,    0.3,  # 2015
          0.1,    0.7,    0.5,   -0.1,   -0.3,    0.3,   -0.4,  # 2015
          0.2,    0.2,   -0.9,   -0.9,   -0.1,    0.0,    0.7,  # 2015
          0.0,   -0.9,   -0.9,    0.4,    0.4,    0.5,    1.6,  # 2015
         -0.5,   -0.5,    1.0,   -1.2,   -0.2,   -0.1,    0.8,  # 2015
          0.4,   -0.1,   -0.1,    0.3,    0.4,    0.1,    0.5,  # 2015
          0.5,   -0.3,   -0.4,   -0.4,   -0.3,   -0.8,          # 2015
         10.3,   18.1,  -26.6,   -8.7,   -3.3,  -27.4,    2.1,  # 2017
        -14.1,    3.4,   -5.5,    8.2,   -0.7,   -0.4,  -10.1,  # 2017
          1.8,   -0.7,    0.2,   -1.3,   -9.1,    5.3,    4.1,  # 2017
          2.9,   -4.3,   -5.2,   -0.2,    0.5,    0.6,   -1.3,  # 2017
          1.7,   -0.1,   -1.2,    1.4,    3.4,    3.9,    0.0,  # 2017
         -0.3,   -0.1,    0.0,   -0.7,   -2.1,    2.1,   -0.7,  # 2017
         -1.2,    0.2,    0.3,    0.9,    1.6,    1.0,    0.3,  # 2017
         -0.2,    0.8,   -0.5,    0.4,    1.3,   -0.2,    0.1,  # 2017
         -0.3,   -0.6,   -0.6,   -0.8,    0.1,    0.2,   -0.2,  # 2017
          0.2,    0.0,   -0.3,   -0.6,    0.3,    0.5,    0.1,  # 2017
         -0.2,    0.5,    0.4,   -0.2,    0.1,   -0.3,   -0.4,  # 2017
          0.3,    0.3,    0.0, zeros(115)...                    # 2017
         ]


text:
pcountprimes(to; segment_length = 1024 * 32, threads = Threads.nthreads()) = pcountprimes(1, to, segment_length = segment_length, threads = threads)

"""
Get segment `i` when dividing `from:to` into `n` pieces
"""
function segment(from, to, i, n)
    div, rem = divrem(length(from:to), n)
    start = from + (div + 1) * min(i - 1, rem) + div * max(0, i - 1 - rem)
    stop = from + (div + 1) * min(i, rem) + div * max(0, i - rem) - 1
    start, stop
end

function pcountprimes(from, to; segment_length = 1024 * 32, threads = Threads.nthreads())

    from = max(2, from)

    counts = zeros(Int, threads)

    Threads.@sync for i in Base.OneTo(threads)
        Threads.@spawn begin
            start, stop = segment(from, to, i, threads)
            counts[i] = countprimes(start, stop, segment_length)
        end
    end

    return sum(counts)
end


text:
## ---------------------------------------------------
function _key_by(key)
    _toml_order = ["name", "uuid", "authors", "version", "deps", "compat", "extras", "targets"]
    idx = findfirst(isequal(key), _toml_order)
    isnothing(idx) ? 9223372036854775807 : idx
end

function _save_project(projfile, projdict)
    open(projfile, "w") do io
        TOML.print(io, projdict; sorted=true, by=_key_by)
    end
end

_load_project(projfile) = TOML.parsefile(projfile)


text:
module SerializationTests

using BenchmarkTools
using Test

eq(x::T, y::T) where {T<:Union{values(BenchmarkTools.SUPPORTED_TYPES)...}} =
    all(i->eq(getfield(x, i), getfield(y, i)), 1:fieldcount(T))
eq(x::T, y::T) where {T} = isapprox(x, y)

function withtempdir(f::Function)
    d = mktempdir()
    try
        cd(f, d)
    finally
        rm(d, force=true, recursive=true)
    end
    nothing
end

@testset "Successful (de)serialization" begin
    b = @benchmarkable sin(1)
    tune!(b)
    bb = run(b)

    withtempdir() do
        tmp = joinpath(pwd(), "tmp.json")

        BenchmarkTools.save(tmp, b.params, bb)
        @test isfile(tmp)

        results = BenchmarkTools.load(tmp)
        @test results isa Vector{Any}
        @test length(results) == 2
        @test eq(results[1], b.params)
        @test eq(results[2], bb)
    end

    # Nested BenchmarkGroups
    withtempdir() do
        tmp = joinpath(pwd(), "tmp.json")

        g = BenchmarkGroup()
        g["a"] = BenchmarkGroup()
        g["b"] = BenchmarkGroup()
        g["c"] = BenchmarkGroup()
        BenchmarkTools.save(tmp, g)

        results = BenchmarkTools.load(tmp)[1]
        @test results isa BenchmarkGroup
        @test all(v->v isa BenchmarkGroup, values(results.data))
    end
end

@testset "Deprecated behaviors" begin
    b = @benchmarkable sin(1)
    tune!(b)
    bb = run(b)

    @test_throws ArgumentError BenchmarkTools.save("x.jld", b.params)
    @test_throws ArgumentError BenchmarkTools.save("x.txt", b.params)
    @test_throws ArgumentError BenchmarkTools.save("x.json")
    @test_throws ArgumentError BenchmarkTools.save("x.json", 1)

    withtempdir() do
        tmp = joinpath(pwd(), "tmp.json")
        @test_logs (:warn, r"Naming variables") BenchmarkTools.save(tmp, "b", b.params)
        @test isfile(tmp)
        results = BenchmarkTools.load(tmp)
        @test length(results) == 1
        @test eq(results[1], b.params)
    end

    @test_throws ArgumentError BenchmarkTools.load("x.jld")
    @test_throws ArgumentError BenchmarkTools.load("x.txt")
    @test_throws ArgumentError BenchmarkTools.load("x.json", "b")
end

@testset "Error checking" begin
    withtempdir() do
        tmp = joinpath(pwd(), "tmp.json")
        open(tmp, "w") do f
            print(f, """
            {"never":1,"gonna":[{"give":3,"you":4,"up":5}]}
            """)
        end
        try
            BenchmarkTools.load(tmp)
            error("madness")
        catch err
            # This function thows a bunch of errors, so test for this specifically
            @test occursin("Unexpected JSON format", err.msg)
        end
    end

    @test_throws ArgumentError BenchmarkTools.recover([1])
end

end # module


text:
function parse(io::IO)
    d = Decoder()
    decode(d, JSON.parse(io))
end

function parse(s::AbstractString)
    d = Decoder()
    decode(d, JSON.parse(s))
end


text: include("testPlotModel.jl")

max_stars_repo_name: JuliaInv/jInvVisPyPlot

max_issues_repo_name: JuliaInv/jInvVisPyPlot

max_forks_repo_name: JuliaInv/jInvVisPyPlot


text:
#!/usr/bin/env julia6
using ReadGlobal, Decomp

function makeModel(Model::String)

    N, Fil = getnfilter()

    outp = "Model_$Model"
    isdir(outp) || mkdir(outp)

    tm11 = "$(outp)/Tm11"
    tm12 = "$(outp)/Tm12"
    tm13 = "$(outp)/Tm13"
    tm22 = "$(outp)/Tm22"
    tm23 = "$(outp)/Tm23"
    tm33 = "$(outp)/Tm33"

    t11 = "$(Fil)T11_N$N"
    t12 = "$(Fil)T12_N$N"
    t13 = "$(Fil)T13_N$N"
    t22 = "$(Fil)T22_N$N"
    t23 = "$(Fil)T23_N$N"
    t33 = "$(Fil)T33_N$N"

    d11 = "$(Fil)D11_N$N"
    d12 = "$(Fil)D12_N$N"
    d13 = "$(Fil)D13_N$N"
    d22 = "$(Fil)D22_N$N"
    d23 = "$(Fil)D23_N$N"
    d33 = "$(Fil)D33_N$N"

    p11 = "$(Fil)P11_N$N"
    p12 = "$(Fil)P12_N$N"
    p13 = "$(Fil)P13_N$N"
    p22 = "$(Fil)P22_N$N"
    p23 = "$(Fil)P23_N$N"
    p33 = "$(Fil)P33_N$N"

    if Model in ("I", "II", "III")
        if Model == "I"
            doinchunks(propmodel,input=(t11,t22,t33,t12,t13,t23,d11,d22,d33,d12,d13,d23),
                          output=(tm11,tm22,tm33,tm12,tm13,tm23,"$(outp)/alpha","$(outp)/ratio"))
        elseif Model == "II"
            doinchunks(inphmodel,input=(t11,t22,t33,t12,t13,t23,d11,d22,d33,d12,d13,d23),
                          output=(tm11,tm22,tm33,tm12,tm13,tm23,"$(outp)/alpha0","$(outp)/alpha1","$(outp)/alpha2","$(outp)/ratio"))
        elseif Model == "III"
            doinchunks(propmodel,input=(t11,t22,t33,t12,t13,t23,p11,p22,p33,p12,p13,p23),
                          output=(tm11,tm22,tm33,tm12,tm13,tm23,"$(outp)/alpha","$(outp)/ratio"))
        end

    elseif Model in ("V", "VI")
        if Model == "V"
            doinchunks(mixmodel,input=(t11,t22,t33,t12,t13,t23,
                          "Model_I/Tm11","Model_I/Tm22","Model_I/Tm33","Model_I/Tm12","Model_I/Tm13","Model_I/Tm23",
                          "Model_III/Tm11","Model_III/Tm22","Model_III/Tm33","Model_III/Tm12","Model_III/Tm13","Model_III/Tm23"),
                           output=(tm11,tm22,tm33,tm12,tm13,tm23,"$(outp)/ratio"))
        else
            doinchunks(mixmodel,input=(t11,t22,t33,t12,t13,t23,
                          "Model_II/Tm11","Model_II/Tm22","Model_II/Tm33","Model_II/Tm12","Model_II/Tm13","Model_II/Tm23",
                          "Model_III/Tm11","Model_III/Tm22","Model_III/Tm33","Model_III/Tm12","Model_III/Tm13","Model_III/Tm23"),
                           output=(tm11,tm22,tm33,tm12,tm13,tm23,"$(outp)/ratio"))
        end
    end

    return 0
end

for f::String in ARGS
    makeModel(f)
end