text:
module World

using Random

include("coordinates.jl")
#using .Coordinatesm

include("resources.jl")
include("organism.jl")

include("display.jl")

struct Parameters
    time::Float64
    energy::Int
    cells::Int
    org_max_energie::Int
    res_max_count::Int
    org_max_new::Int
    split_boundaries::Tuple{Float64,Float64}
    function Parameters(args...)
        max_res=3000
        max_energy=max_res*50
        new(
            floatmin(),     #start time
            max_energy,           #world energy
            100000,          #number of world cells on the sphere
            100,            #maximun possible energy of an organism
            max_res,           #maximum number of resources
            5,              #maximum number of new organisms per step
            (0.0,10.0)      #arbitrary boundaries for split decision
        )
    end
end

struct DerivedParameters
    possible_coordinates::PossibleCoordinates
    function DerivedParameters(wp::Vector{Parameters})
        new(calc_possible_coordinates(wp[1].cells))
    end
end

struct StateHelpers
    res_pol_positions2resource::Dict{Coordinates,Array{Resource,1}}
    res_index_positions2resource::Dict{Coordinates2dIndex,Array{Resource,1}}
    resource2expunge::Array{Resource,1}
    function StateHelpers()
        res_pol_positions2resource=Dict{Coordinates,Array{Resource,1}}()
        res_index_positions2resource=Dict{Coordinates2dIndex,Array{Resource,1}}()
        resource2expunge=Array{Resource,1}()
        new(res_pol_positions2resource,
            res_index_positions2resource,
            resource2expunge,
            )
    end
end

@enum StateFloatIndices currentTime=1 maxLifeTime randomizeDelta
@enum StateIntIndices freeEnergy=1 step maxsteps cells res_id org_ancestor_count max_generations
struct State
    display::Vector{Display}
    wp::Vector{Parameters}
    floats::Array{Float64,1}
    ints::Array{Int,1}
    resources::Array{Resource,1}
    organisms::Array{Organism,1}
    cur_max_ancestor::Array{Organism,1}
    helper::Vector{StateHelpers}
    function State(wp::Vector{Parameters})
        display=Vector{Display}(undef,1)
        display[1]=display_initialize()
        floats=Array{Float64,1}(undef,length(instances(StateFloatIndices)))
        ints=Array{Int,1}(undef,length(instances(StateIntIndices)))
        resources=Array{Resource,1}(undef,0)
        organisms=Array{Organism,1}(undef,0)
        cur_max_ancestor=Array{Organism,1}(undef,1)
        helper=Vector{StateHelpers}(undef,1)
        helper[1]=StateHelpers()
        new(display, wp, floats, ints, resources, organisms, cur_max_ancestor, helper)
    end
    function State(ws1::State, ws2::State)
        copyto!(ws1.wp,ws2.wp)
        copyto!(ws1.floats,ws2.floats)
        copyto!(ws1.ints,ws2.ints)
        empty!(ws1.resources)
        append!(ws1.resources,ws2.resources)
        empty!(ws1.organisms)
        append!(ws1.organisms,ws2.organisms)
        copyto!(ws1.cur_max_ancestor,ws2.cur_max_ancestor)
        copyto!(ws1.helper,ws2.helper)
    end
end

function set!(ws::State,floatIndex::StateFloatIndices,value::Float64)
    ws.floats[Int(floatIndex)]=value
end
function set!(ws::State,intIndex::StateIntIndices,value::Int)
    ws.ints[Int(intIndex)]=value
end
function get(ws::State,floatIndex::StateFloatIndices)
    ws.floats[Int(floatIndex)]
end
function get(ws::State,intIndex::StateIntIndices)
    ws.ints[Int(intIndex)]
end
function do_timestep!(ws::State)
    set!(ws,currentTime,nextfloat(get(ws,currentTime)))
    get(ws,currentTime)
end
function get_time(ws::State)
    get(ws,currentTime)
end
function set_maxsteps!(ws::State,value::Int)
    set!(ws,maxsteps,value)
end
function get_maxsteps(ws::State)
    get(ws,maxsteps)
end
function set_step!(ws::State,value::Int)
    set!(ws,step,value)
end
function get_step(ws::State)
    get(ws,step)
end
function do_step!(ws::State)
    set_step!(ws,get_step(ws)+1)
end

function clean(display)
    display_close(display)
end

function init!(param::Parameters,derived_param::DerivedParameters,ws::State)
    empty!(ws.resources)
    empty!(ws.organisms)
    ws.helper[1]=StateHelpers()
    set!(ws,currentTime,param.time)
    set!(ws,maxLifeTime,0.0)
    set!(ws,randomizeDelta,0.0)
    set!(ws,freeEnergy,param.energy)
    set!(ws,cells,param.cells)
    set!(ws,step,0)
    set!(ws,res_id,0)
    set!(ws,org_ancestor_count,0)
    set!(ws,max_generations,0)
    energy=get(ws,freeEnergy)
    while energy > 500
        #orgOrRes=rand(1:100)
        #if orgOrRes >= 5
            if ! add_single_resource(ws,derived_param)
                break
            end
            energy=get(ws,freeEnergy)
        #else
        #    if ! add_single_organism(ws,derived_param)
        #        break
        #    end
        #end
    end
end

function add_single_resource(ws::State,derived_param::DerivedParameters)
    energy=get(ws,freeEnergy)
    id=get(ws,res_id)
    res=create_random_resource(id,derived_param.possible_coordinates)
    energy-=res.type.energy
    if energy>=0
        push!(ws.resources,res)
        if haskey(ws.helper[1].res_pol_positions2resource,res.position)
            push!(ws.helper[1].res_pol_positions2resource[res.position],res)
        else
            ws.helper[1].res_pol_positions2resource[res.position]=[res]
        end
        set!(ws,res_id,id+1)
        set!(ws,freeEnergy,energy)
        draw_resource(ws.display[1],res)
        return true
    end
    return false
end

function test_resources(ws::State,wdp::DerivedParameters,N)
    #remove first N resources
    N = N>length(ws.resources) ? length(ws.resources) : N
    i=shuffle(1:length(ws.resources))[1:N]
    for res in ws.resources[i]
        energy=get(ws,freeEnergy)
        energy+=res.free_energy
        set!(ws,freeEnergy,energy)
        res.free_energy=0
        push!(ws.helper[1].resource2expunge,res)
        wipe_resource(ws.display[1],res)
    end
end

function expunge_resources(ws::State,img_width::Int,img_height::Int)
    for res in ws.helper[1].resource2expunge
        #res_pol_positions2resource::Dict{Coordinates,Array{Resource,1}}
        p=res.position
        if haskey(ws.helper[1].res_pol_positions2resource,p)
            resources=ws.helper[1].res_pol_positions2resource[p]
            index1=length(resources)
            depleted=zeros(Int,index1)
            index2=1
            for res2 in resources
                if res.id == res2.id
                    depleted[index1]=index2
                end
                index1-=1
                index2+=1
            end
            unique!(depleted)
            for index in depleted
                if index>0
                    deleteat!(resources,index)
                end
            end
        end
        #res_index_positions2resource::Dict{Coordinates2dIndex,Array{Resource,1}}
        pos2d_index=polar_to_2d_index(p,img_width,img_height)
        if haskey(ws.helper[1].res_index_positions2resource,pos2d_index)
            resources=ws.helper[1].res_index_positions2resource[pos2d_index]
            index1=length(resources)
            depleted=zeros(Int,index1)
            index2=1
            for res2 in resources
                if res.id == res2.id
                    depleted[index1]=index2
                end
                index1-=1
                index2+=1
            end
            unique!(depleted)
            for index in depleted
                if index>0
                    deleteat!(resources,index)
                end
            end
        end
        #resources::Array{Resource,1}
        index1=length(ws.resources)
        depleted=zeros(Int,index1)
        index2=1
        for res2 in ws.resources
            if res.id == res2.id
                depleted[index1]=index2
            end
            index1-=1
            index2+=1
        end
        unique!(depleted)
        for index in depleted
            if index>0
                deleteat!(ws.resources,index)
            end
        end
    end
    empty!(ws.helper[1].resource2expunge)
end

function add_single_organism(ws::State,derived_param::DerivedParameters)
    energy=get(ws,freeEnergy)
    cur_time=get(ws,currentTime)
    ancestor_count=get(ws,org_ancestor_count)
    org = Organism(derived_param.possible_coordinates,cur_time,ws.wp[1].org_max_energie)
    energy-=org.energy
    if energy>=0
        push!(ws.organisms,org)
        set!(ws,freeEnergy,energy)
        set!(ws,org_ancestor_count,ancestor_count+1)
        draw_org(ws.display[1],org)
        return true
    end
    return false
end

function split_organisms(ws::State,derived_param::DerivedParameters)
    siblings=Array{Organism,1}(undef,0)
    for org in ws.organisms
        if org.energy>1
            life_time=get(ws,currentTime)-org.birth_time
            split=0
            split_ok=true
            try
                split=org.trigger_split_function(life_time,org.trigger_split_expression_parameter)
            catch
                split_ok=false
            end
            if split_ok
                if split > ws.wp[1].split_boundaries[1] && split < ws.wp[1].split_boundaries[2]
                    split_organism!(org,ws,derived_param,siblings)
                end
            else
                org.energy=0
            end
        end
    end
    for org in siblings
        draw_org(ws.display[1],org)
        push!(ws.organisms,org)
    end
end

function split_organism!(org::Organism,ws::State,derived_param::DerivedParameters,siblings::Array{Organism,1})
    if org.energy>1
        new_energy1=floor(Int,org.energy/2)
        if new_energy1<=0
            println("split: new energy1: ",new_energy1)
        end
        new_energy2=org.energy-new_energy1
        if new_energy2<=0
            println("split: new energy2: ",new_energy2)
        end
        org.energy=new_energy1
        cur_time=get(ws,currentTime)
        sibling=Organism(org, true, get(ws,randomizeDelta) )
        sibling.energy=new_energy2
        sibling.is_sibling=true
        sibling.generation=org.generation+1
        sibling.current_position=get_random_coordinates(derived_param.possible_coordinates)
        if( get(ws,max_generations) < sibling.generation )
            set!(ws,max_generations,sibling.generation)
            ws.cur_max_ancestor[1]=deepcopy(org)
        end
        push!(siblings,sibling)
    end
end

function digest_external_energy_org(ws::State)
    res2wipe=Dict{Resource,Bool}()
    for organism in ws.organisms
        if haskey(ws.helper[1].res_pol_positions2resource,organism.current_position)
            external_energy=0
            for res in ws.helper[1].res_pol_positions2resource[organism.current_position]
                external_energy+=res.free_energy
            end
            if external_energy>0
                new_energy=0
                new_energy_ok=true
                try
                    new_energy=floor(Int,organism.energy_function(external_energy,organism.energy_expression_parameter))
                catch
                    new_energy_ok=false
                end
                if new_energy_ok
                    if new_energy>external_energy
                        new_energy=external_energy
                    end
                    if new_energy<0
                        new_energy=0
                    end
                    organism.energy+=new_energy
                    while new_energy>0
                        for res in ws.helper[1].res_pol_positions2resource[organism.current_position]
                            if res.free_energy>0
                                if res.free_energy>=new_energy
                                    res.free_energy-=new_energy
                                    new_energy=0
                                else
                                    new_energy-=res.free_energy
                                    res.free_energy=0
                                end
                            end
                            if res.free_energy==0
                                res2wipe[res]=true
                            end
                            if new_energy==0
                                break
                            end
                        end
                    end
                else
                    organism.energy=0
                end
            end
        end
    end
    for res in keys(res2wipe)
        wipe_resource(ws.display[1],res)
        push!(ws.helper[1].resource2expunge,res)
    end
end

function expunge_organisms(ws::State)
    index1=length(ws.organisms)
    dead=zeros(Int,index1)
    index2=1
    for organism in ws.organisms
        if organism.energy==0
            dead[index1]=index2
            cur_time=get(ws,currentTime)
            lifeTime=cur_time-organism.birth_time
            if lifeTime>get(ws,maxLifeTime)
                set!(ws,maxLifeTime,lifeTime)
            end
        end
        index1-=1
        index2+=1
    end
    for index in dead
        if index>0
            deleteat!(ws.organisms,index)
        end
    end
