text:
# ===========================================================================
# Conservative linearization and range enclosures
# ===========================================================================

# ------------------------------------------
# Enclosing the range of 1/2 * v' * Hi * v
# ------------------------------------------

abstract type AbstractRangeEnclosureMethod end

struct IntervalArithmeticEnclosure <: AbstractRangeEnclosureMethod end

# range enclosure of the Hessian using IA
function _enclose_hessian(Hsym::Vector{<:AbstractMatrix}, var, X0::IntervalBox, c, alg::AbstractRangeEnclosureMethod=IntervalArithmeticEnclosure())
    d = var - c
    ex = [1/2 * d' * Hi * d for Hi in Hsym]
    dict = [vi => ci for (vi, ci) in zip(var, X0.v)]
    y = Symbolics.substitute.(ex, Ref(dict))
    return IntervalBox(Symbolics.value.(y))
end

# fallback for generic LazySet
function _enclose_hessian(Hsym::Vector{<:AbstractMatrix}, var, X0::LazySet, c, alg::AbstractRangeEnclosureMethod=IntervalArithmeticEnclosure())
    B0 = box_approximation(X0)
    return convert(IntervalBox, B0)
end

# ------------------------------------------
# Conservative linearization
# ------------------------------------------

_new_variables(n) = [Num(SymbolicUtils.Sym{Real}(:x, i)) for i in 1:n]

function _linearize(ivp::InitialValueProblem,
                    alg::AbstractRangeEnclosureMethod=IntervalArithmeticEnclosure())

    n = statedim(ivp)
    var = _new_variables(n)  # @variables x[1:n]

    J = _jacobian_function(ivp, var)
    H = _hessian_function(ivp, var)
    Hsym = _hessian_expression(ivp, var)

    # linearization center
    X0 = initial_states(ivp)
    c = center(X0)

    # evaluations
    f = ivp.s.f
    ff = eval(build_function(f, var)[1]) # TODO generalize with VectorField?
    fc = ff(c)
    Jc = J(c)
    A = Jc
    b = fc - Jc*c

    # input set
    B = Matrix(1.0I, n, n)
    UH = _enclose_hessian(Hsym, var, X0, c, alg)
    UH = convert(Hyperrectangle, UH)
    ivplin = @ivp(x' = A*x + b + B*u, x(0) ∈ X0, x ∈ Universe(n), u ∈ UH)

    return ivplin
end


text:
"""
Implementation of Koza-type (tree-based) Genetic Programming

The constructor takes following keyword arguments:

- `populationSize`: The size of the population
- `initialization`: A strategy for population initialization (default: `:half`)
  - Possible values: `:grow`, `:full`, `:half`
- `evaluate`: A collection of fitness cases required for evaluation of a fitness function

"""
@kwdef struct TreeGP{T} <: AbstractOptimizer
    populationSize::Integer = 50
    terminals::Vector{Symbol} = [:x, :rand]
    functions::Vector{Symbol} = [:+, :-, :*, :/, :^]
    mindepth::Int = 0
    maxdepth::Int = 10
    initialization::Symbol = :grow
    evaluate::AbstractVector{T} = zeros(T, 0)
end

"""
    initial_population(method::TreeGP, expr::{Expr,Nothing})

Initialize a random population of expressions derived from `expr`.
"""
function initial_population(method::GP, expr::{Expr,Nothing}))
    n = population_size(method)

    return [genrandexpr(method.terminals, method.functions,
                        mindepth=method.mindepth, maxdepth=method.maxdepth)
            for i in 1:n]
end

arity(sym) = 2
function genrandexpr(terms, funcs; maxdepth=2, mindepth=1, method=:grow)
    # println("depth: $maxdepth => $mindepth")
    tl = length(terms)
    fl = length(funcs)
    expr = if (maxdepth <= 0 || ( method == :grow && rand() < tl/(tl+fl) )) && mindepth <= 0
        termsym = rand(terms)
        term = termsym == :rand ? rand() : termsym
        # println("term: $term")
        term
    else
        func = rand(funcs)
        # println("func: $func, ")
        args = Any[]
        for i in 1:arity(func)
            arg = genrandexpr(terms, funcs,
                              maxdepth=maxdepth-1, mindepth=mindepth-1, method=method)
            # println("arg_$i: $arg ")
            push!(args, arg)
        end
        Expr(:call, func, args...)
    end
    return expr
end

terminals = [:x, :rand, 1, 2, 3, 4]
functions = [:+, :-, :*, :/, :^]

ex = genrandexpr(terminals, functions, maxdepth=3, mindepth=2)
exfun = eval(Expr(:->, :x, ex))
exfun(1)

depth(ex) = isa(ex, Expr) ? maximum( depth(e) for e in ex.args )+1 : 0
children(ex) = isa(ex, Expr) ? sum( children(e) for e in ex.args) : 1

children(ex)
depth(ex)
dump(ex)

function subtree(root, idx)
    m, left, right = if !isa(root, Expr) || length(root.args) < 2
        1, root, root
    else
        left, right = root.args[2:3] # left & right nodes
        children(left)+1, left, right
    end
    if m == idx
        return root
    elseif m > idx
        return subtree(left, idx)
    else
        return subtree(right, idx - m)
    end
end
randsubtree(ex) = subtree(ex, rand(1:children(ex)))

[subtree(ex, i) for i in 1:children(ex)]

ext: jl


text:
using OrdinaryDiffEq, ForwardDiff, Statistics, LinearAlgebra

# which md curve to plot
which_dir = 1

## define dynamics of differential equation
function f(du, u, p, t)
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2] # prey
    du[2] = -p[3] * u[2] + p[4] * u[1] * u[2] # predator
end

u0 = [1.0;1.0] # initial populations
tspan = (0.0, 10.0) # time span to simulate over
t = collect(range(0, stop=10., length=200)) # time points to measure
p = [1.5,1.0,3.0,1.0] # initial parameter values
nom_prob = ODEProblem(f, u0, tspan, p) # package as an ODE problem
nom_sol = solve(nom_prob, Tsit5()) # solve

## Model features of interest are mean prey population, and max predator population (over time)
function features(p)
    prob = remake(nom_prob; p=p)
    sol = solve(prob, Vern9(); saveat=t)
    return [mean(sol[1,:]), maximum(sol[2,:])]
end

nom_features = features(p)

## loss function, we can take as l2 difference of features vs nominal features
function loss(p)
    prob = remake(nom_prob; p=p)
    p_features = features(p)
    loss = sum(abs2, p_features - nom_features)
    return loss
end

## gradient of loss function
function lossgrad(p, g)
    g[:] = ForwardDiff.gradient(p) do p
        loss(p)
    end
    return loss(p)
end

## package the loss and gradient into a DiffCost structure
cost = DiffCost(loss, lossgrad)

"""
We evaluate the hessian once only, at p.
Why? to find locally insensitive directions of parameter perturbation
The small eigenvalues of the Hessian are one easy way of defining these directions
"""
hess0 = ForwardDiff.hessian(loss, p)
ev(i) = -eigen(hess0).vectors[:,i]

init_dir = ev(which_dir); momentum = 1. ; span = (-15., 15.)
curve_prob = MDCProblem(cost, p, init_dir, momentum, span)

# rr = map(1:2) do i
# curve_prob_orig = curveProblem(cost, p, init_dir, momentum, span)

cb = [
    Verbose([CurveDistance(0.1:1:10), HamiltonianResidual(2.3:4:10)]),
    ParameterBounds([1,3], [-10.,-10.], [10.,10.])
    ]

cb = [
    Verbose([CurveDistance(0.1:1:10), HamiltonianResidual(2.3:4:10)])
    ]

# don't make the user do (curve_prob...) for Verbose
# make MomentumReadjustment

@time mdc = evolve(curve_prob, Tsit5; mdc_callback=cb);

# return cost_trajectory(mdc, mdc.sol.t) |> mean, cost_trajectory(mdc2, mdc2.sol.t) |> mean
# end

# function sol_at_p(p)
#     prob = remake(nom_prob; p=p)
#     sol = solve(prob, Tsit5())
# end

# p1 = plot(mdc; pnames=[L"p_1" L"p_2" L"p_3" L"p_4"])

# cost_vec = [mdc.cost(el) for el in eachcol(trajectory(mdc))]
# p2 = plot(distances(mdc), log.(cost_vec), ylabel="log(cost)", xlabel="distance", title="cost over MD curve");

# mdc_plot = plot(p1, p2, layout=(2, 1), size=(800, 800))

# nominal_trajectory = plot(sol_at_p(mdc(0.)[:states]), label=["prey" "predator"])
# perturbed_trajectory = plot(sol_at_p(mdc(-15.)[:states]), label=["prey" "predator"])
# traj_comparison = plot(nominal_trajectory, perturbed_trajectory, layout=(2, 1), xlabel="time", ylabel="population")

# Lessons


text:
@everywhere using Random
@everywhere using DataFrames
@everywhere using JSON
@everywhere using CSV
@everywhere using Combinatorics
@everywhere include("Gillespie.jl")

using Dates
function runif(min::Float64,max::Float64, n::Int64 = 1)
    r = (max - min).*rand(n) .+ min
    return r
end

function logunif(min::Float64,max::Float64, n::Int64 = 1)
    exps = (max - min).*rand(n) .+ min
    r = 10.0 .^ exps
    return r
end

function Sample_MoBlueParameters(total_time::Float64, n_sims::Int64, seed::Int64)
    Random.seed!(seed)
    seeds = rand(seed:(seed + n_sims), n_sims)
    total_mass = 1000000
    # You are trying a very extreme set of parameters, very low "stable_backward" and higher k_f, to see if this traps the mo36 structure
    sleep_time = logunif(-1.0, 1.0, n_sims)
    k_f = 0.3*1e-3 #logunif(-6.0,0.0, n_sims)
    k_d_stable = 0.01# logunif(-6.0,0.0, n_sims)
    dimerization_ratio = logunif(-3.0,2.0, n_sims)
    mo6_enhance = 0.0#logunif(-1.0, 2.0, n_sims)
    mo36_enhance = 100.0#logunif(0.0, 3.0, n_sims)
    #nano_stuct_enhance =0.0#
    nano_stuct_enhance = 100.0#logunif(1.0, 3.0, n_sims)
    @sync @distributed for i=1:n_sims
        sleep(sleep_time[i])
        #i_seed = Int64(floor(10.0^(10^sleep_time[i])))
        Run_MoBlue(seeds[i], total_mass, total_time, k_f, k_d_stable, dimerization_ratio[i], mo6_enhance, mo36_enhance, nano_stuct_enhance, nano_stuct_enhance)
    end

end

max_stars_repo_name: croningp/clusterautocatalyticsets

max_issues_repo_name: croningp/clusterautocatalyticsets

max_forks_repo_name: croningp/clusterautocatalyticsets


text:
using Ejemplo
using Test

@testset "Ejemplo.jl" begin
    @test f_xy(1,5) == 6# Write your tests here.
end


text:
struct Font
    surface::Cairo.CairoSurface
    charDict::Dict{Char, Tuple{Int, Int, Int, Int}}
    charSpacing::Int
    lineSpacing::Int

    Font(filename::String, fontString::String, charSpacing::Int=1, lineSpacing::Int=1) = new(open(Cairo.read_from_png, filename), generateFontDict(fontString), charSpacing, lineSpacing)
end

stringLineAlignOffsets = Dict{Symbol, Function}(
    :left => (widest, width) -> 0,
    :center => (widest, width) -> round(Int, (widest - width) / 2),
    :right => (widest, width) -> widest - width,
)

function getQuadTuple(row::Int, col::Int, width::Int=3, height::Int=5)
    return (8 * row - 8, 8 * col - 8, width, height)
end

# Use first occurance of a symbol as the target quad
function generateFontDict(fontString::String)
    fontMap = Dict{Char, Tuple{Int, Int, Int, Int}}()

    for (i, row) in enumerate(split(fontString, "\n"))
        for (j, c) in enumerate(row)
            get!(fontMap, c) do
                getQuadTuple(j, i, 3, 5)
            end
        end
    end

    return fontMap
end

# Multiline string
function lineWidth(font::Font, s::AbstractString, replaceMissing::Char=' ')
    return isempty(s) ? 0 : sum([get(font.charDict, c, font.charDict[replaceMissing])[3] for c in s]) + (length(s) - 1) * font.charSpacing
end

function lineHeight(font::Font, s::AbstractString, replaceMissing::Char=' ')
    return isempty(s) ? 0 : maximum([get(font.charDict, c, font.charDict[replaceMissing])[4] for c in s])
end

function Base.size(font::Font, s::AbstractString, replaceMissing::Char=' ')
    lines = split(s, "\n")

    width = maximum(lineWidth.(Ref(font), lines))
    height = sum(lineHeight.(Ref(font), lines)) + (length(lines) - 1) * font.lineSpacing

    return width, height
end