end

function housekeeping_organisms(ws::State)
    for organism in ws.organisms
        if organism.energy>0
            new_energy=0
            new_energy_ok=true
            try
                new_energy=floor(Int,organism.housekeeping_function(organism.energy,organism.housekeeping_expression_parameter))
            catch
                new_energy_ok=false
            end
            if new_energy_ok
                if new_energy>=organism.energy
                    new_energy=rand(1:organism.energy)-1
                end
                if new_energy<0
                    new_energy=0
                end
                world_energy=get(ws,freeEnergy)
                world_energy+=(organism.energy-new_energy)
                set!(ws,freeEnergy,world_energy)
                organism.energy=new_energy
            else
                organism.energy=0
            end
        end
        if organism.energy==0
            wipe_org(ws.display[1],organism)
        end
    end
end

function move_organisms(ws::State,wdp::DerivedParameters)
    for organism in ws.organisms
        if organism.energy>0
            move=true
            new_ϕ=0
            new_θ=0
            try
                new_ϕ = organism.ϕ_function(organism.current_position.ϕ,organism.ϕ_expression_parameter)
                new_θ = organism.θ_function(organism.current_position.θ,organism.θ_expression_parameter)
            catch
                move=false
            end
            if move
                organism.target_position = Coordinates(new_ϕ,new_θ)
                new_target_position = align_to_possible_coordinates(organism.target_position,wdp.possible_coordinates)
                if new_target_position.ϕ<0.0 || new_target_position.θ<0.0
                    println(organism.target_position," => ",new_target_position)
                end
                organism.target_position = new_target_position
                if organism.current_position!=organism.target_position
                    wipe_org(ws.display[1],organism)
                    organism.current_position = organism.target_position
                    draw_org(ws.display[1],organism)
                end
            else #remove organism
                organism.energy=0
            end
        end
    end
end

end


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 0eba536c-e941-11eb-153a-d7d86c5bab3d
begin
	import Pkg
	Pkg.activate(joinpath(@__DIR__))
	using SimLib, Plots, PlutoUI
	plotly()
	cd(joinpath(@__DIR__, ".."))
end

# ╔═╡ 2d6947b6-9306-48e3-b0fc-631c551b3f13
location = SaveLocation()

# ╔═╡ e13e2d66-25e5-47a9-8243-f66c84eecc53
pdd = PositionDataDescriptor(:noisy_chain_pbc, 1, 6, 30, [0.5, 1.0, 1.5], location)

# ╔═╡ a6dfa177-2d3c-410e-95ae-3b6eb1656de9
begin
	posdata = nothing
	with_terminal() do
		posdata = load_or_create(pdd);
	end
end

# ╔═╡ 4e4e8c35-af74-4a83-a3a6-18be7de97120
@bind i Slider(1:10)

# ╔═╡ 774f28cb-a09d-441b-b770-b8005f23ca6c
begin
	scatter(dropdims(posdata[:,:,i,:];dims=1), 1:posdata.system_size)
end

# ╔═╡ 19fbe000-a5cc-4804-839e-2dcf2d4f8f9b
posdata[:,:,1,1]

# ╔═╡ Cell order:
# ╠═0eba536c-e941-11eb-153a-d7d86c5bab3d
# ╠═2d6947b6-9306-48e3-b0fc-631c551b3f13
# ╠═e13e2d66-25e5-47a9-8243-f66c84eecc53
# ╠═a6dfa177-2d3c-410e-95ae-3b6eb1656de9
# ╠═4e4e8c35-af74-4a83-a3a6-18be7de97120
# ╠═774f28cb-a09d-441b-b770-b8005f23ca6c
# ╠═19fbe000-a5cc-4804-839e-2dcf2d4f8f9b

ext: jl


text:
# This file is a part of JuliaFEM
# License is MIT: see https://github.com/JuliaFEM/FEModels.jl/blob/master/LICENSE

using Base.Test
using FEModels

@testset "download file" begin
    fn = fe_download("block-and-cylinder", "mesh_sparse.inp", "50dbcbbcf4504c1a8746bc2a0655b69de4604276")
    @test isfile(fn)
end


text:
#using Morbit
#using Test

# Two standard (easy) test problems with known Pareto Data
# Many combinations of constaints, heterogenity, descent method and sampling algorithm are used

#using Logging
#global_logger(SimpleLogger(stdout, Logging.Debug))
#%%
@testset "1In1OutProblem_SteepestDescent" begin

    x0 = [5.0];
    f1(x) = x[1]^2;

    # I) unconstrained
    opt_settings = AlgoConfig(
        max_iter = 30
    );
    ## 1) treat objective as cheap
    mop = MixedMOP();
    add_objective!( mop, f1, :cheap );
    x,fx = optimize(mop, x0; algo_config = opt_settings);

    @test x ≈ [ 0.0 ] atol=1e-2

    ## 2) treat objective as expensive
    opt_settings = AlgoConfig(
        max_iter = 30
    );
    mop = MixedMOP();
    add_objective!( mop, f1, :expensive )
    x,fx = optimize(mop, x0; algo_config = opt_settings)

    @test x ≈ [ 0.0 ] atol = 1e-2

    #II) constrained

    ## 1) treat objective as cheap
    opt_settings = AlgoConfig(
        max_iter = 10
    );
    mop = MixedMOP( [-8.0], [8.0] );
    add_objective!( mop, f1, :cheap)
    x,fx = optimize(mop, x0; algo_config = opt_settings)

    @test x ≈ [ 0.0 ] atol=0.1

    ## 2) treat objective as expensive
    opt_settings = AlgoConfig(
        max_iter = 30
    );
    mop = MixedMOP( [-8.0], [8.0] );
    add_objective!( mop, f1, :expensive);
    x,fx = optimize(mop, x0; algo_config = opt_settings);

    @test x ≈ [ 0.0 ] atol=.1
end

#%%
@testset "2In2OutProblem_Steepest_Descent" begin
    x0 = [ π , π^2 ]
    ub = 10 .* ones(2);
    lb = -ub;

    g1(x) = sum( (x .- 1.0).^2 );
    g2(x) = sum( (x .+ 1.0).^2 );

    # unconstrained, cheap
    @testset "unbound_cheap_cheap" begin
        opt_settings = AlgoConfig(
            Δ_0= 0.2,
            max_iter = 30,
        );
        mop = MixedMOP();
        add_objective!(mop, g1, :cheap);
        add_objective!(mop, g2, :cheap);
        x,fx = optimize(mop, x0; algo_config = opt_settings);

        @test x[1] ≈ x[2] atol = .1
    end

    # unconstrained, expensive
    @testset "unbound_exp_exp" begin
        opt_settings = AlgoConfig(
            max_iter = 30
        );
        mop = MixedMOP();
        add_objective!(mop, g1, :expensive);
        add_objective!(mop, g2, :expensive);
        x,fx = optimize(mop, x0; algo_config = opt_settings);

        @test x[1] ≈ x[2] atol = .1
    end

    # unconstrained, heterogenous
    @testset "unbound_exp_cheap" begin
        opt_settings = AlgoConfig(
            max_iter = 30
        );
        mop = MixedMOP();
        add_objective!(mop, g1, :expensive);
        add_objective!(mop, g2, :cheap);
        x,fx = optimize(mop, x0; algo_config = opt_settings);

        @test x[1] ≈ x[2] atol = .1
    end

    # constrained, cheap
    @testset "bound_cheap_cheap" begin
        opt_settings = AlgoConfig(
            max_iter = 30
        );
        mop = MixedMOP(lb, ub);
        add_objective!(mop, g1, :cheap);
        add_objective!(mop, g2, :cheap);
        x,fx = optimize(mop, x0; algo_config = opt_settings);

        @test x[1] ≈ x[2] atol = .1
    end

    # constrained, expensive
    @testset "bound_exp_exp" begin
        opt_settings = AlgoConfig(
            max_iter = 50,
        );
        mop = MixedMOP(lb,ub);
        add_objective!(mop, g1, :expensive);
        add_objective!(mop, g2, :expensive);
        x,fx = optimize(mop, x0; algo_config = opt_settings);

        @test x[1] ≈ x[2] atol = .1
    end

    # constrained, heterogenous
    @testset "bound_exp_cheap" begin
        opt_settings = AlgoConfig(
            max_iter = 50,
            Δ_max = .4,
            _eps_crit = 0.1,
            strict_acceptance_test = true,
        );
        mop = MixedMOP(lb,ub);
        add_objective!(mop, g1, RbfConfig(kernel=:multiquadric));
        add_objective!(mop, g2, :cheap);
        x,fx = optimize(mop, x0; algo_config = opt_settings);

        @test x[1] ≈ x[2] atol = .1
    end
end
#%%
#=
@testset "direct_search" begin
    # 1D1D
    x0 = [5.0];
    f1(x) = x[1]^2;

    for lb_ub ∈ [ (Float32[],Float32[]), ([-5.0], [5.0]), ([-10.0], [10.0])]
        lb, ub = lb_ub
        for type ∈ [:cheap, :expensive ]
            for ideal_point ∈ [Float32[], [0.0]]
                opt_settings = AlgoConfig(
                    max_iter = 20,
                    descent_method = :direct_search,
                    ideal_point = ideal_point,
                )
                mop = MixedMOP(lb,ub)
                Morbit.scale( mop, x0)
                add_objective!(mop, f1, type)
                x, fx = optimize(mop, x0; algo_config = opt_settings)
                @test x[end] ≈ 0.0 atol = .1
            end
        end
    end

    # 2D2D
    #=
    x0 = [ π , π^2 ]
    ub = 10 .* ones(2);
    lb = -ub;

    g1(x) = sum( (x .- 1.0).^2 );
    g2(x) = sum( (x .+ 1.0).^2 );
    for lb_ub ∈ [ ( [],[] ), ([-10.0,-100], [10.0, 10.0])]
        lb,ub = lb_ub
        for type_tuple ∈ [ (:cheap, :cheap), (:expensive, :expensive), (:expensive, :cheap) ]
            type1, type2 = type_tuple
            for ideal_point ∈ [ [], zeros(2) ]
                opt_settings = AlgoConfig(
                    max_iter = 50,
                    rbf_kernel = :multiquadric,
                    rbf_shape_parameter = cs -> let Δ = cs.iter_data.Δ; return 1/(10*Δ) end,
                    Δ_0= 0.5,
                    all_objectives_descent = true,
                    use_max_points = type_tuple ==  (:expensive, :expensive) ? true : false,
                    max_model_points = 12,
                )
                mop = MixedMOP(lb,ub)
                add_objective!(mop, g1, type1)
                add_objective!(mop, g2, type2)
                x, fx = optimize(mop, x0; algo_config = opt_settings)
                @test x[1] ≈ x[2] atol = .5
            end
        end
    end
    =#
end
=#


text:
sun = Star("Kerbol", 261600000., 1.1723328e18, (254, 198, 20));
kerbin =    CelestialBody("Kerbin", 600000.,    3.5316000e12,       21549.425,          0.,             Orbit(13599840256.,     0.,             0.,                     0.,                     0.,             3.14,   0.,     sun),       (6,   154, 244))
mun =       CelestialBody("Mun",    200000.,    6.5138398e10,       138984.38,          deg2rad(230.),  Orbit(12000000.,        0.,             0.,                     0.,                     0.,             1.7,    0.,     kerbin),    (146, 150, 146))
minmus =    CelestialBody("Minmus", 60000.,     1765800026.,        40400.,             deg2rad(230.),  Orbit(47000000.,        0.,             deg2rad(6.),            deg2rad(78.),           deg2rad(38.),   0.9,    0.,     kerbin),    (160, 254, 176))
moho =      CelestialBody("Moho",   250000.,    168.609378655e9,    1210000.,           deg2rad(190.),  Orbit(5263138304.,      0.200000003,    deg2rad(7.),            deg2rad(70.),           deg2rad(15.),   3.14,   0.,     sun),       (210, 178, 110))
eve =       CelestialBody("Eve",    700000.,    8171.730229211e9,   80500.,             0.,             Orbit(9832684544.,      0.01,           deg2rad(2.099999905),   deg2rad(15.000000000),  deg2rad(0.),    3.14,   0.,     sun),       (200, 160, 240))
duna =      CelestialBody("Duna",   320000.,    301.363211975e9,    65517.859375,       deg2rad(90.),   Orbit(20726155264.,     0.051,          deg2rad(0.06),          deg2rad(135.5),         0.,             3.14,   0.,     sun),       (254, 66,  30 ))
ike =       CelestialBody("Ike",    130000.,    18568368573.,       65517.862134808,    0.,             Orbit(3200000.,         0.03,           deg2rad(0.2),           0.,                     0.,             1.7,    0.,     duna),      (168, 164, 150))
jool =      CelestialBody("Jool",   6000000.,   282528004209995.,   36000.,             0.,             Orbit(68773560320.,     0.05,           deg2rad(1.304),         deg2rad(52.),           deg2rad(0.),    0.1,    0.,     sun),       (126, 190, 2  ))
laythe =    CelestialBody("Laythe", 500000.,    1962000029236.,     52980.87905938,     deg2rad(90.),   Orbit(27184000.,        0.,             0.,                     0.,                     0.,             3.14,   0.,     jool),      (90,  126, 154))
vall =      CelestialBody("Vall",   300000.,    207481499474.,      105962.088893924,   0.,             Orbit(43152000.,        0.,             0.,                     0.,                     0.,             0.9,    0.,     jool),      (196, 254, 248))
bop =       CelestialBody("Bop",    65000.,     2486834944.,        544507.428516654,   deg2rad(230.),  Orbit(128500000.,       0.235,          deg2rad(15.),           deg2rad(10.),           deg2rad(25.),   0.9,    0.,     jool),      (226, 202, 118))
tylo =      CelestialBody("Tylo",   600000.,    2825280042100.,     211926.35802123,    0.,             Orbit(68500000.,        0.,             deg2rad(0.025),         deg2rad(0.),            deg2rad(0.),    3.14,   0.,     jool),      (216, 220, 214))
gilly =     CelestialBody("Gilly",  13000.,     8289450.,           28255.,             deg2rad(5.),    Orbit(31500000.,        0.55,           deg2rad(12.),           deg2rad(80.),           deg2rad(10.),   0.9,    0.,     eve),       (230, 218, 166))
pol =       CelestialBody("Pol",    44000.,     721702080.,         901902.623531173,   deg2rad(25.),   Orbit(179890000.,       0.17085,        deg2rad(4.25),          deg2rad(2.0),           deg2rad(15.),   0.9,    0.,     jool),      (214, 182, 10 ))
dres =      CelestialBody("Dres",   138000.,    21484488600.,       34800.,             deg2rad(25.),   Orbit(40839348203.,     0.145,          deg2rad(5.),            deg2rad(280.),          deg2rad(90.),   3.14,   0.,     sun),       (198, 202, 200))
eeloo =     CelestialBody("Eeloo",  210000.,    74410814527.,       19460.,             deg2rad(25.),   Orbit(90118820000.,     0.26,           deg2rad(6.15),          deg2rad(50.),           deg2rad(260.),  3.14,   0.,     sun),       (244, 244, 218))


text:
# export PrerecordedDriver

# type PrerecordedDriver <: DriverModel{NextState}
#     trajdata::Trajdata # log we pull from
#     time_start::Float64   # starting time [s]
#     t_current::Float64
#     timestep::Float64
#     nextstate::VehicleState
# end
# function PrerecordedDriver(trajdata::Trajdata, frame_start::Int, egoid::Int, timestep::Float64=get_timestep(trajdata))
#     time_start = get_time(trajdata, frame_start)
#     t_current = NaN
#     nextstate = get_state(trajdata, egoid, frame_start)
#     PrerecordedDriver(trajdata, time_start, t_current, timestep, nextstate)
# end

# get_name(::PrerecordedDriver) = "PrerecordedDriver"
# function reset_hidden_state!(model::PrerecordedDriver)
#     model.t_current = model.time_start
#     model
# end
# function observe!(model::PrerecordedDriver, scene::Scene, roadway::Roadway, egoid::Int)
#     model.t_current += model.timestep

#     i = searchsortedfirst(model.trajdata.frames, model.t_current, by=frame_or_time->isa(frame_or_time, TrajdataFrame) ? frame_or_time.t : frame_or_time) # index of first frame ≥ t_current

#     if i == 1
#         lo, hi = 1, 2
#     elseif i ≤ nframes(model.trajdata)
#         lo, hi = i-1, i
#     else # t_current > all times in model.trajdata.frames
#         lo, hi = i-2, i-1
#     end

#     t = model.t_current
#     t_lo = get_time(model.trajdata, lo)
#     t_hi = get_time(model.trajdata, hi)

#     γ = (t - t_lo) / (t_hi - t_lo)

#     if (γ == 0.0 && !iscarinframe(model.trajdata, egoid, lo)) ||
#        (γ == 1.0 && !iscarinframe(model.trajdata, egoid, hi)) ||
#        (!iscarinframe(model.trajdata, egoid, lo) || !iscarinframe(model.trajdata, egoid, hi))

#         model.nextstate = VehicleState()
#     else
#         state_lo = get_state(model.trajdata, egoid, lo)
#         state_hi = get_state(model.trajdata, egoid, hi)
#         model.nextstate = lerp(state_lo, state_hi, γ, model.trajdata.roadway)
#     end

#     model
# end

# Base.rand(model::PrerecordedDriver) = NextState(model.nextstate)
# Distributions.pdf(model::PrerecordedDriver, a::NextState) = a.s == model.nextstate ? Inf : 0.0
# Distributions.logpdf(model::PrerecordedDriver, a::NextState) = a.s == model.nextstate ? Inf : -Inf


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 8149551a-f4c6-11ea-252d-29f2baaad175
using Printf

# ╔═╡ 34cfe268-f4c7-11ea-1ec0-61e696684aa8
using Flux, Flux.Data.MNIST, Statistics

# ╔═╡ 34d0a516-f4c7-11ea-1e64-0d13a3840193
using Flux: onehotbatch, onecold, crossentropy, @epochs

# ╔═╡ 34d16fd2-f4c7-11ea-2f25-6362afceb48c
using Base.Iterators: partition

# ╔═╡ 34d9308c-f4c7-11ea-3401-43a143aa6f38
using BSON: @load, @save

# ╔═╡ 34da4742-f4c7-11ea-3add-c94710249ab8
using CuArrays

# ╔═╡ 34e01046-f4c7-11ea-3ec3-2528feadb9f0
using Random

# ╔═╡ e2020640-f4c6-11ea-32a6-e34d72594529
function prepare_dataset(;train=true)
    train_or_test = ifelse(train,:train,:test)
    imgs = MNIST.images(train_or_test)
    X = hcat(float.(vec.(imgs))...)
    labels = MNIST.labels(train_or_test)
    Y = onehotbatch(labels,0:9)
    return X, Y
end

# ╔═╡ e741e54e-f4c6-11ea-2e13-21ad98e74417
function define_model(;hidden)
    mlp = Chain(Dense(28^2,hidden,relu),
                Dense(hidden,hidden,relu),
                Dense(hidden,10),
                softmax)
    return mlp
end

# ╔═╡ ec9bf1e2-f4c6-11ea-1c77-87e5bc4d9d4a
function split_dataset_random(X, Y)
    divide_ratio=0.9
    shuffled_indices = shuffle(1:size(Y)[2])
    divide_idx = round(Int,0.9*length(shuffled_indices))
    train_indices = shuffled_indices[1:divide_idx]
    val_indices = shuffled_indices[divide_idx:end]
    train_X = X[:,train_indices]
    train_Y = Y[:,train_indices]
    val_X = X[:,val_indices]
    val_Y = Y[:,val_indices]
    return train_X, train_Y, val_X, val_Y
end

# ╔═╡ f7c15526-f4c6-11ea-0db3-13fff316d2c9
function train()
    println("Start to train")
    epochs = 10
    X, Y = prepare_dataset(train=true)
    train_X, train_Y, val_X,val_Y = split_dataset_random(X, Y)
    model = define_model(hidden=100) |> gpu
    loss(x,y)= crossentropy(model(x),y)
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    batchsize = 64
    train_dataset = gpu.([(train_X[:,batch] ,train_Y[:,batch]) for batch in partition(1:size(train_Y)[2],batchsize)])
    val_dataset = gpu.([(val_X[:,batch] ,val_Y[:,batch]) for batch in partition(1:size(val_Y)[2],batchsize)])

    callback_count = 0
    eval_callback = function callback()
        callback_count += 1
        if callback_count == length(train_dataset)
            println("action for each epoch")
            total_loss = 0
            total_acc = 0
            for (vx, vy) in val_dataset
                total_loss += loss(vx, vy)
                total_acc += accuracy(vx, vy)
            end
            total_loss /= length(val_dataset)
            total_acc /= length(val_dataset)
            @show total_loss, total_acc
            callback_count = 0
            pretrained = model |> cpu
            @save "pretrained.bson" pretrained
            callback_count = 0
        end
        if callback_count % 50 == 0
            progress = callback_count / length(train_dataset)
           @printf("%.3f\n", progress)
        end
    end
    optimizer = ADAM(params(model))

    @epochs epochs Flux.train!(loss, train_dataset, optimizer, cb = eval_callback)

    pretrained = model |> cpu
    weights = Tracker.data.(params(pretrained))
    @save "pretrained.bson" pretrained
    @save "weights.bson" weights
    println("Finished to train")
end

# ╔═╡ fd61af12-f4c6-11ea-254e-05caae039eb4
function predict()
    println("Start to evaluate testset")
    println("loading pretrained model")
    @load "pretrained.bson" pretrained
    model = pretrained |> gpu
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    println("prepare dataset")
    X, Y = prepare_dataset(train=false)
    X = X |> gpu
    Y = Y |> gpu
    @show accuracy(X, Y)
    println("Done")
end

# ╔═╡ 02de4586-f4c7-11ea-37e5-5fbc6d9ed260
function predict2()
    println("Start to evaluate testset")
    println("loading pretrained model")
    @load "weights.bson" weights
    model = define_model(hidden=100)
    Flux.loadparams!(model, weights)
    model = model |> gpu
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    println("prepare dataset")
    X, Y = prepare_dataset(train=false)
    X = X |> gpu
    Y = Y |> gpu
    @show accuracy(X, Y)
    println("Done")
end

# ╔═╡ 0a717b74-f4c7-11ea-3599-9384c8513b2a
train()

# ╔═╡ 1cb86220-f4c7-11ea-2e3d-571d34765494
predict()

# ╔═╡ 227c0586-f4c7-11ea-10db-bb45821ebb67
predict2()

# ╔═╡ Cell order:
# ╠═8149551a-f4c6-11ea-252d-29f2baaad175
# ╠═34cfe268-f4c7-11ea-1ec0-61e696684aa8
# ╠═34d0a516-f4c7-11ea-1e64-0d13a3840193
# ╠═34d16fd2-f4c7-11ea-2f25-6362afceb48c
# ╠═34d9308c-f4c7-11ea-3401-43a143aa6f38
# ╠═34da4742-f4c7-11ea-3add-c94710249ab8
# ╠═34e01046-f4c7-11ea-3ec3-2528feadb9f0
# ╠═e2020640-f4c6-11ea-32a6-e34d72594529
# ╠═e741e54e-f4c6-11ea-2e13-21ad98e74417
# ╠═ec9bf1e2-f4c6-11ea-1c77-87e5bc4d9d4a
# ╠═f7c15526-f4c6-11ea-0db3-13fff316d2c9
# ╠═fd61af12-f4c6-11ea-254e-05caae039eb4
# ╠═02de4586-f4c7-11ea-37e5-5fbc6d9ed260
# ╠═0a717b74-f4c7-11ea-3599-9384c8513b2a
# ╠═1cb86220-f4c7-11ea-2e3d-571d34765494
# ╠═227c0586-f4c7-11ea-10db-bb45821ebb67

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
using LinearAlgebra: ishermitian, issymmetric, mul!
using SparseArrays: AbstractSparseMatrix
using SimpleValueGraphs: OneEdgeValGraph
using Graphs: DefaultDistance