function drawString(ctx::Cairo.CairoContext, font::Font, s::AbstractString, x::Int=0, y::Int=0; replaceMissing::Char=' ', align::Symbol=:left, tint::Union{colorTupleType, Nothing}=colors.canvas_font_color)
    shouldSave = x != 0 || y != 0

    if shouldSave
        Cairo.save(ctx)
        Cairo.translate(ctx, x, y)
    end

    accX = 0
    accY = 0

    lineAllignFunc = get(stringLineAlignOffsets, align, :left)
    widest, height = size(font, s, replaceMissing)

    for line in split(s, "\n")
        width = lineWidth(font, line, replaceMissing)
        accX = lineAllignFunc(widest, width)

        for c in line
            quadX, quadY, width, height = get(font.charDict, c, font.charDict[replaceMissing])
            drawImage(ctx, font.surface, accX, accY, quadX, quadY, width, height, tint=tint)

            accX += width + font.charSpacing
        end

        accY += lineHeight(font, line, replaceMissing) + font.lineSpacing
    end

    if shouldSave
        Cairo.restore(ctx)
    end
end

const pico8FontString = raw"""

 !"#$% ´()*+,-./
0123456789:;<=>?
@abcdefghijklmno
pqrstuvwxyz[\]^_
`ABCDEFGHIJKLMNO
PQRSTUVWXYZ{|}~
"""

const pico8Font = missing


text:
Base.@irrational _sqrt3by4pi 0.4886025119029199 √(3/(4big(pi)))
Base.@irrational _sqrt3by2pi 0.690988298942671 √(3/(2big(pi)))

"""
    SphericalHarmonics.Pole <: Real

Supertype of [`NorthPole`](@ref) and [`SouthPole`](@ref).
"""
abstract type Pole <: Real end
"""
    SphericalHarmonics.NorthPole <: SphericalHarmonics.Pole

The angle ``\\theta = 0`` in spherical polar coordinates.
"""
struct NorthPole <: Pole end
"""
    SphericalHarmonics.SouthPole <: SphericalHarmonics.Pole

The angle ``\\theta = π`` in spherical polar coordinates.
"""
struct SouthPole <: Pole end

Base.cos(::NorthPole) = one(Float64)
Base.cos(::SouthPole) = -one(Float64)
Base.sin(::Pole) = zero(Float64)

for DT in [:zero, :one]
    @eval Base.$DT(::Type{<:Pole}) = $DT(Float64)
    @eval Base.$DT(x::Pole) = $DT(typeof(x))
end

Base.promote_rule(::Type{<:Pole}, T::Type{<:Real}) = promote_type(Float64, T)

# Return the value of θ corresponding to the poles
Base.Float64(::SouthPole) = Float64(pi)
Base.Float64(::NorthPole) = zero(Float64)

Base.float(p::Pole) = Float64(p)


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Visualizer for the (generic) DecisionTrees.jl
However, this visualizer is not generic. [FIXME]
Produces a TikzQTrees
"""
module DecisionTreeVis

export decisiontreevis, get_tree, get_metric

using DecisionTrees
using ExprSearch, DerivationTrees
using DataFrameSets
using TreeToJSON
using TikzQTrees
using RLESUtils, LatexUtils

get_tree() = error("get_tree not defined")
get_metric() = error("get_metric not defined")

function decisiontreevis{T}(dtree::DecisionTree, Dl::DFSetLabeled{T}, fileroot::AbstractString,
    limit_members::Int64, fmt_pretty::Format, fmt_natural::Format; plotpdf::Bool=true)

    get_depth(tree::DecisionTree) = get_depth(tree.root)
    get_depth(node::DTNode) = node.depth
    get_children(tree::DecisionTree) = get_children(tree.root)
    get_children(node::DTNode) = node.children

    get_name(tree::DecisionTree) = get_name(tree.root)
    function get_name(node::DTNode)
        id_text = "node id=" * string(node.id)
        members = getmeta(Dl, node.members)[:encounter_id] #show encounter numbers
        sort!(members)
        if length(members) <= limit_members
            members_text = "members=" * join(members, ",")
        else
            members_text = "members=" * join(members[1:limit_members], ",") *
                ", and $(length(members)-limit_members) more."
        end
        label = "label=$(node.label)"
        confidence = "confidence=" * string(signif(node.confidence, 3))
        if node.split_rule != nothing
            tree = get_tree(node.split_rule)
            expr = string(get_expr(tree))
            pretty = string("\\bf{", pretty_string(tree, fmt_pretty), "}")
            natural = pretty_string(tree, fmt_natural, true)
            score = "fitness (lower is better)=" * string(signif(get_metric(node.split_rule), 4))
            text = join([id_text, members_text, label, confidence, expr,
                pretty, natural, score], "\\\\")
        else
            text = join([id_text, members_text, label, confidence], "\\\\")
        end
        text
    end

    viscalls = VisCalls(get_name, get_children, get_depth)
    write_json(dtree, viscalls, "$(fileroot)_decisiontree.json")
    if plotpdf
        plottree("$(fileroot)_decisiontree.json", level_dist_cm=5,
            outfileroot="$(fileroot)_decisiontree")
    end
end

#these probably don't belong here
get_tree(result::SearchResult) = result.tree
get_metric(result::SearchResult) = result.fitness

end #module


text:
@testset "Basic polynomial functionalities" begin
    @testset "convert monomial to polynomial" begin
        x = semi_elementary_monomial((1,2,3),(0,0,0))
        p = to_polynomial(x)
        s = semi_elementary_polynomial(3)
        push!(s,semi_elementary_monomial((0,1,2),(0,0,1)),1)
        @test s == p
    end
end

max_stars_repo_name: yushoteke/SymmetricPolynomials

max_issues_repo_name: yushoteke/SymmetricPolynomials

max_forks_repo_name: yushoteke/SymmetricPolynomials


text:
module Supergrid

export runmodel, buildmodel, readresults, saveresults, analyzeresults, listresults, loadresults, makesets, makeparameters,
        fix_timezone_error, chart_energymix_scenarios

using JuMP, CPLEX, Gurobi, GLPKMathProgInterface, GLPK, Clp, Parameters, AxisArrays, Plots, JLD2, Statistics

include("helperfunctions.jl")
include("types.jl")
include("inputdata.jl")
include("builddispatch.jl")
include("jumpmodel.jl")
include("dispatchmodel.jl")
include("output.jl")
include("iewruns.jl")

defaultoptions() = Dict(
    :regionset => :Europe8,             # :Eurasia21, :Europe8
    :inputdatasuffix => "",             # e.g. "_landx2" to read solar input data "GISdata_solar2018_Europe8_landx2.mat"
    :runname => "",                     # change the run name without changing run parameters (e.g. if you modify the code)
    :islandindexes => [],               # superregion groupings, defaults to [1:8, 9:15, 16:21] for eurasia21, [] for europe8
    :carbontax => 0.0,                  # €/ton CO2
    :carboncap => 1.0,                  # global cap in kg CO2/kWh elec  (BAU scenario: ~0.5 kgCO2/kWh elec)
    :discountrate => 0.05,
    :maxbioenergy => 0.05,              # max share of biofuel of annual regional electricity demand (assuming CCGT, less if GT)
    :nuclearallowed => true,
    :globalnuclearlimit => Inf,         # maximum total nuclear capacity in all regions (GW)
    :hydroinvestmentsallowed => false,
    :transmissionallowed => :all,       # :none, :islands, :all
    :hours => 1,                        # 1,2,3 or 6 hours per period
    :solarwindarea => 1,                # area multiplier for GIS solar & wind potentials
    :datayear => 2017,                  # year of the ERA5 input data (produced by GlobalEnergyGIS.jl)
    :selectdays => 1,
    :skipdays => 0,
    :solver => :cplex,
    :threads => 4,
    :showsolverlog => true,
    :rampingconstraints => false,
    :rampingcosts => false,
    :disabletechs => [],
    :disableregions => [],
    :sspscenario => "ssp2-34",          # SSP scenario for synthetic demand
    :sspyear => 2050,                   # SSP year for synthetic demand
    :datafolder => "",                  # Full path to GIS input data. Set to "" to use the folder in HOMEDIR/.GlobalEnergyGIS_config.
    :resultsfile => "results.jld2"      # use "" to skip saving the results in the database
)

function autorunname(options)
    name = ""
    for (key,value) in setdiff(options, defaultoptions())
        name *= "$key=$value, "
    end
    if isempty(name)
        name = "default"
    else
        name = name[1:end-2]
    end
    name
end

function buildmodel(; optionlist...)
    options, hourinfo, sets, params = buildsetsparams(; optionlist...)
    return buildvarsmodel(options, hourinfo, sets, params)
end

function buildsetsparams(; optionlist...)
    println("\nReading input data...")
    options = merge(defaultoptions(), optionlist)
    hourinfo = HourSampling(options)
    @time sets = makesets(hourinfo, options)
    @time params = makeparameters(sets, options, hourinfo)
    return options, hourinfo, sets, params
end

function buildvarsmodel(options, hourinfo, sets, params)
    println("\nBuilding model...")
    modelname = initjumpmodel(options)
    print("  - variables:   ")
    @time vars = makevariables(modelname, sets)
    print("  - extra bounds:")
    @time setbounds(sets, params, vars, options)
    print("  - constraints: ")
    @time constraints = makeconstraints(modelname, sets, params, vars, hourinfo, options)
    print("  - objective:   ")
    @time makeobjective(modelname, sets, vars)

    return ModelInfo(modelname, sets, params, vars, constraints, hourinfo, options)
end

# BASIC USAGE: (carbon tax 50 €/ton CO2, 3-hour time periods)
# m, annualelec, capac, tcapac, chart = runmodel(carboncap=50, hours=3, [more options]...);
function runmodel(; name="", group="", optionlist...)       # carbon tax in €/ton CO2
    model = buildmodel(; optionlist...)

    #writeMPS(model, "model3.mps")
    if model.options[:solver] == :cplex
        println("\nSolving model using CPLEX version $(CPLEX.version())...")
    else
        println("\nSolving model...")
    end

    status = solve(model.modelname)
    println("\nSolve status: $status")

    println("\nReading results...")
    results = readresults(model, status)

    filename = model.options[:resultsfile]

    if !isempty(filename)
        if isempty(name)
            name = autorunname(model.options)
        end
        println("\nSaving results to disk...")
        saveresults(results, name, resultsfile=filename, group=group)
    end

    annualelec, capac, tcapac, chart = analyzeresults(results)

    if status != :Optimal
        @warn "The solver did not report an optimal solution. It could still be fine, but examine the log."
    end

    return results, annualelec, capac, tcapac, chart
end

end #module

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""Contains example usage of the features in the OpticSim.jl package."""
module Examples
using ..OpticSim
using ..OpticSim.Vis
using ..OpticSim.Geometry
using ..OpticSim.Emitters
using ..OpticSim.GlassCat
using ..OpticSim.Repeat

using StaticArrays
using DataFrames: DataFrame
using Images
using Unitful
using Plots
using LinearAlgebra
import Luxor

include("docs_examples.jl")
include("other_examples.jl")
include("repeating_structure_examples.jl")
include("eyemodels.jl")

end #module Examples
export Examples


text:
function solve(nums, use_growing=false)
    counter = zeros(Int, maximum(nums) + 1)
    grow = copy(counter)
    for num in nums
        counter[num + 1] += 1
    end
    s = 0
    while length(counter) > 1
        front, back = counter[1] + grow[1], counter[end] + grow[end]
        if front < back
            counter[2] += counter[1]
            popfirst!(counter)
            if use_growing
                popfirst!(grow)
                grow[1] = front
            end
            s += front
        else
            counter[end-1] += counter[end]
            pop!(counter)
            if use_growing
                pop!(grow)
                grow[end] = back
            end
            s += back
        end
    end
    return s
end

function main()
    nums = (x -> parse(Int, x)).(split(readlines()[1], ','))
    println(solve(nums, false))
    println(solve(nums, true))
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
using Test
using OptionalArgChecks

@testset "argcheck" begin
    function f(x)
        @mark argcheck x<2 && if x == 1
            return "foo"
        else
            error("Test")
        end
        return x + 3
    end

    for x in 0:2
        @test @skip argcheck f(x) == x + 3
    end

    g(x) = @argcheck x > 0

    @test_throws ArgumentError g(-1)
    @test @unsafe_skipargcheck(g(-1)) === nothing
    @test_throws ArgumentError @skip(argcheck, g(-1)) === nothing
    @test g(1) === nothing

    g(x) = @check x > 0

    @test_throws Exception g(-1)
    @test @unsafe_skipargcheck(g(-1)) === nothing
    @test_throws Exception @skip(argcheck, g(-1)) === nothing
    @test g(1) === nothing

    outer(x) = inner(x)
    inner(x) = (@argcheck x; return x)
    @test @unsafe_skipargcheck(inner(true)) == true
    @test @unsafe_skipargcheck(inner(false)) == false
    @test @unsafe_skipargcheck(inner(false), recursive=true) == false
    @test @unsafe_skipargcheck(inner(false), recursive=false) == false

    @test @unsafe_skipargcheck(outer(true)) == true
    @test @unsafe_skipargcheck(outer(false)) == false
    @test @unsafe_skipargcheck(outer(false), recursive=true) == false
    @test_throws ArgumentError @unsafe_skipargcheck(outer(false), recursive=false)
end

@testset "mark skip" begin
    function simple()
        @mark return1 begin
            return 1
        end
        return 2
    end

    @test simple() === 1
    @test @skip(does_not_exist, simple()) === 1
    @test @skip(return1, simple()) === 2

    @test_throws LoadError eval(:(@skip(2 + 3, simple())))

    function indirect()
        simple()
    end

    @test indirect() === 1
    @test @skip(does_not_exist, indirect()) === 1
    @test @skip(return1, indirect()) === 2

    function complex()
        ret = Int[]
        push!(ret,1)
        @mark two begin
            push!(ret,2)
        end
        push!(ret, 3)
        @mark four begin
            push!(ret, 4)
        end
        ret
    end

    @test complex() == 1:4
    @test @skip(does_not_exist, complex()) == 1:4
    @test @skip(two, complex()) == [1,3,4]
    @test @skip(four, complex()) == [1,2,3]
    @test @skip([two, four], complex()) == [1,3]

    function nested()
        ret = Int[]
        @mark nested begin
            push!(ret, 1)
            @mark nested begin
                push!(ret,2)
            end
            push!(ret,3)
        end
        push!(ret, 4)
        ret
    end

    @test nested() == 1:4
    @test @skip(nested, nested()) == [4]

    function complex_nested()
        ret = Char[]
        push!(ret,'a')
        @mark bcdg begin
            push!(ret,'b')
            @mark cdh begin
                push!(ret,'c')
                @mark bcdg begin
                    push!(ret,'d')
                end
            end
        end
        push!(ret,'e')
        @mark fgh begin
            push!(ret,'f')
            @mark bcdg begin
                push!(ret,'g')
            end
            @mark cdh begin
                push!(ret,'h')
            end
        end
        ret
    end

    @test @skip([], complex_nested()) == 'a':'h'
    @test @skip([bcdg], complex_nested()) == collect("aefh")
    @test @skip([cdh], complex_nested()) == collect("abefg")
    @test @skip([fgh], complex_nested()) == 'a':'e'
    @test @skip([bcdg, cdh], complex_nested()) == collect("aef")
    @test @skip([fgh, cdh], complex_nested()) == collect("abe")
    @test @skip([bcdg, cdh, fgh], complex_nested()) == collect("ae")
end

@testset "recursion through function calls" begin
    function inner()
        ret = Int[]
        push!(ret, 2)
        @mark three begin
            push!(ret, 3)
        end
        ret
    end

    function outer()
        ret = Int[]
        push!(ret, 1)
        x = inner()
        append!(ret, x)
        @mark four begin
            push!(ret, 4)
        end
        push!(ret, 5)
        ret
    end

    @test outer() == 1:5
    @test @skip(four, outer(), recursive=true) == [1,2,3,5]
    @test @skip(four, outer(), recursive=false) == [1,2,3,5]
    @test @skip(three, outer(), recursive=true) == [1,2,4,5]
    @test @skip(three, outer(), recursive=false) == 1:5
    @test @skip([three, four], outer()) == [1,2,5]
    @test @skip([three, four], outer(), recursive=false) == [1,2,3,5]
end

using Documenter
DocMeta.setdocmeta!(OptionalArgChecks, :DocTestSetup, :(using OptionalArgChecks); recursive=true)
doctest(OptionalArgChecks; manual = false)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct TopologyParameters <: SwaggerModel
    targetResourceGroupName::Any # spec type: Union{ Nothing, String } # spec name: targetResourceGroupName

    function TopologyParameters(;targetResourceGroupName=nothing)
        o = new()
        validate_property(TopologyParameters, Symbol("targetResourceGroupName"), targetResourceGroupName)
        setfield!(o, Symbol("targetResourceGroupName"), targetResourceGroupName)
        o
    end
end # type TopologyParameters

const _property_map_TopologyParameters = Dict{Symbol,Symbol}(Symbol("targetResourceGroupName")=>Symbol("targetResourceGroupName"))
const _property_types_TopologyParameters = Dict{Symbol,String}(Symbol("targetResourceGroupName")=>"String")
Base.propertynames(::Type{ TopologyParameters }) = collect(keys(_property_map_TopologyParameters))
Swagger.property_type(::Type{ TopologyParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_TopologyParameters[name]))}
Swagger.field_name(::Type{ TopologyParameters }, property_name::Symbol) =  _property_map_TopologyParameters[property_name]

function check_required(o::TopologyParameters)
    (getproperty(o, Symbol("targetResourceGroupName")) === nothing) && (return false)
    true
end

function validate_property(::Type{ TopologyParameters }, name::Symbol, val)
end


text:
# this script does a Monte Carlo that applies IV to solve inconsistency of OLS estimator
# when there is measurement error of the regressors
using Econometrics, Distributions, LinearAlgebra, Random

function GIVmoments(theta, y, x, inst)
	e = y - x*theta
	ms = e.*inst
end

function main()
# do the Monte Carlo
n = 100
sig = 1
reps = 1000
results = zeros(reps,3)
ystar = zeros(n)
for rep = 1:reps
	x = randn(n) # an exogenous regressor
	e = randn(n) # the error term
	# generate the dep var
	for t = 2:n
	  ystar[t] = 0.0 + 0.9*ystar[t-1] + 1.0*x[t] + e[t]
	end
    # add measurement error
	y = ystar + sig*randn(n)
	# now do GMM, using the data with meas. error in both dep. var. and regressor
	data = [y lag(y,1) x lags(x,2)];
    data = data[3:end,:] # drop first obsns, missing due to lags
    y = data[:,1]
    ylag = data[:,2]
    x = data[:,3]
    xlags = data[:,4:end]
    regressors = [ones(size(y)) ylag x]
	instruments = [ones(size(y)) x xlags]
    xhat = instruments*(instruments\regressors)
    thetahat = xhat\y
    ms = GIVmoments(thetahat, y, regressors, instruments)
    weight = inv(cov(ms))
    moments = theta-> GIVmoments(theta, y, regressors, instruments)
	junk, objvalue, junk, junk, junk = gmm(moments, thetahat, weight)
   	objvalue = size(y,1)*objvalue
	df = size(ms,2) - size(thetahat,1)
	reject = objvalue .> quantile.(Ref(Chisq(df)),[0.9 0.95 0.99]) # does the test reject at these signif levels?
    results[rep,:] = reject
end
cnames = ["10%", "5%", "1%"]
println("rejection frequencies, nominal 10% 5% and 1%:")
prettyprint(mean(results,dims=1),cnames)
return
end
main()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
import Distributions: logpdf
import Random: rand
import Base: length
export Gaussian, normal_logpdf, StandardGaussian, GaussianPair, logpdf
export standard_to_natural, natural_to_standard
export standard_to_expectation, expectation_to_standard
export expectation_to_natural, natural_to_expectation
export standard_to_unconstrained, unconstrained_to_standard
export natural_to_unconstrained, unconstrained_to_natural
export expectation_to_unconstrained, unconstrained_to_expectation

"""
    Gaussian{Tm, Tchol_Σ<:Cholesky}

A (multivariate) Gaussian with mean vector `m` and cholesky factorisation `chol_Σ` of the
covariance matrix.
"""
struct Gaussian{Tm<:AbstractVector{<:Real}, Tchol_Σ<:Cholesky}
    m::Tm
    chol_Σ::Tchol_Σ
end

length(x::Gaussian) = length(x.m)

rand(rng::AbstractRNG, x::Gaussian) = x.m + x.chol_Σ.U' * randn(rng, length(x.m))
function rand(rng::AbstractRNG, x::Gaussian, S::Int)
    return x.m .+ x.chol_Σ.U' * randn(rng, length(x.m), S)
end

# Returns a `Vector` of length `size(Y, 2)` containing the logpdf of each column of `Y`
# under `x`.
function logpdf(x::Gaussian, Y::AbstractMatrix{<:Real})
    @assert length(x) == size(Y, 1)
    tmp = length(x) * log(2π) + logdet(x.chol_Σ)
    return .-(tmp .+ diag_Xt_invA_X(x.chol_Σ, Y .- x.m)) / 2
end
logpdf(x::Gaussian, y::AbstractVector{<:Real}) = logpdf(x, reshape(y, :, 1))[1]

# A vectorised logpdf for the normal distribution because Zygote's broadcasting sucks
# if you don't do things correctly :S Even differnetiating through this kind of sucks,
# which is a real shame.
function normal_logpdf(y, m, σ)
    return .-(log(2π) .+ 2 .* log.(σ) .+ ((y .- m) ./ σ).^2) ./ 2
end

"""
    StandardGaussian

`D`-dimensional multivariate Gaussian with zero mean and identity covariance.
"""
struct StandardGaussian
    D::Int
end

length(X::StandardGaussian) = X.D

rand(rng::AbstractRNG, X::StandardGaussian) = randn(rng, length(X))
rand(rng::AbstractRNG, X::StandardGaussian, S::Int) = randn(rng, length(X), S)

function logpdf(X::StandardGaussian, y::AbstractVector{<:Real})
    return -(X.D * log(2π) + sum(abs2, y)) / 2
end

function logpdf(X::StandardGaussian, Y::AbstractMatrix{<:Real})
    return .-(X.D * log(2π) .+ vec(sum(abs2, Y;  dims=1))) / 2
end

"""
    GaussianPair{Tx1<:Gaussian, Tx2<:Gaussian}

A pair of independent Gaussians `x1` and `x2`.
"""
struct GaussianPair{Tx1<:Gaussian, Tx2<:Gaussian}
    x1::Tx1
    x2::Tx2
end

length(x::GaussianPair) = length(x.x1) + length(x.x2)

function rand(rng::AbstractRNG, x::GaussianPair, S::Int)
    return (x1=rand(rng, x.x1, S), x2=rand(rng, x.x2, S))
end
rand(rng::AbstractRNG, x::GaussianPair) = (x1=rand(rng, x.x1), x2=rand(rng, x.x2))

function logpdf(x::GaussianPair, Y::NamedTuple{(:x1, :x2)})
    return logpdf(x.x1, Y.x1) + logpdf(x.x2, Y.x2)
end

#
# Constrained parametrisations.
#

"""
    natural_to_standard(θ₁::AbstractVector{<:Real}, θ₂::AbstractMatrix{<:Real})

Convert from natural to standard parametrisation.
"""
function natural_to_standard(θ₁::AbstractVector{<:Real}, θ₂::Symmetric{<:Real})
    S = Symmetric(inv(cholesky(-θ₂)) ./ 2)
    return S * θ₁, S
end