@testset "matrices" begin

    @testset "AdjacencyMatrix" begin

        for G      in (ValGraph, ValDiGraph),
            V      in TEST_VERTEX_TYPES_SMALL,
            E_VALS in TEST_EDGEVAL_TYPES_SMALL,
            gs     in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

            g = DummyValGraph(G{V, Tuple{}, E_VALS}(gs.graph;
                edgeval_init=(s, d) -> rand_sample(E_VALS)
            ))

            @testset "g::$(typeof(g))" begin

                a = adjacency_matrix(g)

                @testset "constructor" begin
                    @test typeof(a) == AdjacencyMatrix{typeof(g)}
                    @test a isa AbstractSparseMatrix{Bool, Int}
                end

                @testset "size(a)" begin
                     @test size(a) == (nv(g), nv(g))
                end

                @testset "getindex" begin
                    for u in 1:nv(g), v in 1:nv(g)
                        @test a[u, v] == has_edge(g, u, v)
                    end
                end

                if !is_directed(g)
                    @testset "ishermitian for undirected" begin
                        @test ishermitian(a)
                    end

                    @testset "issymmetric for undirected" begin
                        @test issymmetric(a)
                    end

                    @testset "adjoint is same matrix for undirected" begin
                        @test adjoint(a) === a
                    end

                    @testset "transpose is same matrix for undirected" begin
                        @test transpose(a) === a
                    end
                end

                @testset "ishermitian" begin
                    @test ishermitian(a) == ishermitian(Matrix(a))
                end

                @testset "issymmetric" begin
                    @test issymmetric(a) == issymmetric(Matrix(a))
                end

                @testset "adjoint" begin
                    @test adjoint(a) == adjoint(Matrix(a))
                end

                @testset "transpose" begin
                    @test transpose(a) == transpose(Matrix(a))
                end

            end
        end

    end

    @testset "mul! with AdjacencyMatrix and Vector" begin

        g1 = DummyValGraph(ValGraph{Int8}(6))
        add_edge!(g1, 2, 3)
        add_edge!(g1, 2, 4)
        add_edge!(g1, 5, 5)
        v1 = zeros(Float64, 6)
        mul!(v1, adjacency_matrix(g1), [1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
        @test v1 == [0.0, 7.0, 2.0, 2.0, 5.0, 0.0]

        g2 = DummyValGraph(ValDiGraph{Int16}(6))
        add_edge!(g2, 1, 1)
        add_edge!(g2, 2, 3)
        add_edge!(g2, 3, 4)
        add_edge!(g2, 4, 3)
        add_edge!(g2, 4, 5)
        v2 = zeros(Int16, 6)
        mul!(v2, adjacency_matrix(g2), [1, 2, 3, 4, 5, 6])
        @test v2 == [1, 3, 4, 8, 0, 0]
    end

    @testset "ValMatrix" begin
        for G      in (ValGraph, ValDiGraph),
            V      in TEST_VERTEX_TYPES_SMALL,
            E_VALS in TEST_EDGEVAL_TYPES_SMALL,
            gs     in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

            g = DummyValGraph(G{V, Tuple{}, E_VALS}(gs.graph;
                edgeval_init=(s, d) -> rand_sample(E_VALS)
            ))

            @testset "g::$(typeof(g))" begin
                for key in allkeys_for_E_VALS(E_VALS)

                    @testset "key = $key, zerovalue=$zv" for zv in (nothing, rand_sample(fieldtype(E_VALS, key)))

                        M = ValMatrix(g, key, zv)

                        @testset "size" begin
                            @test size(M) == (nv(g), nv(g))
                        end

                        @testset "getindex" begin
                            @test all(Iterators.product(vertices(g), vertices(g))) do (s, d)
                                M[s, d] == (has_edge(g, s, d) ? get_edgeval(g, s, d, key) : zv)
                            end
                        end

                        if eltype(M) <: Union{Number, Missing}
                            @testset "ishermitian" begin
                                 @test ishermitian(M) == ishermitian(Matrix(M))
                            end
                        end

                        if eltype(M)  <: Union{Number, Missing}
                            @testset "issymmetric" begin
                                @test issymmetric(M) == issymmetric(Matrix(M))
                            end
                        end

                        if eltype(M) <: Real && !is_directed(g)
                            @testset "adjoint for real, symmetric" begin
                                @test adjoint(M) === M
                            end
                        end

                        if eltype(M) <: Union{Number, Missing}
                            @testset "adjoint" begin
                                @test adjoint(M) == adjoint(Matrix(M))
                            end
                        end

                        if !is_directed(g)
                            @testset "transpose for symmetric" begin
                                @test transpose(M) === M
                            end
                        end

                        if eltype(M) <: Union{Number, Missing}
                            @testset "transpose" begin
                                @test transpose(M) == transpose(Matrix(M))
                            end
                        end

                        @testset "weights" begin

                            @testset "weights(g, $key; zerovalue=$zv)" begin
                                Mw = weights(g, key, zerovalue=zv)
                                @test typeof(M) == typeof(Mw)
                                @test M == Mw
                            end

                            if fieldtype(E_VALS, key) isa Number
                                @testset "weights(g, $key)" begin
                                    Mw = weights(g, key)
                                    @test M isa ValMatrix{fieldtype(E_VALS, key), typeof(g), key}
                                end
                            end

                            if g isa OneEdgeValGraph
                                @testset "weights(g; zerovalue=$zv)" begin
                                    Mw = weights(g, zerovalue=zv)
                                    @test M == Mw
                                end

                                if fieldtype(E_VALS, key) isa Number
                                    @testset "weights(g)" begin
                                        Mw = weights(g)
                                        @test M isa ValMatrix{fieldtype(E_VALS, key), typeof(g), key}
                                    end
                                end
                            end
                        end

                    end
                end
            end
        end
    end

    @testset "weights for value graphs without values" for
        G      in (ValGraph, ValDiGraph),
        V      in TEST_VERTEX_TYPES_SMALL,
        E_VALS in (Tuple{}, NamedTuple{(), Tuple{}}),
        gs     in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

        g = DummyValGraph(G{V, Tuple{}, E_VALS}(gs.graph;
            edgeval_init=(s, d) -> rand_sample(E_VALS))
        )

        @test weights(g) == DefaultDistance(nv(g))
    end

    @testset "weights for graph with two edge values and no key specified" begin

        @test_throws ArgumentError weights(DummyValGraph(ValGraph(3, edgeval_types=(a=Int, b=String))))
        @test_throws ArgumentError weights(DummyValGraph(ValDiGraph(2, edgeval_types=(Int, String))); zerovalue=nothing)
    end

    @testset "convert AdjacencyMatrix to SparseMatrixCSC" begin

        g1 = ValGraph{Int8}(4)
        add_edge!(g1, 1, 2)
        add_edge!(g1, 2, 3)
        add_edge!(g1, 3, 3)
        m1 = AdjacencyMatrix(g1)
        @test SparseMatrixCSC(m1) isa SparseMatrixCSC{Bool, Int}
        @test SparseMatrixCSC(m1) == [0 1 0 0; 1 0 1 0; 0 1 1 0; 0 0 0 0]

        g2 = ValDiGraph{Int16}(4; edgeval_types=(a=String, b=Int))
        add_edge!(g2, 1, 1, ("", 1))
        add_edge!(g2, 1, 2, ("", 2))
        add_edge!(g2, 1, 3, ("", 3))
        add_edge!(g2, 2, 1, ("", 3))
        add_edge!(g2, 3, 2, ("", 0))
        m2 = AdjacencyMatrix(g2)
        @test SparseMatrixCSC(m2) isa SparseMatrixCSC{Bool, Int}
        @test SparseMatrixCSC(m2) == [1 1 1 0; 1 0 0 0; 0 1 0 0; 0 0 0 0]

        g3 = ValOutDiGraph{UInt32}(4; vertexval_types=(String, Int), vertexval_init=undef)
        add_edge!(g3, 1, 1)
        add_edge!(g3, 1, 3)
        add_edge!(g3, 2, 3)
        add_edge!(g3, 3, 1)
        m3 = AdjacencyMatrix(g3)
        @test SparseMatrixCSC(m3) isa SparseMatrixCSC{Bool, Int}
        @test SparseMatrixCSC(m3) == [1 0 1 0; 0 0 1 0; 1 0 0 0; 0 0 0 0]

        g4 = DummyValGraph(ValGraph(0, graphvals=(1, 2, "xyz")))
        m4 = AdjacencyMatrix(g4)
        @test SparseMatrixCSC(m4) isa SparseMatrixCSC{Bool, Int}
        @test SparseMatrixCSC(m4) == Matrix{Bool}(undef, 0, 0)

    end

    @testset "convert ValMatrix to SparseMatrixCSC" begin

        g1 = ValGraph{Int8}(4, edgeval_types=(Int,))
        add_edge!(g1, 1, 2, (12,))
        add_edge!(g1, 2, 3, (23,))
        add_edge!(g1, 3, 3, (33,))
        m1 = ValMatrix(g1, 1, 0)
        @test SparseMatrixCSC(m1) isa SparseMatrixCSC{Int, Int}
        @test SparseMatrixCSC(m1) == [0 12 0 0; 12 0 23 0; 0 23 33 0; 0 0 0 0]

        g2 = ValDiGraph{Int16}(4; edgeval_types=(a=String, b=Int))
        add_edge!(g2, 1, 1, ("", 11))
        add_edge!(g2, 1, 2, ("", 12))
        add_edge!(g2, 1, 3, ("", 13))
        add_edge!(g2, 2, 1, ("", 21))
        add_edge!(g2, 3, 2, ("", 32))
        m2 = ValMatrix(g2, :b, 0)
        @test SparseMatrixCSC(m2) isa SparseMatrixCSC{Int, Int}
        @test SparseMatrixCSC(m2) == [11 12 13 0; 21 0 0 0; 0 32 0 0; 0 0 0 0]

        g3 = ValOutDiGraph{UInt32}(4; edgeval_types=(Int, Float64))
        add_edge!(g3, 1, 1, (11, 1.0))
        add_edge!(g3, 1, 3, (13, 1.3))
        add_edge!(g3, 2, 3, (23, 2.3))
        add_edge!(g3, 3, 1, (31, 3.1))
        m3 = ValMatrix(g3, 2, 0.0)
        @test SparseMatrixCSC(m3) isa SparseMatrixCSC{Float64, Int}
        @test SparseMatrixCSC(m3) == [1.0 0 1.3 0; 0 0 2.3 0; 3.1 0 0 0; 0 0 0 0]

        g4 = ValGraph(0, graphvals=(1, 2, "xyz"), edgeval_types=(Float64, String))
        m4 = ValMatrix(g4, 1, 0.0)
        @test SparseMatrixCSC(m4) isa SparseMatrixCSC{Float64, Int}
        @test SparseMatrixCSC(m4) == Matrix{Float64}(undef, 0, 0)

    end

end # testset


text:
using SIIPExamples
using InfrastructureSystems
using Logging

using Test

pkgpath = dirname(dirname(pathof(SIIPExamples)))
testpath = joinpath(pkgpath, "test")

exclude = ["04_bar_stack_plots.jl", "08_US-system-simulations.jl", "US_system.jl"]

logger = InfrastructureSystems.configure_logging(console_level = Logging.Error)

for (root, dirs, files) in walkdir(testpath)
    if root != testpath
        @testset "Testing Files in $root" begin
            for file in files
                if file in exclude
                    @info "skipping $file"
                else
                    fname = joinpath(root, file)
                    @testset "$file" begin
                        @test try
                            include(fname)
                            true
                        finally
                        end
                    end
                end
            end
        end
    end
end


text:
using Distributed
addprocs(4)
nprocs()

using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using Random
using LinearAlgebra
using DelimitedFiles
using GLM
using MendelGWAS
using CSV
using RCall
R"library(glmnet)"
R"require(pscl)"
R"require(boot)"

function iht_lasso_marginal(n::Int64, p::Int64, d::UnionAll, l::Link)
    #construct snpmatrix, covariate files, and true model b
    mafs = clamp!(0.5rand(p), 0.01, 0.5) #set lowest minor allele freq to be 0.01
    x = simulate_random_snparray(n, p, "tmp.bed", mafs=mafs)
    xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
    z = ones(n, 1) # the intercept
    x_float = [convert(Matrix{Float64}, x, center=true, scale=true) z] #Float64 version of x

    # simulate response, true model b, and the correct non-0 positions of b
    true_b = zeros(p)
    true_b[1:10] .= collect(0.1:0.1:1.0)
    shuffle!(true_b)
    correct_position = findall(x -> x != 0, true_b)

    #simulate phenotypes (e.g. vector y)
    if d == Normal || d == Poisson || d == Bernoulli
        prob = linkinv.(l, xbm * true_b)
        clamp!(prob, -20, 20)
        y = [rand(d(i)) for i in prob]
    elseif d == NegativeBinomial
        μ = linkinv.(l, xbm * true_b)
        prob = 1 ./ (1 .+ μ ./ nn)
        y = [rand(d(nn, i)) for i in prob] #number of failtures before nn success occurs
    elseif d == Gamma
        μ = linkinv.(l, xbm * true_b)
        β = 1 ./ μ # here β is the rate parameter for gamma distribution
        y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
    end
    y = Float64.(y)

    #specify path and folds
    num_folds = 3
    folds = rand(1:num_folds, size(x, 1));

    #run glmnet via Rcall
    @rput x_float y folds num_folds #make variables visible to R
    R"lasso_cv_result = cv.glmnet(x_float, y, nfolds = num_folds, foldid = folds, family='poisson')"
    R"lasso_beta_tmp = glmnet(x_float, y, lambda=lasso_cv_result$lambda.min, family='poisson')$beta"
    R"lasso_beta = as.vector(lasso_beta_tmp)"
    @rget lasso_cv_result lasso_beta #pull result from R to Julia
    lasso_k_est = count(!iszero, lasso_beta)

    #find non-zero entries returned by best lasso model as largest k estimate
    path = collect(1:50);

    #run IHT's cross validation routine
    mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, use_maf=false, debias=false, showinfo=false, parallel=true);
    iht_k_est = argmin(mses)
    iht_result = L0_reg(x, xbm, z, y, 1, iht_k_est, d(), l, debias=false, init=false, use_maf=false, show_info=false)
    iht_beta = iht_result.beta

    #Now run MendelGWAS
    make_bim_fam_files(x, y, "tmp") #create .bim and .bed files for MendelGWAS

    #create a control file for MendelGWAS
    open("tmp_control.txt", "w") do f
        write(f, "plink_input_basename = tmp \n")
        write(f, "plink_field_separator = '\t' \n")
        write(f, "output_table = tmp_table.txt \n")
        write(f, "regression = poisson \n")
        write(f, "regression_formula = Trait ~ \n")
    end

    #run marginal analysis
    GWAS("tmp_control.txt");

    # calculate false positive/negatives based on bonferroni correction
    p_values = CSV.read("tmp_table.txt", delim=',', header=true)[:Pvalue]
    significance = 0.05 / Float64(p)
    passing_snps_position = findall(p_values .<= significance)
    k = 10
    marginal_found = [correct_position[snp] in passing_snps_position for snp in 1:k]

    #clean up
    rm("tmp.bed", force=true)
    rm("tmp.bim", force=true)
    rm("tmp.fam", force=true)
    rm("tmp_table.txt", force=true)
    rm("tmp_control.txt", force=true)
    rm("Mendel_Output.txt", force=true)

    # to test if MendelGWAS is correct, fit using GLM.jl too (MendelGWAS is correct)
#     data = DataFrame(X=zeros(n), Y=y)
#     placeholder = zeros(n)
#     pvalues = zeros(p)
#     for i in 1:p
#         copyto!(placeholder, @view(x[:, i]), center=true, scale=true)
#         data.X .= placeholder
#         result = glm(@formula(Y ~ X), data, Poisson(), LogLink())
#         pvalues[i] = coeftable(result).cols[4][2].v
#     end
#     passing_snps_position2 = findall(pvalues .<= significance)
#     marginal_found2 = [correct_position[snp] in passing_snps_position2 for snp in 1:k]

    #Fit 0-inflated Poisson using R package because MendelGWAS has tons of false positives
    y = Int64.(y)
    @rput y n p
    R"
    zip_pvalues = matrix(0.0, p, 1)
    tmp = matrix(0, n, 1)
    for (i in 1:p) {
        tmp[,] = x_float[, i]
        zip_pvalues[i] = summary(zeroinfl(y ~ tmp))$coefficient[2]$zero[2,4]
    }
    "
    @rget zip_pvalues
    zip_pvalues = reshape(zip_pvalues, (10000,))
    zip_passing_snps = findall(zip_pvalues .<= significance)
    zip_found = [correct_position[snp] in zip_passing_snps for snp in 1:k]

    #show lasso and IHT's reconstruction result
    compare_model = DataFrame(
        true_β  = true_b[correct_position],
        IHT_β   = iht_beta[correct_position],
        lasso_β = lasso_beta[correct_position],
        marginal_found = marginal_found,
        zero_inf_Pois = zip_found)
    @show compare_model

    #compute true/false positives/negatives for IHT and lasso
    iht_tp = count(!iszero, iht_beta[correct_position])
    iht_fp = iht_k_est - iht_tp
    iht_fn = k - iht_tp
    lasso_tp = count(!iszero, lasso_beta[correct_position])
    lasso_fp = lasso_k_est - lasso_tp
    lasso_fn = k - lasso_tp
    marginal_tp = count(!iszero, true_b[passing_snps_position])
    marginal_fp = length(passing_snps_position) - marginal_tp
    marginal_fn = k - marginal_tp
#     juliaglm_tp = count(!iszero, true_b[passing_snps_position2])
#     juliaglm_fp = length(passing_snps_position2) - juliaglm_tp
#     juliaglm_fn = k - juliaglm_tp
    zip_tp = count(!iszero, true_b[zip_passing_snps])
    zip_fp = length(zip_passing_snps) - zip_tp
    zip_fn = k - zip_tp

    println("IHT true positives = $iht_tp")
    println("IHT false positives = $iht_fp")
    println("IHT false negatives = $iht_fn")
    println("LASSO true positives = $lasso_tp")
    println("LASSO false positives = $lasso_fp")
    println("LASSO false negatives = $lasso_fn")
    println("marginal true positives = $marginal_tp")
    println("marginal false positives = $marginal_fp")
    println("marginal false negatives = $marginal_fn")
#     println("julia glm true positives = $juliaglm_tp")
#     println("julia glm false positives = $juliaglm_fp")
#     println("julia glm false negatives = $juliaglm_fn" * "\n")
    println("zero inflated poisson true positives = $zip_tp")
    println("zero inflated poisson  false positives = $zip_fp")
    println("zero inflated poisson  false negatives = $zip_fn" * "\n")

    return iht_tp, iht_fp, iht_fn, lasso_tp, lasso_fp, lasso_fn,
        marginal_tp, marginal_fp, marginal_fn, zip_tp, zip_fp, zip_fn
end

function run()
    #simulat data with k true predictors, from distribution d and with link l.
    n = 1000
    p = 10000
    d = Poisson
    l = canonicallink(d())

    #run function above, saving results in 4 vectors
    total_runs = 50
    iht_true_positives = zeros(total_runs)
    iht_false_positives = zeros(total_runs)
    iht_false_negatives = zeros(total_runs)
    lasso_true_positives = zeros(total_runs)
    lasso_false_positives = zeros(total_runs)
    lasso_false_negatives = zeros(total_runs)
    marginal_true_positives = zeros(total_runs)
    marginal_false_positives = zeros(total_runs)
    marginal_false_negatives = zeros(total_runs)
    zip_true_positives = zeros(total_runs)
    zip_false_positives = zeros(total_runs)
    zip_false_negatives = zeros(total_runs)
    for i in 1:10
        println("current run = $i")

        #set random seed
        Random.seed!(i)

        iht_tp, iht_fp, iht_fn, lasso_tp, lasso_fp, lasso_fn,
            marginal_tp, marginal_fp, marginal_fn, zip_tp, zip_fp,
                zip_fn = iht_lasso_marginal(n, p, d, l)
        iht_true_positives[i] = iht_tp
        iht_false_positives[i] = iht_fp
        iht_false_negatives[i] = iht_fn
        lasso_true_positives[i] = lasso_tp
        lasso_false_positives[i] = lasso_fp
        lasso_false_negatives[i] = lasso_fn
        marginal_true_positives[i] = marginal_tp
        marginal_false_positives[i] = marginal_fp
        marginal_false_negatives[i] = marginal_fn
        zip_true_positives[i] = zip_tp
        zip_false_positives[i] = zip_fp
        zip_false_negatives[i] = zip_fn
    end

    return iht_true_positives, iht_false_positives, iht_false_negatives,
            lasso_true_positives, lasso_false_positives, lasso_false_negatives,
            marginal_true_positives, marginal_false_positives, marginal_false_negatives,
            zip_true_positives, zip_false_positives, zip_false_negatives
end

function run2()
	iht_true_positives, iht_false_positives, iht_false_negatives,
    lasso_true_positives, lasso_false_positives, lasso_false_negatives,
    marginal_true_positives, marginal_false_positives, marginal_false_negatives,
    zip_true_positives, zip_false_positives, zip_false_negatives = run()

    poisson_iht_true_positives = sum(iht_true_positives)
	poisson_iht_false_positives = sum(iht_false_positives)
	poisson_iht_false_negatives = sum(iht_false_negatives)
	poisson_lasso_true_positives = sum(lasso_true_positives)
	poisson_lasso_false_positives = sum(lasso_false_positives)
	poisson_lasso_false_negatives = sum(lasso_false_negatives)
	poisson_marginal_true_positives = sum(marginal_true_positives)
	poisson_marginal_false_positives = sum(marginal_false_positives)
	poisson_marginal_false_negatives = sum(marginal_false_negatives)
	poisson_zip_true_positives = sum(zip_true_positives)
	poisson_zip_false_positives = sum(zip_false_positives)
	poisson_zip_false_negatives = sum(zip_false_negatives)

	result = [poisson_iht_true_positives poisson_iht_false_positives poisson_iht_false_negatives;
        poisson_lasso_true_positives poisson_lasso_false_positives poisson_lasso_false_negatives;
        poisson_marginal_true_positives poisson_marginal_false_positives poisson_marginal_false_negatives;
        poisson_zip_true_positives poisson_zip_false_positives poisson_zip_false_negatives]

    writedlm("result1", result)
end

run2()


text:
struct HighsModel
   colcost
   collower
   colupper
   rowlower
   rowupper
   astart
   aindex
   avalue
end

struct HighsSolution
   colvalue
   coldual
   rowvalue
   rowdual
end

struct HighsBasis
   colbasisstatus
   rowbasisstatus
end

function Highs_call(model)
   n_col = convert(Cint, size(model.colcost, 1))
   n_row = convert(Cint, size(model.rowlower, 1))
   n_nz = convert(Cint, size(model.aindex, 1))

   colcost = convert(Array{Cdouble}, model.colcost)
   collower = convert(Array{Cdouble}, model.collower)
   colupper = convert(Array{Cdouble}, model.colupper)

   rowlower = convert(Array{Cdouble}, model.rowlower)
   rowupper = convert(Array{Cdouble}, model.rowupper)
   matstart = convert(Array{Cint}, model.astart)
   matindex = convert(Array{Cint}, model.aindex)
   matvalue = convert(Array{Cdouble}, model.avalue)

   solution = HighsSolution(Array{Cdouble, 1}(undef, n_col), Array{Cdouble, 1}(undef, n_col), Array{Cdouble, 1}(undef, n_row),  Array{Cdouble, 1}(undef, n_row))
   basis = HighsBasis(Array{Cint, 1}(undef, n_col), Array{Cint, 1}(undef, n_row))

   modelstatus = convert(Cint, 0)

   status = ccall((:Highs_call, "libhighs.so"), Cint, (Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ref{Cint}),
   n_col, n_row, n_nz, colcost, collower, colupper, rowlower, rowupper, matstart, matindex, matvalue, solution.colvalue, solution.coldual, solution.rowvalue, solution.rowdual, basis.colbasisstatus, basis.rowbasisstatus, modelstatus)

   return status, solution, basis, modelstatus
end

function Highs_create()
   return ccall((:Highs_create, "libhighs.so"), Ptr{Cvoid}, () )
end

function Highs_destroy(highs)
   ccall((:Highs_destroy, "libhighs.so"), Cvoid, (Ptr{Cvoid},), highs)
end

function Highs_run(highs)
   return ccall((:Highs_run, "libhighs.so"), Cint, (Ptr{Cvoid},), highs)
end

function Highs_readModel(highs, filename)
   name = convert(Cstring, pointer(filename))
   return ccall((:Highs_readModel, "libhighs.so"), Cint, (Ptr{Cvoid}, Cstring), highs, name)
end

function Highs_writeModel(highs, filename)
   name = convert(Cstring, pointer(filename))
   return ccall((:Highs_writeModel, "libhighs.so"), Cint, (Ptr{Cvoid}, Cstring), highs, name)
end

function Highs_passLp(highs, model)
   n_col = convert(Cint, size(model.colcost, 1))
   n_row = convert(Cint, size(model.rowlower, 1))
   n_nz = convert(Cint, size(model.aindex, 1))

   colcost = convert(Array{Cdouble}, model.colcost)
   collower = convert(Array{Cdouble}, model.collower)
   colupper = convert(Array{Cdouble}, model.colupper)

   rowlower = convert(Array{Cdouble}, model.rowlower)
   rowupper = convert(Array{Cdouble}, model.rowupper)
   matstart = convert(Array{Cint}, model.astart)
   matindex = convert(Array{Cint}, model.aindex)
   matvalue = convert(Array{Cdouble}, model.avalue)

   return ccall((:Highs_passLp, "libhighs.so"), Cint, (Ptr{Cvoid},Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},Ptr{Cint},Ptr{Cdouble}),
   highs, n_col, n_row, n_nz, colcost, collower, colupper, rowlower, rowupper, matstart, matindex, matvalue)
end

function Highs_setOptionValue(highs, option, value)
   opt = convert(Cstring, pointer(option))
   val = convert(Cstring, pointer(val))
   return ccall((:Highs_setOptionValue, "libhighs.so"), Cint, (Ptr{Cvoid}, Cstring, Cstring), highs, opt, val)
end

function Highs_getNumCols(highs)
   return ccall((:Highs_getNumCols, "libhighs.so"), Cint, ())
end

function Highs_getNumRows(highs)
   return ccall((:Highs_getNumRows, "libhighs.so"), Cint, ())
end

function Highs_getNumNz(highs)
   return ccall((:Highs_getNumNz, "libhighs.so"), Cint, ())
end

function Highs_getSolution(highs)
   numcol = Highs_getNumCols(highs)
   numrow = Highs_getNumRows(highs)

   solution = HighsSolution(Array{Cdouble, 1}(undef, numcol), Array{Cdouble, 1}(undef, numcol), Array{Cdouble, 1}(undef, numrow),  Array{Cdouble, 1}(undef, numrow))

   ccall((:Highs_getSolution, "libhighs.so"), Cvoid, (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), highs, solution.colvalue, solution.coldual, solution.rowvalue, solution.rowdual)

   return solution
end

function Highs_getBasis(highs)
   numcol = Highs_getNumCols(highs)
   numrow = Highs_getNumRows(highs)

   basis = HighsBasis(Array{Cint, 1}(undef, numcol), Array{Cint, 1}(undef, numrow))

   ccall((:Highs_getBasis, "libhighs.so"), Cvoid, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}), highs, basis.colbasisstatus, basis.rowbasisstatus)

   return basis
end

function Highs_getObjectiveValue(highs)
   return ccall((:Highs_getObjectiveValue, "libhighs.so"), Cdouble, (Ptr{Cvoid},), highs)
end

function Highs_getIterationCount(highs)
   return ccall((:Highs_getIterationCount, "libhighs.so"), Cint, (Ptr{Cvoid},), highs)
end

function Highs_addRow(highs, lower, upper, indices, values)
   n_new_nz = convert(Cint, size(indices, 1))

   lo = convert(Cdouble, lower)
   hi = convert(Cdouble, upper)

   idx = convert(Array{Cint}, indices)
   val = convert(Array{Cdouble}, values)

   return ccall((:Highs_addRow, "libhighs.so"), Cint, (Ptr{Cvoid}, Cdouble, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}), highs, lo, hi, n_new_nz, idx, val)
end

function Highs_addRows(highs, lower, upper, starts, indices, values)
   n_new_rows = convert(Cint, size(lower, 1))
   n_new_nz = convert(Cint, size(indices, 1))

   lo = convert(Array{Cdouble}, lower)
   hi = convert(Array{Cdouble}, upper)

   st = convert(Array{Cint}, starts)
   idx = convert(Array{Cint}, indices)
   val = convert(Array{Cdouble}, values)

   return ccall((:Highs_addRows, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
   highs, n_new_rows, lo, hi, n_new_nz, st, idx, val)
end

function Highs_addCol(highs, cost, lower, upper, indices, values)
   n_new_nz = convert(Cint, size(indices, 1))

   cc = convert(Cdouble, cost)
   lo = convert(Cdouble, lower)
   hi = convert(Cdouble, upper)

   idx = convert(Array{Cint}, indices)
   val = convert(Array{Cdouble}, values)

   return ccall((:Highs_addCol, "libhighs.so"), Cint, (Ptr{Cvoid}, Cdouble, Cdouble, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}),
   highs, cc, lo, hi, n_new_nz, idx, val)
end

function Highs_addCols(highs, costs, lower, upper, starts, indices, values)
   n_new_cols = convert(Cint, size(lower, 1))
   n_new_nz = convert(Cint, size(indices, 1))

   co = convert(Array{Cdouble}, costs)
   lo = convert(Array{Cdouble}, lower)
   hi = convert(Array{Cdouble}, upper)

   st = convert(Array{Cint}, starts)
   idx = convert(Array{Cint}, indices)
   val = convert(Array{Cdouble}, values)

   return ccall((:Highs_addCols, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
   highs, n_new_cols, co, lo, hi, n_new_nz, st, idx, val)
end

function Highs_changeObjectiveSense(highs, sense)
   sns = convert(Cint, sense)

   return ccall((:Highs_changeObjectiveSense, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint), highs, sns)
end

function Highs_changeColCost(highs, colidx, cost)
   col = convert(Cint, colidx)
   cst = convert(Cdouble, cost)

   return ccall((:Highs_changeColCost, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Cdouble), highs, colidx, cost)
end

function Highs_changeColsCostBySet(highs, set, cost)
   num_set_entries = convert(Cint, size(set, 1))

   st = convert(Array{Cint}, set)
   cst = convert(Array{Cdouble}, cost)

   return ccall((:Highs_changeColsCostBySet, "libhighs.so"), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}), highs, st, cst)
end

function Highs_changeColsCostByMask(highs, mask, cost)
   msk = convert(Array{Cint}, mask)
   cst = convert(Array{Cdouble}, cost)

   return ccall((:Highs_changeColsCostByMask, "libhighs.so"), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}), highs, msk, cst)
end

function Highs_changeColBounds(highs, col, lower, upper)
   colidx = convert(Cint, col)

   lo = convert(Cdouble, lower)
   hi = convert(Cdouble, upper)

   return ccall((:Highs_changeColBounds, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Cdouble, Cdouble), highs, colidx, lo, hi)
end

function Highs_changeColsBoundsByRange(highs, from, to, lower, upper)
   f = convert(Cint, from)
   t = convert(Int33, to)

   lo = convert(Array{Cdouble}, lower)
   hi = convert(Array{Cdouble}, upper)

   return ccall((:Highs_changeColsBoundsByRange, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}),
   highs, f, t, lo, hi)
end

function Highs_changeColsBoundsBySet(highs, set, lower, upper)
   nset = convert(Cint, size(set, 1))

   st = convert(Array{Cint}, set)

   lo = convert(Array{Cdouble}, lower)
   hi = convert(Array{Cdouble}, upper)

   return ccall((:Highs_changeColsBoundsBySet, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}),
   highs, nset, st, lo, hi)
end

function Highs_changeColsBoundsByMask(highs, mask, lower, upper)
   msk = convert(Array{Cint}, mask)

   lo = convert(Array{Cdouble}, lower)
   hi = convert(Array{Cdouble}, upper)

   return ccall((:Highs_changeColsBoundsByMask, "libhighs.so"), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}),
   highs, msk, lo, hi)
end

function Highs_changeRowBounds(highs, row, lower, upper)
   idx = convert(Cint, row)
   lo = convert(Cdouble, lower)
   hi = convert(Cdouble, upper)

   return ccall((:Highs_changeRowBounds, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Cdouble, Cdouble),
   highs, idx, lo, hi)
end

function Highs_changeRowsBoundsBySet(highs, set, lower, upper)
   nst = convert(Cint, size(set, 1))
   st = convert(Array{Cint}, set)

   lo = convert(Array{Cdouble}, lower)
   hi = convert(Array(Cdouble), upper)

   return ccall((:Highs_changeRowsBoundsBySet, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}),
   highs, nst, st, lo, hi)
end

function Highs_changeRowsBoundsByMask(highs, mask, lower, upper)
   msk = convert(Array{Cint}, mask)

   lo = convert(Array{Cdouble}, lower)
   hi = convert(Array{Cdouble}, upper)

   return ccall((:Highs_changeRowsBoundsByMask, "libhighs.so"), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}),
   highs, msk, lo, hi)
end

#=

int Highs_getColsByRange(
    void *highs,          //!< HiGHS object reference
    const int from_col,   //!< The index of the first column to
                          //!< get from the model
    const int to_col,     //!< One more than the last column to get
                          //!< from the model
    int* num_col,          //!< Number of columns got from the model
    double *costs,        //!< Array of size num_col with costs
    double *lower,        //!< Array of size num_col with lower bounds
    double *upper,        //!< Array of size num_col with upper bounds
    int* num_nz,           //!< Number of nonzeros got from the model
    int *matrix_start,    //!< Array of size num_col with start
                          //!< indices of the columns
    int *matrix_index,    //!< Array of size num_nz with row
                          //!< indices for the columns
    double *matrix_value  //!< Array of size num_nz with row
                          //!< values for the columns
);

int Highs_getColsBySet(
    void *highs,                //!< HiGHS object reference
    const int num_set_entries,  //!< The number of indides in the set
    const int *set,             //!< Array of size num_set_entries with indices
                                //!< of columns to get
    int* num_col,                //!< Number of columns got from the model
    double *costs,              //!< Array of size num_col with costs
    double *lower,              //!< Array of size num_col with lower bounds
    double *upper,              //!< Array of size num_col with upper bounds
    int* num_nz,                 //!< Number of nonzeros got from the model
    int *matrix_start,          //!< Array of size num_col with start indices
                                //!< of the columns
    int *matrix_index,          //!< Array of size num_nz with row indices
                                //!< for the columns
    double *matrix_value        //!< Array of size num_nz with row values
                                //!< for the columns
);

int Highs_getColsByMask(
    void *highs,          //!< HiGHS object reference
    const int *mask,      //!< Full length array with 1 => get; 0 => not
    int* num_col,          //!< Number of columns got from the model
    double *costs,        //!< Array of size num_col with costs
    double *lower,        //!< Array of size num_col with lower bounds
    double *upper,        //!< Array of size num_col with upper bounds
    int* num_nz,           //!< Number of nonzeros got from the model
    int *matrix_start,    //!<  Array of size num_col with start
                          //!<  indices of the columns
    int *matrix_index,    //!<  Array of size num_nz with row indices
                          //!<  for the columns
    double *matrix_value  //!<  Array of size num_nz with row values
                          //!<  for the columns
);

int Highs_getRowsByRange(
    void *highs,          //!< HiGHS object reference
    const int from_row,   //!< The index of the first row to get from the model
    const int to_row,     //!< One more than the last row get from the model
    int* num_row,          //!< Number of rows got from the model
    double *lower,        //!< Array of size num_row with lower bounds
    double *upper,        //!< Array of size num_row with upper bounds
    int* num_nz,           //!< Number of nonzeros got from the model
    int *matrix_start,    //!< Array of size num_row with start indices of the
                          //!< rows
    int *matrix_index,    //!< Array of size num_nz with column indices for the
                          //!< rows
    double *matrix_value  //!< Array of size num_nz with column values for the
                          //!< rows
);

int Highs_getRowsBySet(
    void *highs,                //!< HiGHS object reference
    const int num_set_entries,  //!< The number of indides in the set
    const int *set,             //!< Array of size num_set_entries with indices
                                //!< of rows to get
    int* num_row,                //!< Number of rows got from the model
    double *lower,              //!< Array of size num_row with lower bounds
    double *upper,              //!< Array of size num_row with upper bounds
    int* num_nz,                 //!< Number of nonzeros got from the model
    int *matrix_start,          //!< Array of size num_row with start indices
                                //!< of the rows
    int *matrix_index,          //!< Array of size num_nz with column indices
                                //!< for the rows
    double *matrix_value        //!< Array of size num_nz with column
                                //!< values for the rows
);

int Highs_getRowsByMask(
    void *highs,          //!< HiGHS object reference
    const int *mask,      //!< Full length array with 1 => get; 0 => not
    int* num_row,          //!< Number of rows got from the model
    double *lower,        //!< Array of size num_row with lower bounds
    double *upper,        //!< Array of size num_row with upper bounds
    int* num_nz,           //!< Number of nonzeros got from the model
    int *matrix_start,    //!< Array of size num_row with start indices
                          //!< of the rows
    int *matrix_index,    //!< Array of size num_nz with column indices
                          //!< for the rows
    double *matrix_value  //!< Array of size num_nz with column
                          //!< values for the rows
);

=#

function Highs_deleteColsByRange(highs, from, to)
   f = convert(Cint, from)
   t = convert(Cint, to)

   return ccall((:Highs_deleteColsByRange, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Cint), highs, f, t)
end

function Highs_deleteColsBySet(highs, set)
   nset = convert(Cint, size(set, 1))

   st = convert(Array{Cint}, set)

   return ccall((:Highs_deleteColsBySet, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}), highs, nset, st)
end

function Highs_deleteColsByMask(highs, mask)
   msk = convert(Array{Cint}, mask)

   return ccall((:Highs_deleteColsByMask, "libighs.so"), Cint, (Ptr{Cvoid}, Ptr{Cint}), highs, msk)
end

function Highs_deleteRowsByRange(highs, from, to)
   f = convert(Cint, from)
   t = convert(Cint, to)

   return ccall((:Highs_deleteRowsByRange, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Cint), highs, f, t)
end

function Highs_deleteRowsBySet(highs, set)
   nset = convert(Cint, size(set, 1))

   st = convert(Array{Cint}, set)

   return ccall((:Highs_deleteRowsBySet, "libhighs.so"), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}), highs, nset, st)
end

function Highs_deleteRowsByMask(highs, mask)
   msk = convert(Array{Cint}, mask)

   return ccall((:Highs_deleteRowsByMask, "libhighs.so"), Cint, (Ptr{Cvoid}, Ptr{Cint}), highs, msk)
end

max_stars_repo_name: ycollet/HiGHS

max_issues_repo_name: ycollet/HiGHS

max_forks_repo_name: ycollet/HiGHS


text:
struct ListNode{T}
    value::T
    next::Union{ListNode{T}, Nothing}
end

mutable struct LinkedList{T}
    head::Union{ListNode{T}, Nothing}
end

LinkedList(T::Type) = LinkedList{T}(nothing)

Base.iterate(ll::LinkedList) = ll.head === nothing ? nothing : (ll.head.value, ll.head)
Base.iterate(ll::LinkedList{T}, state::ListNode{T}) where T =
    state.next === nothing ? nothing : (state.next.value, state.next)

function Base.getindex(ll::LinkedList, idx::Integer)
    idx < 1 && throw(BoundsError("$idx is less than 1"))
    for v in ll
        idx -= 1
        idx == 0 && return v
    end
    throw(BoundsError("index beyond end of linked list"))
end

function Base.pushfirst!(ll::LinkedList{T}, items::T...) where T
    for item in reverse(items)
        ll.head = ListNode{T}(item, ll.head)
    end
    ll
end

Base.show(io::IO, ll::LinkedList{T}) where T =
    print(io, "LinkedList{$T}[" * join(ll, ", ") * "]")

Base.eltype(ll::LinkedList{T}) where T = T

Base.length(ll::LinkedList) = count(v -> true, ll)

Base.firstindex(ll::LinkedList) = 1
Base.lastindex(ll::LinkedList) = length(ll)

max_stars_repo_path: Chapter06/03. Defining your own types – linked list/ll.jl


text:
module M1_1
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    3.38058845865336e-11*l1.^2 - 8.57996768441923e-11*l1.*l2 + 1.77343591248946e-11*l1.*l3 + 6.96600978936586e-14*l1 + 5.44544999423272e-11*l2.^2 - 2.25334113136583e-11*l2.*l3 - 8.83600620216189e-14*l2 + 2.33986794159514e-12*l3.^2 + 1.82331574873355e-14*l3 + 3.59116069801844e-17
end
end


text:
# Algorithm 3.2
# Newton's interpolatory divided difference.

function dividedDiff(x, y)
	n = length(x)
	out = zeros(n, n)
	for i ∈ 1:n
		out[i, 1] = y[i]
	end
	for i ∈ 2:n
		for j ∈ 2:i
			out[i, j] = (out[i, j - 1] - out[i - 1, j - 1]) / (x[i] - x[i - j + 1])
		end
	end
	return out
end

println(dividedDiff([1, 1.3, 1.6, 1.9, 2.2], [-0.76, -0.54, -0.57, -0.57, 0.3]))


text:
module QuasiArrays
using Base, LinearAlgebra, LazyArrays
import Base: getindex, size, axes, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
                Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, copy, copyto!, zero,
                map, eachindex, eltype, first, last, firstindex, lastindex, in, reshape, all,
                isreal, iszero, isempty, empty, isapprox, fill!, getproperty
import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinline_meta,
                DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type,
                _maybetail, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,
                index_shape, to_shape, unsafe_length, @nloops, @ncall, unalias, _unaliascopy,
                to_index, to_indices, _to_subscript_indices, _splatmap, dataids
import Base: ViewIndex, Slice, IdentityUnitRange, ScalarIndex, RangeIndex, view, viewindexing, ensure_indexable, index_dimsum,
                check_parent_index_match, reindex, _isdisjoint, unsafe_indices, _unsafe_ind2sub,
                _ind2sub, _sub2ind, _ind2sub_recurse, _lookup, SubArray,
                parentindices, reverse, ndims, checkbounds,
                promote_shape, maybeview, checkindex, checkbounds_indices,
                throw_boundserror, rdims, replace_in_print_matrix, show,
                hcat, vcat, hvcat
import Base: *, /, \, +, -, ^, inv
import Base: exp, log, sqrt,
          cos, sin, tan, csc, sec, cot,
          cosh, sinh, tanh, csch, sech, coth,
          acos, asin, atan, acsc, asec, acot,
          acosh, asinh, atanh, acsch, asech, acoth
import Base: Array, Matrix, Vector

import Base.Broadcast: materialize, materialize!, BroadcastStyle, AbstractArrayStyle, Style, broadcasted, Broadcasted, Unknown,
                        newindex, broadcastable, preprocess, _eachindex, _broadcast_getindex,
                        DefaultArrayStyle, axistype, throwdm, instantiate, combine_eltypes, eltypes

import LinearAlgebra: transpose, adjoint, checkeltype_adjoint, checkeltype_transpose, Diagonal,
                        AbstractTriangular, pinv, inv, promote_leaf_eltypes, power_by_squaring,
                        integerpow, schurpow, tr

import LazyArrays: MemoryLayout, UnknownLayout, Mul, ApplyLayout, BroadcastLayout, ⋆,
                    lmaterialize, _lmaterialize, InvOrPInv, ApplyStyle, AbstractLazyLayout, LazyLayout,
                    FlattenMulStyle, IdentityMulStyle, MulAddStyle, LazyArrayApplyStyle,
                    Applied, flatten, _flatten, arguments, _mat_mul_arguments, _vec_mul_arguments,
                    rowsupport, colsupport, tuple_type_memorylayouts, applylayout, broadcastlayout,
                    LdivApplyStyle, most, InvLayout, PInvLayout, sub_materialize,
                    _mul, rowsupport, DiagonalLayout, adjointlayout, transposelayout, conjlayout,
                    sublayout, call, combine_mul_styles, result_mul_style, LazyArrayStyle, lazy_getindex

import Base.IteratorsMD

export AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector, materialize,
       QuasiArray, QuasiMatrix, QuasiVector, QuasiDiagonal, Inclusion,
       QuasiAdjoint, QuasiTranspose, ApplyQuasiArray, ApplyQuasiMatrix, ApplyQuasiVector,
       BroadcastQuasiArray, BroadcastQuasiMatrix, BroadcastQuasiVector

if VERSION < v"1.3-"
    """
    broadcast_preserving_zero_d(f, As...)

    Like [`broadcast`](@ref), except in the case of a 0-dimensional result where it returns a 0-dimensional container

    Broadcast automatically unwraps zero-dimensional results to be just the element itself,
    but in some cases it is necessary to always return a container — even in the 0-dimensional case.
    """
    function broadcast_preserving_zero_d(f, As...)
        bc = broadcasted(f, As...)
        r = materialize(bc)
        return length(axes(bc)) == 0 ? fill!(similar(bc, typeof(r)), r) : r
    end
    broadcast_preserving_zero_d(f) = fill(f())
    broadcast_preserving_zero_d(f, as::Number...) = fill(f(as...))
else
    import Base.Broadcast: broadcast_preserving_zero_d
end

abstract type AbstractQuasiArray{T,N} end
AbstractQuasiVector{T} = AbstractQuasiArray{T,1}
AbstractQuasiMatrix{T} = AbstractQuasiArray{T,2}
AbstractQuasiVecOrMat{T} = Union{AbstractQuasiVector{T}, AbstractQuasiMatrix{T}}
const AbstractQuasiOrVector{T} = Union{AbstractVector{T},AbstractQuasiVector{T}}
const AbstractQuasiOrMatrix{T} = Union{AbstractMatrix{T},AbstractQuasiMatrix{T}}
const AbstractQuasiOrArray{T} = Union{AbstractArray{T},AbstractQuasiArray{T}}

cardinality(d) = length(d)

size(A::AbstractQuasiArray) = map(cardinality, axes(A))
axes(A::AbstractQuasiArray) = error("Override axes for $(typeof(A))")

include("indices.jl")
include("abstractquasiarray.jl")
include("multidimensional.jl")
include("subquasiarray.jl")
include("quasireshapedarray.jl")
include("quasibroadcast.jl")
include("abstractquasiarraymath.jl")

include("quasiarray.jl")
include("quasiarraymath.jl")

include("lazyquasiarrays.jl")
include("quasiconcat.jl")

include("matmul.jl")
include("inv.jl")
include("quasiadjtrans.jl")
include("quasidiagonal.jl")
include("dense.jl")

promote_leaf_eltypes(x::AbstractQuasiArray{T}) where {T<:Number} = T
promote_leaf_eltypes(x::AbstractQuasiArray) = mapreduce(promote_leaf_eltypes, promote_type, x; init=Bool)

function isapprox(x::AbstractQuasiArray, y::AbstractQuasiArray;
    atol::Real=0,
    rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x),promote_leaf_eltypes(y),atol),
    nans::Bool=false, norm::Function=norm)
    d = norm(x - y)
    if isfinite(d)
        return d <= max(atol, rtol*max(norm(x), norm(y)))
    else
        # Fall back to a component-wise approximate comparison
        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol, nans=nans), zip(x, y))
    end
end

end


text:
using DataFrames
using CSV
using Plots
#%%
# run(`./pi2`)
x=-1:0.1:5
f(x)=exp(-x)-x
y=f.(x)

plot(x,y)
xaxis!("x")
yaxis!("y")

#%%
using Pkg

Pkg.add("ClusterManagers")

#%%
using ClusterManagers
using Distributed

n_processes=5
ENV["HOSTNAME"]="127.0.0.1"
ClusterManagers.addprocs_htc(n_processes)

#%%
@distributed for i=1:n_processes
       run(`hostname`)
       end

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
using KernelFunctions

SUITE["KernelFunctions"] = BenchmarkGroup()

kernelnames = ["SqExponentialKernel"]
kerneltypes = ["ARD", "ISO"]
kernels = Dict{String,Dict{String,KernelFunctions.Kernel}}()
for k in kernelnames
    kernels[k] = Dict{String,KernelFunctions.Kernel}()
    SUITE["KernelFunctions"][k] = BenchmarkGroup()
    for kt in kerneltypes
        SUITE["KernelFunctions"][k][kt] = BenchmarkGroup()
        kernels[k][kt] = eval(
            Meta.parse(
                "KernelFunctions." *
                k *
                "(" *
                (kt == "ARD" ? "alpha*ones(Float64,dim)" : "alpha") *
                ")",
            ),
        )
    end
end

for k in kernelnames
    for kt in kerneltypes
        SUITE["KernelFunctions"][k][kt]["k(X,Y)"] = @benchmarkable KernelFunctions.kernelmatrix(
            $(kernels[k][kt]), $X, $Y; obsdim=1
        )
        # SUITE["KernelFunctions"][k][kt]["k!(X,Y)"] = @benchmarkable KernelFunctions.kernelmatrix!(KXY,$(kernels[k][kt]),$X,$Y) setup=(KXY=copy($KXY))
        SUITE["KernelFunctions"][k][kt]["k(X)"] = @benchmarkable KernelFunctions.kernelmatrix(
            $(kernels[k][kt]), $X; obsdim=1
        )
        # SUITE["KernelFunctions"][k][kt]["k!(X)"] = @benchmarkable KernelFunctions.kernelmatrix!(KX,$(kernels[k][kt]),$X) setup=(KX=copy($KX))
        # SUITE["KernelFunctions"][k][kt]["kdiag(X)"] = @benchmarkable KernelFunctions.kernelmatrix_diag($(kernels[k][kt]),$X)
        # SUITE["KernelFunctions"][k][kt]["kdiag!(X)"] = @benchmarkable KernelFunctions.kernelmatrix_diag!(kX,$(kernels[k][kt]),$X) setup=(kX=copy($kX))
    end
end
# results = run(SUITE)


text:
const DEFAULT_COLL_NAME = "jingpengw_test"
const DEFAULT_EXP_NAME = "test"
const DEFAULT_CHAN_NAME = "em"
const DEFAULT_DATA_TYPE = UInt8
const DEFAULT_RESOLUTION_LEVEL = 0
const DEFAULT_BOSSAPI_VERSION = "v1"
const DEFAULT_ARRAY_DIMENSION = 3

struct BOSSArray{T, N} <: AbstractArray{T,N}
    urlPrefix       :: String
    collectionName  :: String
    experimentName  :: String
    channelName     :: String
    resolutionLevel :: Int
    headers         :: Dict{String, String}
end

function BOSSArray(
         foo             ::Type{T},
         N               ::Int,
         urlPrefix       ::String,
         collectionName  ::String,
         experimentName  ::String,
         channelName     ::String,
         resolutionLevel ::Int,
         headers         ::Dict{String, String} ) where T
    BOSSArray{T,N}(   urlPrefix, collectionName, experimentName,
                channelName, resolutionLevel, headers)
end

# @generated function (::Type{BOSSArray{T,N}}){T,N}(
#                         T               ::DataType,
#                         N               ::Int,
#                         urlPrefix       ::String,
#                         collectionName  ::String,
#                         experimentName  ::String,
#                         channelName     ::String,
#                         resolutionLevel ::Int,
#                         headers         ::Dict{String, String} )
#     return quote
#         $(Expr(:meta, :inline))
#         BOSSArray{T,N}( urlPrefix, collectionName, experimentName,
#                         channelName, resolutionLevel, headers )
#     end
# end

function BOSSArray(;
                    T               ::DataType  = DEFAULT_DATA_TYPE,
                    N               ::Int       = DEFAULT_ARRAY_DIMENSION,
                    collectionName  ::String    = DEFAULT_COLL_NAME,
                    experimentName  ::String    = DEFAULT_EXP_NAME,
                    channelName     ::String    = DEFAULT_CHAN_NAME,
                    resolutionLevel ::Int       = DEFAULT_RESOLUTION_LEVEL)
    urlPrefix = "$(ENV["INTERN_PROTOCOL"])://$(ENV["INTERN_HOST"])"*
                    "/$(DEFAULT_BOSSAPI_VERSION)"
    headers = Dict("Authorization"  => "Token $(ENV["INTERN_TOKEN"])",
                    "content-type"  => "application/blosc",
                    "Accept"        => "application/blosc")

    BOSSArray{T,N}( urlPrefix, collectionName, experimentName, channelName,
                resolutionLevel, headers )
end


text:
using StackOverflow

mutable struct SiteInfo
    new_active_users::Int
    total_users::Int
    badges_per_minute::Float64
    total_badges::Int
    total_votes::Int
    total_comments::Int
    answers_per_minute::Float64
    questions_per_minute::Float64
    total_answers::Int
    total_accepted::Int
    total_unanswered::Int
    total_questions::Int
end

"""
    getsiteinfo(;site::String = "stackoverflow")

Will return a collection of statistics about the site (StackOverflow).
"""
function getsiteinfo(;site::String = "stackoverflow")

    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/info?site=$(site)")

    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)

    newitem = json_obj["items"]

    for item in newitem

        new_active_users = 0
        total_users = 0
        badges_per_minute = 0
        total_badges = 0
        total_votes = 0
        total_comments = 0
        answers_per_minute = 0
        questions_per_minute = 0
        total_answers = 0
        total_accepted = 0
        total_unanswered = 0
        total_questions = 0

        for (key, value) in item
            if occursin("total_users", key)
                total_users = value
            elseif occursin("total_badges", key)
                total_badges = value
            elseif occursin("total_answers", key)
                total_answers = value
            elseif occursin("total_comments", key)
                total_comments = value
            elseif occursin("questions_per_minute", key)
                questions_per_minute = value
            elseif occursin("total_accepted", key)
                total_accepted = value
            elseif occursin("total_votes", key)
                total_votes = value
            elseif occursin("badges_per_minute", key)
                badges_per_minute = value
            elseif occursin("total_unanswered", key)
                total_unanswered = value
            elseif occursin("new_active_users", key)
                new_active_users = value
            elseif occursin("answers_per_minute", key)
                answers_per_minute = value
            elseif occursin("total_questions", key)
                total_questions = value
            end
        end

        site_info = SiteInfo(new_active_users, total_users, badges_per_minute, total_badges, total_votes, total_comments, answers_per_minute, questions_per_minute, total_answers, total_accepted, total_unanswered, total_questions)

        return site_info
    end
end


text:
#=
    http://cplint.eu/example/inference/dice.swinb

    """
    A six-sided die is repeatedly thrown until the outcome is six. on(T,F)
    means that on the Tth throw the face F came out.
    """

    What is the probability that the die lands on face 1 at time 0?

    Cf ~/blog/dice_6_throws.blog
     ~/psi/dice_6_throws.psi
     ~/webppl/dice_6_throws.wppl

     Also see the related problem rolling_dice4.jl where we investigate
     the probability of the length of this (stopping) sequence.

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

# ix: Throw index
# val: The value to check
@model function dice_6_throws(ix=1,val=1)

    max_len = 100
    function throws(a)
        # Note: If t is a random variable then we got StackOverflowError
        # t ~  DiscreteUniform(1,6) # Don't work!
        # Hence we draw a "plain" random value
        t = rand(DiscreteUniform(1,6))
        if t == 6
            return a
        else
            return throws(vcat(a,t))
        end
    end

    a = throws(Int16[])
    len ~ Dirac(length(a))

    # Probability that the die lands on face val on throw ix
    test ~ Dirac(length(a) >= ix && a[ix] == val)
end

function run_model(ix, val)
    model = dice_6_throws(ix,val)

    num_chns = 4

    # chns = sample(model, Prior(), 10_000)

    # chns = sample(model, MH(), MCMCThreads(), 100_000, num_chns)
    # chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)

    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, MH(), 1_000)

    # chns = sample(model, PG(15), MCMCThreads(), 1_000, num_chns)

    # chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)
    # chns = sample(model, SMC(1000), 10_000)
    chns = sample(model, IS(), 10_000)

    display(chns)
    # display(plot(chns))

    show_var_dist_pct(chns,:test)
    show_var_dist_pct(chns,:len)

end

for i in 1:10
    println("\nix:$i val:1")
    run_model(i,1)
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Jonas August 1/2022

#=
This problem was asked by Google.

Given an array of integers where every integer occurs three times except for one integer, which only occurs once, find and return the non-duplicated integer.

For example, given [6, 1, 3, 3, 3, 6, 6], return 1. Given [13, 19, 13, 13], return 19.

Do this in O(N) time and O(1) space.
=#

#=
Assume that the array is correctly formed (one single among triples; no doubles, quadruples, etc.; only integers).

There are only at most (n-1)/3 + 1 unique numbers.

n = 3*m+1 is odd if m even and vice versa.

=#

#=
If y is the single entry, then the sum of the remaining entries is divisible by three.  Let x\y denote set x with y removed.

y is unique in in x => (sum of z s.t. z in x\y) mod 3 = 0

This is necessary condition for identifying y, but is it sufficient?  No!  Consider x=[13,19,13,13].  Try y=13.  Then sum of x\y is 19+13+13 = 45 which is divisible by 3.  Argh!

But maybe there's some function of the elements that can be used instead?

If a,b are integers, then cube root of aab is integer iff b = a?

Relevance Fermat's last theorem?
x^k + y^k = z^k has no positive integer solutions x,y,z for integer k>2.
Try k = 3 and cube entries?

=#

"Find that integer in input array x that has no match, while all other elements come in three, in no particular order.  Takes O(n) time, O(1) space, where n = length of the array."
function single_among_triples_fast(x::Vector{Int})::Int
    # TODO: allow zeros and negatives
    if any(x .== 0) @error("No zeroes allowed") end
    if any(x .< 0) @error("No negatives allowed") end
    p = 1
    for y in x
        p *= y
    end
    for y in x
        # TODO: CHECK that this always works
        c = (p/y)^(1/3)
        # TODO: what tolerance ensures correct results?
        if abs(c-round(c)) < sqrt(eps(c))    # is c approximately an integer?
            return y
        end
    end
    @error("This line should never be reached for well-formed input.")
end

"Find that integer in input array x that has no match, while all other elements come in three, in no particular order.  Takes O(n^2) time, O(1) space, where n = length of the array."
function single_among_triples_brute_force(x::Vector{Int})::Int
    for i in 1:length(x)
       matched = false
       for j in i+1:length(x)
           if x[j] == x[i]
	       matched = true
	       break
	   end
       end
       if !matched
           return x[i]
       end
    end
    @error("This line should never be reached for well-formed input.")
end

#=
Other False starts:

Not O(1) space:

- Convert to linked list.  Pick first element and scan for match and delete matches if found, then start again.  Must end up with unmatched number at beginning.

x[i[1]] is unique in x[:]
x[i[2]] == x[i[3]] == x[i[4]]
x[i[5]] == x[i[6]] == x[i[7]]
...

x[i[1]] + sum(j = 2:3:n) 3*x[i[j]]

Not O(n) time:

- Sort.  Triples will cluster, so linearly scan for the single.

=#

# TODO: Construct random inputs

using Test

for single_among_triples in [single_among_triples_brute_force, single_among_triples_fast]
  # @test @test single_among_triples([]) == ?  not defined? exception?
  @test single_among_triples([19]) == 19
  @test single_among_triples([13, 19, 13, 13]) == 19
  @test single_among_triples([6, 1, 3, 3, 3, 6, 6]) == 1
end

max_stars_repo_path: Single among triples - Problem 40.jl

max_issues_repo_path: Single among triples - Problem 40.jl

max_forks_repo_path: Single among triples - Problem 40.jl