# forwards-mode rule for natural_to_standard. Avoids some repeated computation.
function frule(
    ::typeof(natural_to_standard),
    θ₁::AbstractVector{<:Real},
    θ₂::Symmetric{<:Real},
)
    C = cholesky(-θ₂)
    S = Symmetric(inv(C) ./ 2)
    m = S * θ₁
    function natural_to_standard_pushforward(Δself, Δθ₁, Δθ₂)
        ΔS = Symmetric((C \ Matrix((C \ Δθ₂')')) ./ 2)
        Δm = ΔS * θ₁ + S * Δθ₁
        return (Δm, ΔS)
    end
    return (m, S), natural_to_standard_pushforward
end

"""
    standard_to_natural(m::AbstractVector{<:Real}, S::Symmetric{<:Real})

Convert from standard to natural parametrisation.
"""
function standard_to_natural(m::AbstractVector{<:Real}, S::Symmetric{<:Real})
    C = inv(cholesky(S))
    return C * m, Symmetric(C ./ (-2))
end

"""
    expectation_to_standard(η₁::AbstractVector{<:Real}, η₂::Symmetric{<:Real})

Convert from expectation to standard parametrisation.
"""
function expectation_to_standard(η₁::AbstractVector{<:Real}, η₂::Symmetric{<:Real})
    return η₁, η₂ - Symmetric(η₁ * η₁')
end

"""
    standard_to_expectation(m::AbstractVector{<:Real}, S::Symmetric{<:Real})

Convert from standard to natural parametrisation.
"""
function standard_to_expectation(m::AbstractVector{<:Real}, S::Symmetric{<:Real})
    return m, Symmetric(S + m * m')
end

function frule(
    ::typeof(standard_to_expectation),
    m::AbstractVector{<:Real},
    S::Symmetric{<:Real},
)
    function standard_to_expectation_pushforward(Δself, Δm, ΔS::Symmetric{<:Real})
        Δη₁ = Δm
        Δη₂ = Symmetric(ΔS + m * Δm' + Δm * m')
        return Δη₁, Δη₂
    end
    return standard_to_expectation(m, S), standard_to_expectation_pushforward
end

"""
    expectation_to_natural(η₁::AbstractVector{<:Real}, η₂::Symmetric{<:Real})

Convert from expectation to natural parameters.
"""
function expectation_to_natural(η₁::AbstractVector{<:Real}, η₂::Symmetric{<:Real})
    return standard_to_natural(expectation_to_standard(η₁, η₂)...)
end

"""
    natural_to_expectation(θ₁::AbstractVector{<:Real}, θ₂::Symmetric{<:Real})

Convert from natural to expectation parametrisation.
"""
function natural_to_expectation(θ₁::AbstractVector{<:Real}, θ₂::Symmetric{<:Real})
    m, S = natural_to_standard(θ₁, θ₂)
    η₁, η₂ = standard_to_expectation(m, S)
    return η₁, η₂
end

function frule(
    ::typeof(natural_to_expectation),
    θ₁::AbstractVector{<:Real},
    θ₂::Symmetric{<:Real},
)
    (m, S), nat_to_std_pushforward = frule(natural_to_standard, θ₁, θ₂)
    (η₁, η₂), std_to_exp_pushforward = frule(standard_to_expectation, m, S)
    function natural_to_expectation_pushforward(Δself, Δθ₁, Δθ₂::Symmetric)
        Δm, ΔS = nat_to_std_pushforward(DoesNotExist(), Δθ₁, Δθ₂)
        return std_to_exp_pushforward(DoesNotExist(), Δm, ΔS)
    end
    return (η₁, η₂), natural_to_expectation_pushforward
end

#
# Unconstrained parametrisations.
#

"""
    standard_to_unconstrained(m::AbstractVector{<:Real}, S::Symmetric{<:Real})

Convert from standard to mean and unconstrained upper-triangular matrix.
"""
function standard_to_unconstrained(m::AbstractVector{<:Real}, S::Symmetric{<:Real})
    return m, to_positive_definite_softplus_inv(chol(S))
end

function frule(
    ::typeof(standard_to_unconstrained),
    m::AbstractVector{<:Real},
    S::Symmetric{<:Real},
)
    # Perform forwards-pass and retain information required to compute pushforward.
    U, pushforward_U = frule(chol, S)
    V, pushforward_V = frule(to_positive_definite_softplus_inv, U)

    function standard_to_unconstrained_pushforward(Δself, Δm, ΔS::Symmetric{<:Real})
        ΔU = pushforward_U(DoesNotExist(), ΔS)
        ΔV = pushforward_V(DoesNotExist(), ΔU)
        return (Δm, ΔV)
    end

    return (m, V), standard_to_unconstrained_pushforward
end

"""
    unconstrained_to_standard(m::AbstractVector{<:Real}, U::UpperTriangular{<:Real})

Convert from mean + unconstrained upper-triangular matrix to standard parametrisation.
"""
function unconstrained_to_standard(m::AbstractVector{<:Real}, U::UpperTriangular{<:Real})
    UA = to_positive_definite_softplus(U)
    return m, Symmetric(UA'UA)
end

"""
    natural_to_unconstrained(θ₁::AbstractVector{<:Real}, θ₂::AbstractMatrix{<:Real})

Convert from natural to mean and unconstrained upper-triangular matrix.
"""
function natural_to_unconstrained(θ₁::AbstractVector{<:Real}, θ₂::Symmetric{<:Real})
    return standard_to_unconstrained(natural_to_standard(θ₁, θ₂)...)
end

function frule(
    ::typeof(natural_to_unconstrained),
    θ₁::AbstractVector{<:Real},
    θ₂::Symmetric{<:Real},
)
    (m_, S), pushforward_standard = frule(natural_to_standard, θ₁, θ₂)
    (m, U), pushforward_unconstrained = frule(standard_to_unconstrained, m_, S)

    function pushforward_natural_to_unconstrained(Δself, Δθ₁, Δθ₂::Symmetric{<:Real})
        Δm_, ΔS = pushforward_standard(DoesNotExist(), Δθ₁, Δθ₂)
        return pushforward_unconstrained(DoesNotExist(), Δm_, ΔS)
    end
    return (m, U), pushforward_natural_to_unconstrained
end

"""
    unconstrained_to_natural(m::AbstractVector{<:Real}, U::AbstractMatrix{<:Real})

Convert from mean and unconstrained upper-triangular matrix to natural.
"""
function unconstrained_to_natural(m::AbstractVector{<:Real}, U::UpperTriangular{<:Real})
    return standard_to_natural(unconstrained_to_standard(m, U)...)
end

"""
    expectation_to_unconstrained(η₁::AbstractVector{<:Real}, η₂::AbstractMatrix{<:Real})

Convert from expecation parameters to mean and unconstrained upper-triangular matrix.
"""
function expectation_to_unconstrained(
    η₁::AbstractVector{<:Real},
    η₂::AbstractMatrix{<:Real},
)
    return standard_to_unconstrained(expectation_to_standard(η₁, η₂)...)
end

"""
    unconstrained_to_expectation(η₁::AbstractVector{<:Real}, η₂::AbstractMatrix{<:Real})

Convert from expecation parameters to mean and unconstrained upper-triangular matrix.
"""
function unconstrained_to_expectation(
    m::AbstractVector{<:Real},
    U::UpperTriangular{<:Real},
)
    return standard_to_expectation(unconstrained_to_standard(m, U)...)
end


text:
struct ProductTable{names, T <: NamedTuple{names}, T1 <: AbstractVector{<:NamedTuple}, T2 <: AbstractVector{<:NamedTuple}} <: AbstractMatrix{T}
    t1::T1
    t2::T2
    # TODO: Extra acceleration indices, like foreign keys / 1-to-1 relationships / etc
end

function ProductTable(t1::AbstractVector{NamedTuple{names1,T1}}, t2::AbstractVector{NamedTuple{names2,T2}}) where {names1, T1, names2, T2}
    names = (names1..., names2...)
    T = NamedTuple{names, _cat_types(T1, T2)}
    return ProductTable{names, T, typeof(t1), typeof(t2)}(t1, t2)
end

size(table::ProductTable) = (size(table.t1)..., size(table.t2)...)
axes(table::ProductTable) = (axes(table.t1)..., axes(table.t2)...)
@propagate_inbounds function getindex(table::ProductTable, i1::Int, i2::Int)
    return merge(table.t1[i1], table.t2[i2])
end

colnames(::ProductTable{names}) where {names} = names
#columns(t::ProductTable{names}) where {names} = NamedTuple{names}(map(name ->)
#accelerators(t::ProductTable) = (accelerators(t.t1)..., accelerators(t.t2)...)

@inline function (p::Project{names})(t::ProductTable{<:Any, <:Any, <:AbstractVector{<:NamedTuple{n1}}, <:AbstractVector{<:NamedTuple{n2}}}) where {n1, n2, names}
    p1 = Project{_intersect(Val(names), Val(n1))}()
    p2 = Project{_intersect(Val(names), Val(n2))}()
    # TODO ensure no other names left over
    ProductTable(p1(t.t1), p2(t.t2))
end

@inline function (r::Rename{oldnames, newnames})(t::ProductTable{names}) where {oldnames, newnames, names}
    names2 = _rename(Val(oldnames), Val(newnames), Val(names))
    return ProductTable{names2}(r(t.t1), r(t.t2))
end

cross(t1::AbstractVector{<:NamedTuple}, t2::AbstractVector{<:NamedTuple}) = ProductTable(t1, t2)

function similar(t::ProductTable, ::Type{NamedTuple{names, Ts}}, dims::Tuple{Int}) where {names, Ts}
    data = _makevectors(Ts, dims)
    return Table{names}(data)
end

# Filter - seperate case where predicate applies just to one part of the product table
function filter(pred::Predicate{names}, t::ProductTable{<:Any, <:Any, <:AbstractVector{<:NamedTuple{n1}}, <:AbstractVector{<:NamedTuple{n2}}}) where {names, n1, n2}
    if _issubset(names, n1)
        return ProductTable(filter(pred, t.t1), t.t2)
    elseif _issubset(names, n2)
        return ProductTable(t.t1, filter(pred, t.t2))
    end
    return @inbounds t[findall(pred, t)] # Don't use map - requires N^2 memory and most joins are sparse
end

# Map

# First layer of dispatch seperates cases where predicate applies just to one part of the product table
function map(pred::Predicate{names}, t::ProductTable{<:Any, <:Any, <:AbstractVector{<:NamedTuple{n1}}, <:AbstractVector{<:NamedTuple{n2}}}) where {names, n1, n2}
    if _issubset(names, n1)
        return ProductArray((x, y) -> x, map(pred, t.t1), Array{Nothing}(undef, size(t.t2)))
    elseif _issubset(names, n2)
        return ProductArray((x, y) -> y, Array{Nothing}(undef, size(t.t1)), map(pred, t.t2))
    end

    n1_projected = _intersect(Val(names), Val(n1))
    n2_projected = _intersect(Val(names), Val(n2))
    names_projected = (n1_projected..., n2_projected...)

    t_projected = project(t, names_projected)
    index1 = promote_index(project(accelerators(t.t1), n1_projected)...)
    index2 = promote_index(project(accelerators(t.t2), n2_projected)...)

    return @inbounds _map(pred, t_projected, index1, index2)
end

# Seperate the cases where we know we can or can't do some acceleration with the predicate
function _map(pred::Predicate, t::ProductTable, ::NoIndex, ::NoIndex)
    out = Array{Bool}(undef, size(t))

    @inbounds for i in keys(t)
        out[i] = pred(t[i])
    end

    return out
end

# Default to applying the accelerated filter to the second table
function _map(pred::Predicate, t::ProductTable, ::Any, ::Any)
    out = Array{Bool}(undef, size(t))

    @inbounds for i in keys(t.t1)
        x = t.t1[i]
        tmp = map(Predicate(pred, x), t.t2)
        out[i, :] = tmp
    end

    return out
end

# If the second table has no acceleration index, try switch the ordering (TODO is this a good idea? Should this be more predictable to the user?)
function _map(pred::Predicate, t::ProductTable, ::Any, ::NoIndex) # Reverse looping order
    out = Array{Bool}(undef, size(t))

    @inbounds for i in keys(t.t2)
        x = t.t2[i]
        tmp = map(Predicate(pred, x), t.t1)
        out[:, i] = tmp
    end

    return out
end

# TODO implement a sort-merge join algorithm:

# function _map(pred::IsEqual, t::ProductTable, ::SortIndex, ::SortIndex)
#     ...
# end

# findall

# First layer of dispatch seperates cases where predicate applies just to one part of the product table
function findall(pred::Predicate{names}, t::ProductTable{<:Any, <:Any, <:AbstractVector{<:NamedTuple{n1}}, <:AbstractVector{<:NamedTuple{n2}}}) where {names, n1, n2}
    if _issubset(names, n1)
        return ProductArray((x, y) -> CartesianIndex(x, y), findall(pred, t.t1), keys(t.t2))
    elseif _issubset(names, n2)
        return ProductArray((x, y) -> CartesianIndex(x, y), keys(t.t1), findall(pred, t.t2))
    end

    n1_projected = _intersect(Val(names), Val(n1))
    n2_projected = _intersect(Val(names), Val(n2))
    names_projected = (n1_projected..., n2_projected...)

    t_projected = project(t, names_projected)
    index1 = promote_index(project(accelerators(t.t1), n1_projected)...)
    index2 = promote_index(project(accelerators(t.t2), n2_projected)...)

    return @inbounds _findall(pred, t_projected, index1, index2)
end

# Seperate the cases where we know we can or can't do some acceleration with the predicate
function _findall(pred::Predicate, t::ProductTable, ::NoIndex, ::NoIndex)
    return findall(row -> pred(row), t)
end

# Default to applying the accelerated filter to the second table
function _findall(pred::Predicate, t::ProductTable, ::Any, index2::Any)
    out = Vector{CartesianIndex{2}}()

    @inbounds for i in keys(t.t1)
        x = t.t1[i]
        tmp = _findall(Predicate(pred, x), t.t2, index2)
        for j in tmp # TODO: use a lazily-mapped append!, or something
            push!(out, CartesianIndex(i, j))
        end
    end

    return out
end

# If the second table has no acceleration index, try switch the ordering (TODO is this a good idea? Should this be more predictable to the user?)
function _findall(pred::Predicate, t::ProductTable, index1::Any, ::NoIndex) # Reverse looping order
    out = Vector{CartesianIndex{2}}()

    @inbounds for j in keys(t.t2)
        x = t.t2[j]
        tmp = map(Predicate(pred, x), t.t1)
        for i in tmp # TODO: use a lazily-mapped append!, or something
            push!(out, CartesianIndex(i, j))
        end
    end

    return out
end

# Sort-merge join algorithm:

function _findall(pred::Equals, t::ProductTable, index1::SortIndex, index2::SortIndex)
    out = Vector{CartesianIndex{2}}()
    if isempty(t)
        return out
    end

    @inbounds t1 = view(getdata(t.t1)[1], index1.order) # a single column.
    n1 = length(t1)
    i1 = 1
    @inbounds x1 = t1[1]

    @inbounds t2 = view(getdata(t.t2)[1], index2.order) # a single column.
    n2 = length(t2)
    i2 = 1
    @inbounds x2 = t2[1]

    while true
        if isequal(x1, x2)
            i1_last = searchsortedlast(t1, x1)
            i2_last = searchsortedlast(t2, x2)
            for j1 in i1:i1_last
                for j2 in i2:i2_last
                    @inbounds push!(out, CartesianIndex(index1.order[j1], index2.order[j2]))
                end
            end
            i1 = i1_last + 1; if i1 > n1; break; end
            i2 = i2_last + 1; if i2 > n2; break; end
            @inbounds x1 = t1[i1]
            @inbounds x2 = t2[i2]
        elseif isless(x1, x2)
            i1 = searchsortedfirst(t1, x2); if i1 > n1; break; end
            @inbounds x1 = t1[i1]
        else
            i2 = searchsortedfirst(t2, x1)
            @inbounds x2 = t2[i2]
        end
    end

    return out
end

function _findall(pred::Equals, t::ProductTable, index1::UniqueSortIndex, index2::SortIndex)
    out = Vector{CartesianIndex{2}}()
    if isempty(t)
        return out
    end

    @inbounds t1 = view(getdata(t.t1)[1], index1.order) # a single column.
    n1 = length(t1)
    i1 = 1
    @inbounds x1 = t1[1]

    @inbounds t2 = view(getdata(t.t2)[1], index2.order) # a single column.
    n2 = length(t2)
    i2 = 1
    @inbounds x2 = t2[1]

    while true
        if isequal(x1, x2)
            i2_last = searchsortedlast(t2, x2)
            for j2 in i2:i2_last
                @inbounds push!(out, CartesianIndex(index1.order[i1], index2.order[j2]))
            end
            i1 = i1 + 1; if i1 > n1; break; end
            i2 = i2_last + 1; if i2 > n2; break; end
            @inbounds x1 = t1[i1]
            @inbounds x2 = t2[i2]
        elseif isless(x1, x2)
            i1 = searchsortedfirst(t1, x2); if i1 > n1; break; end
            @inbounds x1 = t1[i1]
        else
            i2 = searchsortedfirst(t2, x1); if i2 > n2; break; end
            @inbounds x2 = t2[i2]
        end
    end

    return out
end

function _findall(pred::Equals, t::ProductTable, index1::SortIndex, index2::UniqueSortIndex)
    out = Vector{CartesianIndex{2}}()
    if isempty(t)
        return out
    end

    @inbounds t1 = view(getdata(t.t1)[1], index1.order) # a single column.
    n1 = length(t1)
    i1 = 1
    @inbounds x1 = t1[1]

    @inbounds t2 = view(getdata(t.t2)[1], index2.order) # a single column.
    n2 = length(t2)
    i2 = 1
    @inbounds x2 = t2[1]

    while true
        if isequal(x1, x2)
            i1_last = searchsortedlast(t1, x1)
            for j1 in i1:i1_last
                @inbounds push!(out, CartesianIndex(index1.order[j1], index2.order[i2]))
            end
            i1 = i1_last + 1; if i1 > n1; break; end
            i2 = i2 + 1; if i2 > n2; break; end
            @inbounds x1 = t1[i1]
            @inbounds x2 = t2[i2]
        elseif isless(x1, x2)
            i1 = searchsortedfirst(t1, x2); if i1 > n1; break; end
            @inbounds x1 = t1[i1]
        else
            i2 = searchsortedfirst(t2, x1); if i2 > n2; break; end
            @inbounds x2 = t2[i2]
        end
    end

    return out
end

function _findall(pred::Equals, t::ProductTable, index1::UniqueSortIndex, index2::UniqueSortIndex)
    out = Vector{CartesianIndex{2}}()
    if isempty(t)
        return out
    end

    @inbounds t1 = view(getdata(t.t1)[1], index1.order) # a single column.
    n1 = length(t1)
    i1 = 1
    @inbounds x1 = t1[1]

    @inbounds t2 = view(getdata(t.t2)[1], index2.order) # a single column.
    n2 = length(t2)
    i2 = 1
    @inbounds x2 = t2[1]

    while true
        if isequal(x1, x2)
            @inbounds push!(out, CartesianIndex(index1.order[i1], index2.order[i2]))
            i1 = i1 + 1; if i1 > n1; break; end
            i2 = i2 + 1; if i2 > n2; break; end
            @inbounds x1 = t1[i1]
            @inbounds x2 = t2[i2]
        elseif isless(x1, x2)
            i1 = searchsortedfirst(t1, x2); if i1 > n1; break; end
            @inbounds x1 = t1[i1]
        else
            i2 = searchsortedfirst(t2, x1); if i2 > n2; break; end
            @inbounds x2 = t2[i2]
        end
    end

    return out
end


text:
using LaplacianOpt
using Test

import Memento
import JuMP
import LinearAlgebra
import GLPK
import MathOptInterface

const LOpt = LaplacianOpt
const LA   = LinearAlgebra
const MOI  = MathOptInterface

# Suppress warnings during testing
LOpt.logger_config!("error")

glpk_optimizer = JuMP.optimizer_with_attributes(GLPK.Optimizer, MOI.Silent() => true)

@testset "LaplacianOpt" begin

    include("utility_tests.jl")
    include("lo_model_tests.jl")
    include("relaxations_test.jl")
    include("log_tests.jl")

end


text:
import Clang
import HiGHS_jll

const LIBHIGHS_HEADERS = [
    joinpath(HiGHS_jll.artifact_dir, "include", "interfaces", "highs_c_api.h"),
]

const GEN_DIR = joinpath(@__DIR__, "..", "src", "wrapper")
if !isdir(GEN_DIR)
    mkdir(GEN_DIR)
end

wc = Clang.init(
    headers = LIBHIGHS_HEADERS,
    output_file = joinpath(GEN_DIR, "libhighs_api.jl"),
    common_file = joinpath(GEN_DIR, "libhighs_common.jl"),
    header_wrapped = (root, current) -> root == current,
    header_library = x -> "libhighs",
    clang_diagnostics = true,
)

run(wc)

rm(joinpath(GEN_DIR, "LibTemplate.jl"))


text:
export AbstractSystem, AbstractSystemCache,
    AbstractHomotopy, AbstractHomotopyCache

"""
    AbstractSystem

Representing a system of polynomials.
"""
abstract type AbstractSystem end

"""
    AbstractSystemCache

A cache to avoid allocations for the evaluation of an [`AbstractSystem`](@ref).
"""
abstract type AbstractSystemCache end

"""
    AbstractHomotopy

Representing a homotopy.
"""
abstract type AbstractHomotopy end

"""
    AbstractHomotopyCache

A cache to avoid allocations for the evaluation of an [`AbstractHomotopy`](@ref).
"""
abstract type AbstractHomotopyCache end

include("systems.jl")
include("homotopies.jl")


text:
using CSV, DataFrames
using FluxDash
using FluxDash.Flux
using Random

df = CSV.read("$(@__DIR__)/../datasets/real_estate.csv", DataFrame)

n_row = nrow(df)
df = df[shuffle(1:n_row)[:], :] #shuffle data

labels = names(df)
in_labels = labels[end]
out_labels = labels[end]

X_train, X_test, y_train, y_test = FluxDash.format_nn_data(
    df, in_labels, out_labels; training_percent=75
)

chn = Base.Channel{Tuple{Float64, Int, Chain}}(Inf)

hidden_outs = [2]
hidden_activations =["relu"]
FluxDash.train_nn(
    X_train, y_train, in_labels, out_labels, hidden_outs,
    chn, hidden_activations; training_percent = 75,
    epoc = 4, opt = "adam", eta = 0.1
)

app = FluxDash.make_app()

port = haskey(ENV, "PORT") ? parse(Int64, ENV["PORT"]) : 8050

t = @async FluxDash.Dash.run_server(app, "0.0.0.0", port)
sleep(120.0)

max_stars_repo_name: DEICIDIA/fluxdash

max_issues_repo_name: DEICIDIA/fluxdash

max_forks_repo_name: DEICIDIA/fluxdash


text:
function combinevectors(x::Vector{Vector{T}})::Vector{T} where T

    if isempty(x)
        return Vector{T}(undef, 0)
    end

    N = sum(length(x[i]) for i = 1:length(x))

    y = Vector{T}(undef,N)

    st_ind = 0
    fin_ind = 0
    for i = 1:length(x)
        st_ind = fin_ind + 1
        fin_ind = st_ind + length(x[i]) - 1

        y[st_ind:fin_ind] = x[i]
    end

    return y
end


text:
using Documenter, EEGTest

makedocs()

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-material"),
    repo = "github.com/mirestrepo/EEGTest.jl.git",
    julia  = "0.6",
    osname = "linux"
)


text:
struct SingleDecision <: MOI.AbstractScalarFunction
    decision::MOI.VariableIndex
end

mutable struct SingleKnown <: MOI.AbstractScalarFunction
    known::MOI.VariableIndex
end

struct VectorOfDecisions <: MOI.AbstractVectorFunction
    decisions::Vector{MOI.VariableIndex}
end
MOI.output_dimension(f::VectorOfDecisions) = length(f.decisions)

struct VectorOfKnowns <: MOI.AbstractVectorFunction
    knowns::Vector{MOI.VariableIndex}
end
MOI.output_dimension(f::VectorOfKnowns) = length(f.knowns)

# Base overrides #
# ========================== #
Base.copy(f::SingleDecision) = f
Base.copy(f::SingleKnown) = f
Base.copy(f::VectorOfDecisions) = VectorOfDecisions(copy(f.decisions))
Base.copy(f::VectorOfKnowns) = VectorOfKnowns(copy(f.knowns))

Base.iszero(::SingleDecision) = false
Base.isone(::SingleDecision) = false
Base.iszero(::SingleKnown) = false
Base.isone(::SingleKnown) = false

# JuMP overrides #
# ========================== #
function DecisionRef(model::Model, f::SingleDecision)
    return DecisionRef(model, f.decision)
end

function KnownRef(model::Model, f::SingleKnown)
    return KnownRef(model, f.known)
end

SingleDecision(dref::DecisionRef) = SingleDecision(index(dref))
function JuMP.moi_function(dref::DecisionRef)
    return SingleDecision(dref)
end
function JuMP.moi_function_type(::Type{DecisionRef})
    return SingleDecision
end
JuMP.jump_function_type(::Model, ::Type{SingleDecision}) = DecisionRef
function JuMP.jump_function(model::Model, decision::SingleDecision)
    return DecisionRef(model, decision)
end

SingleKnown(kref::KnownRef) = SingleKnown(index(kref))
# function JuMP.moi_function(kref::KnownRef)
#     return AffineDecisionFunction(convert(MOI.ScalarAffineFunction{Float64}, 0.0),
#                                   convert(MOI.ScalarAffineFunction{Float64}, 0.0),
#                                   MOI.ScalarAffineFunction([ScalarAffineTerm(1.0, index(kref))], 0.0))
# end
# function JuMP.moi_function_type(::Type{KnownRef})
#     return AffineDecisionFunction{Float64}
# end
function JuMP.moi_function(kref::KnownRef)
    return SingleKnown(kref)
end
function JuMP.moi_function_type(::Type{KnownRef})
    return SingleKnown
end
JuMP.jump_function_type(::Model, ::Type{SingleKnown}) = KnownRef
function JuMP.jump_function(model::Model, known::SingleKnown)
    return KnownRef(model, known)
end

VectorOfDecisions(dvars::Vector{DecisionRef}) = VectorOfDecisions(index.(dvars))
function JuMP.moi_function(decisions::Vector{<:DecisionRef})
    return VectorOfDecisions(index.(decisions))
end
function JuMP.moi_function_type(::Type{<:Vector{<:DecisionRef}})
    return VectorOfDecisions
end
JuMP.jump_function_type(::Model, ::Type{VectorOfDecisions}) = Vector{DecisionRef}
function JuMP.jump_function(model::Model, decisions::VectorOfDecisions)
    return map(decisions.decisions) do decision
        DecisionRef(model, decision)
    end
end

VectorOfKnowns(kvars::Vector{KnownRef}) = VectorOfKnowns(index.(kvars))
function JuMP.moi_function(knowns::Vector{<:KnownRef})
    return VectorOfKnowns(index.(knowns))
end
function JuMP.moi_function_type(::Type{<:Vector{<:KnownRef}})
    return VectorOfKnowns
end
JuMP.jump_function_type(::Model, ::Type{VectorOfKnowns}) = Vector{KnownRef}
function JuMP.jump_function(model::Model, knowns::VectorOfKnowns)
    return map(knowns.knowns) do known
        KnownRef(model, known)
    end
end

# MOI Function interface #
# ========================== #
MOI.constant(f::SingleDecision, T::DataType) = zero(T)
MOI.constant(f::SingleKnown, T::DataType) = zero(T)

MOIU.eval_variables(varval::Function, f::SingleDecision) = varval(f.decision)
MOIU.eval_variables(::Function, f::SingleKnown) = varval(f.known)

function MOIU.map_indices(index_map::Function, f::SingleDecision)
    return SingleDecision(index_map(f.decision))
end
function MOIU.map_indices(index_map::Function, f::SingleKnown)
    return SingleKnown(index_map(f.known))
end
function MOIU.map_indices(index_map::Function, f::VectorOfDecisions)
    return VectorOfDecisions(index_map.(f.decisions))
end
function MOIU.map_indices(index_map::Function, f::VectorOfKnowns)
    return VectorOfKnowns(index_map.(f.knowns))
end

function Base.getindex(it::MOIU.ScalarFunctionIterator{VectorOfDecisions},
                       i::Integer)
    return SingleDecision(it.f.decisions[i])
end
function Base.getindex(it::MOIU.ScalarFunctionIterator{VectorOfDecisions},
                       I::AbstractVector)
    return VectorOfDecisions(it.f.decisions[I])
end

function Base.getindex(it::MOIU.ScalarFunctionIterator{VectorOfKnowns},
                       i::Integer)
    return SingleKnown(it.f.knowns[i])
end
function Base.getindex(it::MOIU.ScalarFunctionIterator{VectorOfKnowns},
                       I::AbstractVector)
    return VectorOfKnowns(it.f.knowns[I])
end

MOIU.scalar_type(::Type{VectorOfDecisions}) = SingleDecision
MOIU.scalar_type(::Type{VectorOfKnowns}) = SingleKnown

MOIU.canonicalize!(f::Union{SingleDecision, SingleKnown, VectorOfDecisions, VectorOfKnowns}) = f

function MOIU.filter_variables(keep::Function, f::SingleDecision)
    if !keep(f.decision)
        error("Cannot remove decision from a `SingleDecision` function of the",
              " same decision.")
    end
    return f
end

function MOIU.filter_variables(keep::Function, f::SingleKnown)
    if !keep(f.known)
        error("Cannot remove known from a `SingleKnown` function of the",
              " same known decision.")
    end
    return f
end

function MOIU.filter_variables(keep::Function, f::VectorOfDecisions)
    return VectorOfDecisions(MOIU._filter_variables(keep, f.decisions))
end

function MOIU.filter_variables(keep::Function, f::VectorOfKnowns)
    return VectorOfKnowns(MOIU._filter_variables(keep, f.knowns))
end

function MOIU.vectorize(funcs::AbstractVector{SingleDecision})
    decisions = MOI.VariableIndex[f.decision for f in funcs]
    return VectorOfDecisions(vars)
end

function MOIU.vectorize(funcs::AbstractVector{SingleKnown})
    knowns = MOI.VariableIndex[f.known for f in funcs]
    return VectorOfKnowns(knowns)
end

function MOIU.scalarize(f::VectorOfDecisions, ignore_constants::Bool = false)
    SingleDecision.(f.decisions)
end

function MOIU.scalarize(f::VectorOfKnowns, ignore_constants::Bool = false)
    SingleKnown.(f.knowns)
end


text:
# Posterior samples and mean estimation from a given initial point
using EstimateSDE
using EstimateSDE: Parameters
using Serialization
using Random
using Parameters

## Settings
settings = (
  proc=testprocess(:B),
  N=100,
  Tav=1.0,
  Tkld=1.0,
)

@unpack proc,N,Tav,Tkld = settings
rng = MersenneTwister(12202)
data = rand(rng,proc,N,Tav)
problem = ProblemProb(data;p=EstimateSDE.order(proc),engine=ElimQ)

## Maximum likelihood estimate
mltor = MlEstimator()
init_tors = default_init_tors(rng,[0.5,1.0,2.0]*Tkld)
proc_ml = estimate(mltor,problem.like;init_tors=init_tors,throw_errors=false)
if proc_ml isa WhiteNoise
  error("No error; correct behavior: ML estimate is white noise; no posterior samples are taken.")
end

## Posterior samples
tor = PmeanEstimator(;rng=MersenneTwister(54858979),N=200,Nwarmup=200,nchains=4)
procm,procs,R = @time estimate(tor,problem,proc_ml)

## Save settings, data and results
d = EstimateSDE.resultsdir("a-samples")
fn = joinpath(d,"samples.jld")
res = Dict{Symbol,Any}()
@pack! res = settings,proc_ml,problem,procm,procs,R
serialize(fn,res)
@info "Settings and results saved to $fn"


text:
#=
Filename: aiyagari_household.jl
Author: Victoria Gregory
Date: 8/29/2016

This file defines the Household type (and its constructor)
for setting up an Aiyagari household problem.
=#

using QuantEcon

"""
Stores all the parameters that define the household's
problem.

##### Fields

- `r::Real` : interest rate
- `w::Real` : wage
- `sigma::Real` : risk aversion
- `beta::AbstractFloat` : discount factor
- `z_chain::MarkovChain` : MarkovChain for income
- `a_min::Real` : minimum on asset grid
- `a_max::Real` : maximum on asset grid
- `a_size::Integer` : number of points on asset grid
- `z_size::Integer` : number of points on income grid
- `n::Integer` : number of points in state space: (a, z)
- `s_vals::Array{TF} where TF<:AbstractFloat` : stores all the possible (a, z) combinations
- `s_i_vals::Array{TI} where TI<:Integer` : stores indices of all the possible (a, z) combinations
- `R::Array{TF} where TF<:AbstractFloat` : reward array
- `Q::Array{TF} where TF<:AbstractFloat` : transition probability array
- `u::Function` : utility function
"""
mutable struct Household{TR<:Real, TF<:AbstractFloat, TI<:Integer}
    r::TR
    w::TR
    sigma::TR
    beta::TF
    z_chain::MarkovChain{TF, Array{TF, 2}, Array{TF, 1}}
    a_min::TR
    a_max::TR
    a_size::TI
    a_vals::AbstractVector{TF}
    z_size::TI
    n::TI
    s_vals::Array{TF}
    s_i_vals::Array{TI}
    R::Array{TR}
    Q::Array{TR}
    u::Function
end

"""
Constructor for `Household`

##### Arguments
- `r::Real(0.01)` : interest rate
- `w::Real(1.0)` : wage
- `beta::AbstractFloat(0.96)` : discount factor
- `z_chain::MarkovChain` : MarkovChain for income
- `a_min::Real(1e-10)` : minimum on asset grid
- `a_max::Real(18.0)` : maximum on asset grid
- `a_size::TI(200)` : number of points on asset grid

"""
function Household{TF<:AbstractFloat}(;
                    r::Real=0.01,
                    w::Real=1.0,
                    sigma::Real=1.0,
                    beta::TF=0.96,
                    z_chain::MarkovChain{TF,Array{TF,2},Array{TF,1}}
                        =MarkovChain([0.9 0.1; 0.1 0.9], [0.1; 1.0]),
                    a_min::Real=1e-10,
                    a_max::Real=18.0,
                    a_size::Integer=200)

    # set up grids
    a_vals = linspace(a_min, a_max, a_size)
    z_size = length(z_chain.state_values)
    n = a_size*z_size
    s_vals = gridmake(a_vals, z_chain.state_values)
    s_i_vals = gridmake(1:a_size, 1:z_size)

    # set up Q
    Q = zeros(TF, n, a_size, n)
    for next_s_i in 1:n
        for a_i in 1:a_size
            for s_i in 1:n
                z_i = s_i_vals[s_i, 2]
                next_z_i = s_i_vals[next_s_i, 2]
                next_a_i = s_i_vals[next_s_i, 1]
                if next_a_i == a_i
                    Q[s_i, a_i, next_s_i] = z_chain.p[z_i, next_z_i]
                end
            end
        end
    end

    if sigma == 1   # log utility
        u = x -> log(x)
    else
        u = x -> (x^(1-sigma)-1)/(1-sigma)
    end

    # placeholder for R
    R = fill(-Inf, n, a_size)
    h = Household(r, w, sigma, beta, z_chain, a_min, a_max, a_size,
                  a_vals, z_size, n, s_vals, s_i_vals, R, Q, u)

    setup_R!(h, r, w)

    return h

end

"""
Update the reward array of a Household object, given
a new interest rate and wage.

##### Arguments
- `h::Household` : instance of Household type
- `r::Real(0.01)` : interest rate
- `w::Real(1.0)` : wage
"""
function setup_R!(h::Household, r::Real, w::Real)

    # set up R
    R = h.R
    for new_a_i in 1:h.a_size
        a_new = h.a_vals[new_a_i]
        for s_i in 1:h.n
            a = h.s_vals[s_i, 1]
            z = h.s_vals[s_i, 2]
            c = w * z + (1 + r) * a - a_new
            if c > 0
                R[s_i, new_a_i] = h.u(c)
            end
        end
    end

    h.r = r
    h.w = w
    h.R = R
    h

end


text:
using Gadfly
using DataFrames
using StatsBase

function bc_upd(ϵ::Float64, α::Float64, τ::Float64, n_agents::Int, n_steps::Int)
    bc_ar = Array{Float64,2}(undef, n_agents, n_steps + 1)
    bc_ar[:, 1] = rand(n_agents)
    @views for j in 1:n_steps, i in 1:n_agents
        @inbounds bc_ar[i, j + 1] = (1 - α)*mean(bc_ar[:, j][findall(abs.(bc_ar[:, j] .- bc_ar[:, j][i]) .< ϵ)]) + α*τ
    end
    return bc_ar
end

const ϵ = 0.1
const α = 0.2
const τ = 0.7

res = bc_upd(ϵ, α, τ, 50, 30)

df₀ = res |> rotr90 |> DataFrame
rename!(df₀, [Symbol("$i") for i in 1:size(df₀, 2)])
df = df₀ |> stack
df[!, :steps] = repeat(0:size(df₀, 1) - 1, outer=size(df₀, 2))

plot(df, x=:steps, y=:value, color=:variable, Geom.point, Geom.line,
    Coord.cartesian(xmax=size(df₀, 1)),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    Guide.title("ϵ = $ϵ / α = $α / τ = $τ"),
    Theme(key_position=:none, point_size=1.5pt))

max_stars_repo_name: IgorDouven/Misinformation

max_issues_repo_name: IgorDouven/Misinformation

max_forks_repo_name: IgorDouven/Misinformation


text:
function initial_head(x)
    2 * √x
end

"""
    transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, Β)

Non-steady flow in an unconfined rectangular aquifer, with Dirichlet h(0, t) = 0
on the left edge, and a Neumann Boundary Condition (dh/dx = 0) on the right.
"""
function transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, Β)
    initial_head(x) / 1.0 + (Β * conductivity * initial_head(aquifer_length) * time) / (specific_yield * aquifer_length * aquifer_length)
end

"""
    drawdown_theis(distance, time, discharge, transmissivity, storativity)

Non-steady flow in a confined aquifer, using the well function of Theis.
"""
function drawdown_theis(distance, time, discharge, transmissivity, storativity)
    u = (storativity * distance^2) / (4 * transmissivity * time)
    return discharge / (4 * π  * transmissivity) * expint(u)
end

function homogenous_aquifer(nrow, ncol)
    shape = (nrow, ncol)
    # Domain, geometry
    domain = ones(Bool, shape)
    Δx = fill(10.0, ncol)
    Δy = fill(10.0, nrow)
    indices, reverse_indices = Wflow.active_indices(domain, false)
    connectivity = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
    ncell = connectivity.ncell

    conf_aqf = Wflow.ConfinedAquifer(
        [0.0, 7.5, 20.0],  # head
        fill(10.0, ncell),  # k
        fill(10.0, ncell),  # top
        fill(0.0, ncell),  # bottom
        fill(100.0, ncell),  # area
        fill(0.1, ncell), # specific storage
        fill(1.0, ncell),  # storativity
        fill(0.0, connectivity.nconnection),  # conductance
    )
    unconf_aqf = Wflow.UnconfinedAquifer(
        [0.0, 7.5, 20.0],  # head
        fill(10.0, ncell),   # k
        fill(10.0, ncell),  # top
        fill(0.0, ncell),  # bottom
        fill(100.0, ncell),  # area
        fill(0.15, ncell),  # specific yield
        fill(0.0, connectivity.nconnection),  # conductance
    )
    return (connectivity, conf_aqf, unconf_aqf)
end

@testset "groundwater" begin
    ncol = 2
    nrow = 3
    shape = (ncol, nrow)
    Δx = [10.0, 20.0]
    Δy = [5.0, 15.0, 25.0]
    collect_connections(con, cell_id) = [con.rowval[nzi] for nzi in Wflow.connections(con, cell_id)]

    @testset "unit: connectivity" begin
        @testset "connection_geometry: y" begin
            I = CartesianIndex(1, 1)
            J = CartesianIndex(2, 1)
            @test Wflow.connection_geometry(I, J, Δx, Δy) == (2.5, 7.5, 10.0)
            @test_throws Exception Wflow.connection_geometry(I, I, Δx, Δy)
        end

        @testset "connection_geometry: x" begin
            I = CartesianIndex(1, 1)
            J = CartesianIndex(1, 2)
            @test Wflow.connection_geometry(I, J, Δx, Δy) == (5.0, 10.0, 5.0)
            @test_throws Exception Wflow.connection_geometry(I, I, Δx, Δy)
        end

        @testset "Connectivity 1D(x)" begin
            # +---+---+
            # | 1 | 2 |
            # +---+---+
            domain = ones(Bool, (1, 2))
            indices, reverse_indices = Wflow.active_indices(domain, false)
            conn = Wflow.Connectivity(indices, reverse_indices, Δx, [5.0])
            @test conn.ncell == 2
            @test conn.nconnection == 2
            @test conn.length1 == [5.0, 10.0]
            @test conn.length2 == [10.0, 5.0]
            @test conn.width == [5.0, 5.0]
            @test conn.colptr == [1, 2, 3]
            @test conn.rowval == [2, 1]
            @test collect_connections(conn, 1) == [2]
            @test collect_connections(conn, 2) == [1]
        end

        @testset "Connectivity 1D(y)" begin
            # +---+
            # | 1 |
            # +---+
            # | 2 |
            # +---+
            # | 3 |
            # +---+
            domain = ones(Bool, (3, 1))
            indices, reverse_indices = Wflow.active_indices(domain, false)
            conn = Wflow.Connectivity(indices, reverse_indices, [10.0], Δy)
            @test conn.ncell == 3
            @test conn.nconnection == 4
            @test conn.length1 == [2.5, 7.5, 7.5, 12.5]
            @test conn.length2 == [7.5, 2.5, 12.5, 7.5]
            @test conn.width == [10.0, 10.0, 10.0, 10.0]
            @test conn.colptr == [1, 2, 4, 5]
            @test conn.rowval == [2, 1, 3, 2]
            @test collect_connections(conn, 1) == [2]
            @test collect_connections(conn, 2) == [1, 3]
            @test collect_connections(conn, 3) == [2]
        end

        @testset "Connectivity 2D" begin
            # +---+---+
            # | 1 | 4 |
            # +---+---+
            # | 2 | 5 |
            # +---+---+
            # | 3 | 6 |
            # +---+---+
            domain = ones(Bool, (nrow, ncol))
            indices, reverse_indices = Wflow.active_indices(domain, false)
            conn = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
            @test conn.ncell == 6
            @test conn.nconnection == 14
            @test conn.colptr == [1, 3, 6, 8, 10, 13, 15]
            @test conn.rowval == [2, 4, 1, 3, 5, 2, 6, 1, 5, 2, 4, 6, 3, 5]
            @test collect_connections(conn, 1) == [2, 4]
            @test collect_connections(conn, 2) == [1, 3, 5]
            @test collect_connections(conn, 3) == [2, 6]
            @test collect_connections(conn, 4) == [1, 5]
            @test collect_connections(conn, 5) == [2, 4, 6]
            @test collect_connections(conn, 6) == [3, 5]
        end

        @testset "Connectivity 2D - partially inactive" begin
            # +---+---+
            # | 1 | 3 |
            # +---+---+
            # | X | 4 |
            # +---+---+
            # | 2 | 5 |
            # +---+---+
            domain = ones(Bool, (nrow, ncol))
            domain[2, 1] = false
            indices, reverse_indices = Wflow.active_indices(domain, false)
            conn = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
            @test conn.ncell == 5
            @test conn.nconnection == 8
            @test conn.colptr == [1, 2, 3, 5, 7, 9]
            @test conn.rowval == [3, 5, 1, 4, 3, 5, 2, 4]
            @test collect_connections(conn, 1) == [3]
            @test collect_connections(conn, 2) == [5]
            @test collect_connections(conn, 3) == [1, 4]
            @test collect_connections(conn, 4) == [3, 5]
            @test collect_connections(conn, 5) == [2, 4]
        end
    end

    @testset "unit: aquifer, boundary conditions" begin
        @testset "harmonicmean_conductance" begin
            # harmonicmean_conductance(k1, k2, H1, H2, l1, l2, width)
            @test Wflow.harmonicmean_conductance(10.0, 10.0, 5.0, 5.0, 0.5, 0.5, 1.0) == 50.0
            @test Wflow.harmonicmean_conductance(10.0, 10.0, 0.0, 5.0, 0.5, 0.5, 1.0) == 0.0
            @test Wflow.harmonicmean_conductance(10.0, 10.0, 5.0, 0.0, 0.5, 0.5, 1.0) == 0.0
            # kD of 10 and 20 -> harmonicmean = 1/(1/10 + 1/20)
            @test Wflow.harmonicmean_conductance(10.0, 10.0, 1.0, 2.0, 1.0, 1.0, 1.0) ≈ (6.0 + 2.0 / 3.0)
        end

        nrow = 1
        ncol = 3
        connectivity, conf_aqf, unconf_aqf = homogenous_aquifer(nrow, ncol)
        Wflow.initialize_conductance!(conf_aqf, connectivity)
        Wflow.initialize_conductance!(unconf_aqf, connectivity)
        ncell = connectivity.ncell

        @testset "saturated_thickness-confined" begin
            @test (
                Wflow.saturated_thickness(conf_aqf, 1) ==
                Wflow.saturated_thickness(conf_aqf, 2) ==
                Wflow.saturated_thickness(conf_aqf, 3) ==
                10.0
                )
        end

        @testset "saturated_thickness-unconfined" begin
            @test Wflow.saturated_thickness(unconf_aqf, 1) == 0.0
            @test Wflow.saturated_thickness(unconf_aqf, 2) == 7.5
            @test Wflow.saturated_thickness(unconf_aqf, 3) == 10.0
        end

        @testset "horizontal_conductance" begin
            @test (
                Wflow.horizontal_conductance(1, 2, 1, conf_aqf, connectivity) ==
                Wflow.harmonicmean_conductance(10.0, 10.0, 10.0, 10.0, 5.0, 5.0, 10.0)
            )
        end

        @testset "conductance" begin
            @test Wflow.conductance(conf_aqf, 2, 3, 3) == 100.0
            @test Wflow.conductance(unconf_aqf, 2, 3, 3) == 100.0  # upstream sat. thickness
            @test Wflow.conductance(unconf_aqf, 1, 2, 1) == 75.0  # upstream sat. thickness
        end

        @testset "minimum_head-confined" begin
            original_head = copy(conf_aqf.head)
            conf_aqf.head[1] = -10.0
            @test Wflow.minimum_head(conf_aqf)[1] == -10.0
            conf_aqf.head .= original_head
        end

        @testset "minimum_head-unconfined" begin
            original_head = copy(unconf_aqf.head)
            unconf_aqf.head[1] = -10.0
            @test Wflow.minimum_head(conf_aqf)[1] == 0.0
            unconf_aqf.head .= original_head
        end

        @testset "stable_timestep" begin
            @test Wflow.stable_timestep(conf_aqf) == 0.25
        end

        # Parametrization in setup is as follows:
        # [0.0, 7.5, 20.0],  # head
        # fill(10.0, ncell),  # k
        # fill(10.0, ncell),  # top
        # fill(0.0, ncell),  # bottom

        @testset "flux-confined" begin
            Q = zeros(3)
            Wflow.flux!(Q, conf_aqf, connectivity)
            # kD = 10 * 10 = 100
            # dH = 7.5, 12.5
            @test Q == [750.0, 500.0, -1250.0]  # TODO
        end

        @testset "flux-unconfined" begin
            Q = zeros(3)
            Wflow.flux!(Q, unconf_aqf, connectivity)
            # KD is based on upstream saturated thickness, i.e. 7.5 m and 20.0 m (which is capped to 10.0)
            @test Q == [562.5, 687.5, -1250.0]
        end

        @testset "river" begin
            river = Wflow.River(
                [2.0, 2.0], [100.0, 100.0], [200.0, 200.0], [1.0, 1.0], [0.0, 0.0], [1, 3]
            )
            Q = zeros(3)
            Wflow.flux!(Q, river, conf_aqf)
            # infiltration, below bottom, flux is (stage - bottom) * inf_cond
            @test Q[1] == (2.0 - 1.0) * 100.0
            # drainage, flux is () * exf_cond
            @test Q[3] == (2.0 - 20.0) * 200.0
        end

        @testset "drainage" begin
            drainage = Wflow.Drainage([2.0, 2.0], [100.0, 100.0], [0.0, 0.0], [1, 2])
            Q = zeros(3)
            Wflow.flux!(Q, drainage, conf_aqf)
            @test Q[1] == 0.0
            @test Q[2] == 100.0 * (2.0 - 7.5)
        end

        @testset "headboundary" begin
            headboundary = Wflow.HeadBoundary([2.0, 2.0], [100.0, 100.0], [0.0, 0.0], [1, 2])
            Q = zeros(3)
            Wflow.flux!(Q, headboundary, conf_aqf)
            @test Q[1] == 100.0 * (2.0 - 0.0)
            @test Q[2] == 100.0 * (2.0 - 7.5)
        end

        @testset "recharge" begin
            recharge = Wflow.Recharge([1.0e-3, 1.0e-3, 1.0e-3], [0.0, 0.0, 0.0], [1, 2, 3])
            Q = zeros(3)
            Wflow.flux!(Q, recharge, conf_aqf)
            @test all(Q .== 1.0e-3 * 100.0)
        end

        @testset "well" begin
            well = Wflow.Well([-1000.0], [0.0], [1])
            Q = zeros(3)
            Wflow.flux!(Q, well, conf_aqf)
            @test Q[1] == -1000.0
        end
    end

    @testset "integration: steady 1D" begin
        connectivity, aquifer, _ = homogenous_aquifer(3, 1)
        constanthead = Wflow.ConstantHead([2.0, 4.0], [1, 3])
        gwf = Wflow.GroundwaterFlow(
           aquifer,
           connectivity,
           constanthead,
           Wflow.AquiferBoundaryCondition[],
        )
        # Set constant head (dirichlet) boundaries
        gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head

        Q = zeros(3)
        Δt = 0.25 # days
        for _ in 1:50
            Wflow.update(gwf, Q, Δt)
        end

        @test gwf.aquifer.head ≈ [2.0, 3.0, 4.0]
    end

    @testset "integration: unconfined transient 1D" begin
        nrow = 1
        ncol = 9
        shape = (nrow, ncol)
        conductivity = 200.0
        top = 150.0
        bottom = 0.0
        specific_yield = 0.15
        cellsize = 500.0
        Β = 1.12
        aquifer_length = cellsize * ncol

        # Domain, geometry
        domain = ones(Bool, shape)
        Δx = fill(cellsize, ncol)
        Δy = fill(cellsize, nrow)
        indices, reverse_indices = Wflow.active_indices(domain, false)
        connectivity = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
        ncell = connectivity.ncell
        xc = collect(range(0.0, stop=aquifer_length - cellsize, step=cellsize))
        aquifer = Wflow.UnconfinedAquifer(
            initial_head.(xc),
            fill(conductivity, ncell),
            fill(top, ncell),
            fill(bottom, ncell),
            fill(cellsize * cellsize, ncell),
            fill(specific_yield, ncell),
            fill(0.0, connectivity.nconnection),
        )
        # constant head on left boundary, 0 at 0
        constanthead = Wflow.ConstantHead([0.0], [1])
        gwf = Wflow.GroundwaterFlow(
           aquifer,
           connectivity,
           constanthead,
           Wflow.AquiferBoundaryCondition[],
        )

        Δt = Wflow.stable_timestep(gwf.aquifer)
        Q = zeros(ncell)
        time = 20.0
        nstep = Int(ceil(time / Δt))
        time = nstep * Δt

        for i in 1:nstep
            Wflow.update(gwf, Q, Δt)
            # Gradient dh/dx is positive, all flow to the left
            @test all(diff(gwf.aquifer.head) .> 0.0)
        end

        ϕ_analytical = [transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, Β) for x in xc]
        difference = gwf.aquifer.head .- ϕ_analytical
        # @test all(difference .< ?)  #TODO
    end

    @testset "integration: confined transient radial 2D" begin
        halfnrow = 20
        wellrow = halfnrow + 1
        nrow = halfnrow * 2 + 1
        ncol = nrow
        shape = (nrow, ncol)
        conductivity = 5.0
        top = 10.0
        bottom = 0.0
        transmissivity = (top - bottom) * conductivity
        cellsize = 10.0
        startinghead = top
        specific_storage = 0.015
        storativity = 0.15
        aquifer_length = cellsize * ncol
        discharge = -50.0

        # Domain, geometry
        domain = ones(Bool, shape)
        Δx = fill(cellsize, ncol)
        Δy = fill(cellsize, nrow)
        indices, reverse_indices = Wflow.active_indices(domain, false)
        connectivity = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
        ncell = connectivity.ncell
        aquifer = Wflow.ConfinedAquifer(
            fill(startinghead, ncell),
            fill(conductivity, ncell),
            fill(top, ncell),
            fill(bottom, ncell),
            fill(cellsize * cellsize, ncell),
            fill(specific_storage, ncell),
            fill(storativity, ncell),
            fill(0.0, connectivity.nconnection), # conductance, to be set
        )

        cell_index = reshape(collect(range(1, ncell, step=1)), shape)
        indices = vcat(cell_index[1, :], cell_index[end, :])# , cell_index[:, 1], cell_index[:, end],)
        constanthead = Wflow.ConstantHead(
            fill(10.0, size(indices)),
            indices,
        )
        # Place a well in the middle of the domain
        well = Wflow.Well([discharge], [0.0], [reverse_indices[wellrow, wellrow]])
        gwf = Wflow.GroundwaterFlow(
            aquifer,
            connectivity,
            constanthead,
            [well],
        )

        Δt = Wflow.stable_timestep(gwf.aquifer)
        Q = zeros(ncell)
        time = 20.0
        nstep = Int(ceil(time / Δt))
        time = nstep * Δt

        for i in 1:nstep
            Wflow.update(gwf, Q, Δt)
        end

        # test for symmetry on x and y axes
        ϕ = reshape(gwf.aquifer.head, shape)
        @test ϕ[1:halfnrow, :] ≈ ϕ[end:-1:halfnrow + 2, :]
        @test ϕ[:, 1:halfnrow] ≈ ϕ[:, end:-1:halfnrow + 2]

        # compare with analytical solution
        start = -0.5 * aquifer_length + 0.5 * cellsize
        stop =  0.5 * aquifer_length - 0.5 * cellsize
        X = collect(range(start, stop=stop, step=cellsize))
        ϕ_analytical = [drawdown_theis(x, time, discharge, transmissivity, storativity) for x in X] .+ 10.0
        # compare left-side, since it's symmetric anyway. Skip the well cell, and its first neighbor
        difference = ϕ[1:halfnrow - 1, halfnrow] - ϕ_analytical[1:halfnrow - 1]
        @test all(difference .< 0.02)
    end
end


text:
using LoopVectorization, StrideArraysCore, Test

function issue238_noavx!(output, matrix, input)
  for i in axes(output, 2), v in axes(output, 1)
    res = zero(eltype(output))
    for ii in axes(matrix, 2)
      res += matrix[i, ii] * input[v, ii]
    end
    output[v, i] = res
  end
  return nothing
end
function issue238!(output, matrix, input)
  @turbo inline=true for i in axes(output, 2), v in axes(output, 1)
    res = zero(eltype(output))
    for ii in axes(matrix, 2)
      res += matrix[i, ii] * input[v, ii]
    end
    output[v, i] += res
  end
  return nothing
end
function issue238_v2!(output, matrix, input)
  @turbo inline=true for i in axes(output, 2), v in axes(output, 1)
    res = output[v, i]
    for ii in axes(matrix, 2)
      res += matrix[i, ii] * input[v, ii]
    end
    output[v, i] = res
  end
  return nothing
end
const MAXTESTSIZE = 8 # N^3 loop iterations, meaning N^3 different functions compiled
function n2testloop(output1,output2,output3,output_nonstatic0,output_nonstatic1)
  n1, n3 = size(output1)
  for n2 ∈ 1:MAXTESTSIZE;
    # @show n1, n2, n3
    input  = randn(n1, n2)
    matrix = randn(n3, n2)
    smatrix = StrideArray(matrix, (StaticInt(n1), StaticInt(n2)))
    issue238_noavx!(output_nonstatic0, matrix, input)
    fill!(output_nonstatic1, 0); issue238!(output_nonstatic1, matrix, input)
    fill!(output1,0); issue238!(output1, matrix, input)
    fill!(output2,0); issue238!(output2, smatrix, input)
    fill!(output3,0); issue238_v2!(output3, smatrix, input)
    if !(output1 ≈ output2 ≈ output3 ≈ output_nonstatic0 ≈ output_nonstatic1)
      @show n1, n2, n3
    end
    @test output_nonstatic0 ≈ output1 # static test
    @test output_nonstatic0 ≈ output2 # static test
    @test output_nonstatic0 ≈ output3 # static test
    @test output_nonstatic0 ≈ output_nonstatic1 # nonstatic test
  end
end

@testset "Statically Sized Arrays" begin
  @show @__LINE__
  for n1 ∈ 1:MAXTESTSIZE, n3 ∈ 1:MAXTESTSIZE
    output1 = StrideArray(undef, StaticInt(n1), StaticInt(n3))
    output2 = StrideArray(undef, StaticInt(n1), StaticInt(n3))
    output3 = StrideArray(undef, StaticInt(n1), StaticInt(n3))
    output_nonstatic0 = StrideArray(undef, n1, n3)
    output_nonstatic1 = StrideArray(undef, n1, n3)
    n2testloop(output1,output2,output3,output_nonstatic0,output_nonstatic1)
  end
end


text:
using CUSPARSE
using CUDArt
using Base.Test

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

##############
# test_axpyi #
##############
function test_axpyi!(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    alpha = rand(elty)
    d_y = CUSPARSE.axpyi!(alpha,d_x,d_y,'O')
    #compare
    h_y = to_host(d_y)
    if VERSION >= v"0.5.0-dev+742"
        y[x.nzind] += alpha * x.nzval
    else
        y[x.rowval] += alpha * x.nzval
    end
    @test h_y ≈ y
end

function test_axpyi(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    alpha = rand(elty)
    d_z = CUSPARSE.axpyi(alpha,d_x,d_y,'O')
    #compare
    h_z = to_host(d_z)
    z = copy(y)
    if VERSION >= v"0.5.0-dev+742"
        z[x.nzind] += alpha * x.nzval
    else
        z[x.rowval] += alpha * x.nzval
    end
    @test h_z ≈ z
    d_z = CUSPARSE.axpyi(d_x,d_y,'O')
    #compare
    h_z = to_host(d_z)
    z = copy(y)
    if VERSION >= v"0.5.0-dev+742"
        z[x.nzind] += x.nzval
    else
        z[x.rowval] += x.nzval
    end
    @test h_z ≈ z
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_axpyi!(elty)
    println("axpyi! took ", toq(), " for ", elty)
    tic()
    test_axpyi(elty)
    println("axpyi took ", toq(), " for ", elty)
end


text:
module SemiCtsDistributions

# package code goes here
export rsemicts, hist, summary, pred.perf

include("semicts.jl")

end # module

max_stars_repo_name: psaiku/semictsdistributions

max_issues_repo_name: psaiku/semictsdistributions

max_forks_repo_name: psaiku/semictsdistributions


text:
using Test

include("rna-transcription.jl")

@testset "basic transformations" begin
    @testset "empty rna sequence" begin
        @test to_rna("") == ""
    end

    @testset "rna complement of cytosine is guanine" begin
        @test to_rna("C") == "G"
    end

    @testset "rna complement of guanine is cytosine" begin
        @test to_rna("G") == "C"
    end

    @testset "rna complement of thymine is adenine" begin
        @test to_rna("T") == "A"
    end

    @testset "rna complement of adenine is uracil" begin
        @test to_rna("A") == "U"
    end
end

@testset "rna complement" begin
    @test to_rna("ACGTGGTCTTAA") == "UGCACCAGAAUU"
end

@testset "error handling" begin
    @testset "dna correctly handles invalid input" begin
        @test_throws ErrorException to_rna("U")
    end

    @testset "dna correctly handles completely invalid input" begin
        @test_throws ErrorException to_rna("XXX")
    end

    @testset "dna correctly handles partially invalid input" begin
        @test_throws ErrorException to_rna("ACGTXXXCTTAA")
    end
end


text:
println("Frank sez hi")
using BinDeps

@BinDeps.setup

mpi = library_dependency("juliampi", aliases=["libjuliampi", "libjuliampi.so", "libjuliampi.dylib"], os = :Unix)
builddir = joinpath(BinDeps.depsdir(mpi), "build")
prefix = joinpath(BinDeps.depsdir(mpi), "usr")
src = joinpath(BinDeps.depsdir(mpi), "src")

if isdir(builddir)
    rm(builddir, recursive = true)
end
if isdir(prefix)
    rm(prefix, recursive = true)
end

# Add cmake flags to specify the MPI library location
configure = `cmake`
flags = ["MPI_C_COMPILER",
         "MPI_C_COMPILE_FLAGS",
         "MPI_C_INCLUDE_PATH",
         "MPI_C_LINK_FLAGS",
         "MPI_C_LIBRARIES",
         "MPI_Fortran_COMPILER",
         "MPI_Fortran_COMPILE_FLAGS",
         "MPI_Fortran_INCLUDE_PATH",
         "MPI_Fortran_LINK_FLAGS",
         "MPI_Fortran_LIBRARIES",
         "MPI_INCLUDE_PATH",
         "MPI_LIBRARIES"]
for flag in flags
    try
        val = ENV["JULIA_$flag"]
        global configure = `$configure -D$flag=$val`
    catch
        nothing
    end
end
configure = `$configure -DCMAKE_INSTALL_PREFIX=$src -DCMAKE_LIB_INSTALL_PREFIX=$prefix/lib ..`

provides(SimpleBuild,
    (@build_steps begin
        CreateDirectory(builddir)
        (@build_steps begin
            ChangeDirectory(builddir)
            (@build_steps begin
                configure
                `make`
                `make install`
            end)
        end)
    end), [mpi], os = :Unix)

@BinDeps.install Dict(:mpi => :mpi)


text:
# Autogenerated wrapper script for Qt5Declarative_jll for powerpc64le-linux-gnu-cxx03
export libqt5qml, libqt5qmlmodels, libqt5qmlworkerscript, libqt5quick, libqt5quickparticles, libqt5quickshapes, libqt5quicktest, libqt5quickwidgets

using Qt5Base_jll
JLLWrappers.@generate_wrapper_header("Qt5Declarative")
JLLWrappers.@declare_library_product(libqt5qml, "libQt5Qml.so.5")
JLLWrappers.@declare_library_product(libqt5qmlmodels, "libQt5QmlModels.so.5")
JLLWrappers.@declare_library_product(libqt5qmlworkerscript, "libQt5QmlWorkerScript.so.5")
JLLWrappers.@declare_library_product(libqt5quick, "libQt5Quick.so.5")
JLLWrappers.@declare_library_product(libqt5quickparticles, "libQt5QuickParticles.so.5")
JLLWrappers.@declare_library_product(libqt5quickshapes, "libQt5QuickShapes.so.5")
JLLWrappers.@declare_library_product(libqt5quicktest, "libQt5QuickTest.so.5")
JLLWrappers.@declare_library_product(libqt5quickwidgets, "libQt5QuickWidgets.so.5")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll)
    JLLWrappers.@init_library_product(
        libqt5qml,
        "lib/libQt5Qml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlmodels,
        "lib/libQt5QmlModels.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlworkerscript,
        "lib/libQt5QmlWorkerScript.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick,
        "lib/libQt5Quick.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickparticles,
        "lib/libQt5QuickParticles.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickshapes,
        "lib/libQt5QuickShapes.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktest,
        "lib/libQt5QuickTest.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickwidgets,
        "lib/libQt5QuickWidgets.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MaybeEither
using Test

@testset "MaybeEither" begin
    @testset "Either" begin
        @test left(3) isa Left{Int, Any}
        @test left(3) isa Either{Int, Any}

        @test right("Sdfds", Vector) isa Either{Vector,String}

        @test right(right(0.3, String)) == Right{String,Float64}(0.3)
        @test left(left(0.3, Char)) == Left{Float64,Char}(0.3)

        @test value(left(16.0im)) == 16.0im
        @test value(right(16.0im)) == 16.0im

        @test isleft(left(0.3)) == true
        @test isleft(right(0.3)) == false

        @test isright(left(0.3)) == false
        @test isright(right(1)) == true
        @test isright(Vector) == false

        l = Left{Float32, String}(left(1, Array))
        @test l == Left{String}(1.0f0)
        @test lefttype(l) == Float32
        @test righttype(l) == String

        r = Right{Int, Float32}(right(1, String))
        @test r == Right{Int}(1.0f0)
        @test lefttype(r) == Int
        @test righttype(r) == Float32

        @test Left{ComplexF64}(Left{AbstractFloat,String}(0.9)) isa
                Left{AbstractFloat, ComplexF64}
        @test Right{String}(Right{Int,Number}(0)) isa Right{String,Number}

        @testset "joins" begin
            @test join(left(right(0.4, 'a'))) == right(0.4, Char)

            @test join(left(right(right(4, left(left(3)))), 'a')) ==
                    Right{Left{Left{Int64,Any},Any},Int}(4)

        end
    end

    @testset "Maybe" begin
        j = Just{Integer}(5)
        @test eltype(j) == Integer
        @test j != Just(5)

        m = Just(Just(Just(None("no"))))
        @test join(m) == None(String)
        @test eltype(m) == String
        @test value(m) == nothing

        j = Just("Hello")
        @test map(length, j, Vector) == map(length, j) == Just(5)
        n = None(String)
        @test map(length, n) == n
        @test map(length, n, Int) == n

        @test join(map(_ -> nothing, Just(4))) == None(Int)
    end
end


text:
pint(s) = parse(Int, s)

max = 0
for a = 1:100, b = 1:100
  global max
  n = BigInt(a) ^ BigInt(b)
  digits = split(string(n), "")
  s = sum(pint, digits)
  if s > max
    max = s
  end
end
println(max)

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff