text:
using JuMP
import JSON
using Distributions
include("prob_1.jl")
include("prob.jl")

function run_heuristic(; file = "../data/nesta_case24_ieee_rts_nk.m", k = 4)

    data = PowerModels.parse_file(file)
    buses = Dict([(parse(Int, k), v) for (k, v) in data["bus"]])
    buses = filter((i, bus) -> bus["bus_type"] != 4, buses)
    branches = Dict([(parse(Int, k), v) for (k, v) in data["branch"]])
    branches = filter((i, branch) -> branch["br_status"] == 1 && branch["f_bus"] in keys(buses) && branch["t_bus"] in keys(buses), branches)
    branch_indexes = collect(keys(branches))

    log_p = Dict([(i, log(branches[i]["prob"])) for i in branch_indexes])
    ref = PowerModels.build_ref(data)

    successful_attacks = Attack[]
    for branch in branches
        gm = Graph()
        gm.data = data
        gm.ref = ref
        create_n_1_graph(gm, buses, branches, branch)
        isolated_load_vertices = check_graph_for_isolation(gm, buses, branches)
        for load in isolated_load_vertices
            push!(load.branch_ids, branch[1])
            push!(successful_attacks, load)
        end
    end

    remove_redundant_attacks(successful_attacks)
    time = Dict(i => 0.0 for i in 2:k)

    for i in 2:k
        if det_or_stoch == 1
            attack = create_and_solve_det_milp(data, ref, successful_attacks, k = i)
        else
            tic()
            attack = create_and_solve_stoch_milp(data, ref, successful_attacks, log_p, k = i)
            time[i] = toq()
        end
        push!(successful_attacks, attack)
    end

    f = open("./output_files/$(length(buses))", "w")
    write(f, "case k h_time h_obj_value probability nf nf_e dc dc_e soc soc_e lines\n")

    attack_sol = []
    for i in 1:length(successful_attacks)
        if length(successful_attacks[i].branch_ids) == 1
            continue
        else
            attack_dict = Dict("k" => length(successful_attacks[i].branch_ids),
                               "load" => successful_attacks[i].total_load,
                               "branch_ids" => successful_attacks[i].branch_ids)

            write(f, "$(length(buses)) $(attack_dict["k"]) $(time[attack_dict["k"]]) $(successful_attacks[i].total_load) ")
            actual_prob = 1.0

            # solve nf model
            for index in successful_attacks[i].branch_ids
                actual_prob *= data["branch"][string(index)]["prob"]
                data["branch"][string(index)]["br_status"] = 0
            end
            write(f, "$(actual_prob) ")
            pm = NFPowerModel(data; solver = CplexSolver())
            post_pfls(pm)
            status, solve_time = solve(pm)
            result = PowerModels.build_solution(pm, status, solve_time)
            println("full NF load shed ... $(result["objective"])")
            println("expected NF load shed ... $(actual_prob*result["objective"])")
            println("###################################")
            solution = Dict{AbstractString,Any}("k" => k, "lines" => collect(successful_attacks[i].branch_ids),
                                                "model" => "nf", "time" => time[attack_dict["k"]],
                                                "prob" => actual_prob, "load_shed" => result["objective"],
                                                "expected_load_shed" => actual_prob*result["objective"])
            write(f, "$(solution["load_shed"]) $(solution["expected_load_shed"]) ")

            # solve dc model
            pm = DCPPowerModel(data; solver = CplexSolver())
            post_pfls(pm)
            status, solve_time = solve(pm)
            result = PowerModels.build_solution(pm, status, solve_time)
            println("full DC load shed ... $(result["objective"])")
            println("expected DC load shed ... $(actual_prob*result["objective"])")
            println("###################################")
            solution = Dict{AbstractString,Any}("k" => k, "lines" => collect(successful_attacks[i].branch_ids),
                                                "model" => "nf", "time" => time[attack_dict["k"]],
                                                "prob" => actual_prob, "load_shed" => result["objective"],
                                                "expected_load_shed" => actual_prob*result["objective"])
            write(f, "$(solution["load_shed"]) $(solution["expected_load_shed"]) ")

            # solve soc model
            pm = SOCWRPowerModel(data; solver = CplexSolver())
            post_pfls(pm)
            status, solve_time = solve(pm)
            result = PowerModels.build_solution(pm, status, solve_time)
            println("full SOC load shed ... $(result["objective"])")
            println("expected SOC load shed ... $(actual_prob*result["objective"])")
            println("###################################")
            solution = Dict{AbstractString,Any}("k" => k, "lines" => collect(successful_attacks[i].branch_ids),
                                                "model" => "nf", "time" => time[attack_dict["k"]],
                                                "prob" => actual_prob, "load_shed" => result["objective"],
                                                "expected_load_shed" => actual_prob*result["objective"])
            write(f, "$(solution["load_shed"]) $(solution["expected_load_shed"]) ")

            #= solve ac model
            pm = ACPPowerModel(data; solver = IpoptSolver())
            post_pfls(pm)
            status, solve_time = solve(pm)
            result = PowerModels.build_solution(pm, status, solve_time)
            println("full AC load shed ... $(result["objective"])")
            println("expected AC load shed ... $(actual_prob*result["objective"])")
            println("###################################")
            solution = Dict{AbstractString,Any}("k" => k, "lines" => collect(successful_attacks[i].branch_ids),
                                                "model" => "nf", "time" => time[attack_dict["k"]],
                                                "prob" => actual_prob, "load_shed" => result["objective"],
                                                "expected_load_shed" => actual_prob*result["objective"])
            write(f, "$(solution["load_shed"]) $(solution["expected_load_shed"]) ")
            =#
            write(f, "$(successful_attacks[i].branch_ids)\n")

            push!(attack_sol, attack_dict)
            # println("$(length(successful_attacks[i].branch_ids)) ... $(successful_attacks[i].total_load); branch_ids ... $([i for i in successful_attacks[i].branch_ids]); bus_ids ... $([i for i in successful_attacks[i].bus_ids])")
        end
    end
    #=
    json_dict = Dict{AbstractString,Any}()
    json_dict["attack"] = attack_sol
    json_string = JSON.json(json_dict)
    if det_or_stoch == 1
        write("det_isolation.json", json_string)
    else
        write("stoch_isolation.json", json_string)
    end
    =#
end

file = String(ARGS[1])
k = parse(Int, ARGS[2])
det_or_stoch = parse(Int, ARGS[3])
run_heuristic(file = file, k = k)

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
using BinDeps

@BinDeps.setup

zmq = library_dependency("zmq", aliases = ["libzmq"])

provides(Sources,URI("http://download.zeromq.org/zeromq-3.2.3.tar.gz"),zmq)
provides(BuildProcess,Autotools(libtarget = "src/.libs/libzmq."*BinDeps.shlib_ext),zmq)

@osx_only begin
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")  end
    using Homebrew
    provides( Homebrew.HB, "zeromq", zmq, os = :Darwin )
end

@BinDeps.install

ext: jl


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

using Test

# It needs to be called first to trigger the crash.
include("issue980.jl")

# Tests for solvers are located in MOI.Test.

include("dummy.jl")

# MOI tests not relying on any submodule
@testset "MOI" begin
    include("isbits.jl")
    include("isapprox.jl")
    include("interval.jl")
    include("errors.jl")
    include("functions.jl")
    include("sets.jl")
    include("attributes.jl")
    include("instantiate.jl")
end

# Utilities submodule tests
@testset "MOI.$(submodule)" for submodule in [
    "Bridges",
    "FileFormats",
    "Test",
    "Utilities",
    "Benchmarks",
]
    include("$(submodule)/$(submodule).jl")
end

# Test hygiene of @model macro
include("hygiene.jl")


text:
struct BufferCreateInfo
    mHandleRef::Ref{vk.VkBufferCreateInfo}

    function BufferCreateInfo(;
                pNext = C_NULL, #::Ptr{Cvoid}
                flags = 0, #::VkBufferCreateFlags
                size = 0, #::VkDeviceSize
                usage = 0, #::VkBufferUsageFlags
                sharingMode = vk.VK_SHARING_MODE_EXCLUSIVE, #::VkSharingMode
                queueFamilyIndexCount = 0, #::UInt32
                pQueueFamilyIndices = C_NULL #::Ptr{UInt32}
            )
        this = new(Ref(vk.VkBufferCreateInfo(
            vk.VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO, #sType::VkStructureType
            pNext, #::Ptr{Cvoid}
            flags, #::VkBufferCreateFlags
            size, #::VkDeviceSize
            usage, #::VkBufferUsageFlags
            sharingMode, #::VkSharingMode
            queueFamilyIndexCount, #::UInt32
            pQueueFamilyIndices, #::Ptr{UInt32}
        )))
    end
end
@class BufferCreateInfo

function handleRef(this::BufferCreateInfo)::Ref{vk.VkBufferCreateInfo}
    return this.mHandleRef
end

function stagingBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_TRANSFER_SRC_BIT)
end

function storageBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_STORAGE_BUFFER_BIT |
                                                 vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT |
                                                 vk.VK_BUFFER_USAGE_TRANSFER_SRC_BIT)
end

function arrayBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT |
                                                 vk.VK_BUFFER_USAGE_VERTEX_BUFFER_BIT)
end

function indexBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT |
                                                 vk.VK_BUFFER_USAGE_INDEX_BUFFER_BIT)
end

function indexArrayBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT |
                                                 vk.VK_BUFFER_USAGE_VERTEX_BUFFER_BIT |
                                                 vk.VK_BUFFER_USAGE_INDEX_BUFFER_BIT)
end

function downloadBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT)
end

function uniformBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT |
                                                 vk.VK_BUFFER_USAGE_UNIFORM_BUFFER_BIT)
end

function raytracingBuffer(size::Csize_t = Csize_t(0))::BufferCreateInfo
    return BufferCreateInfo(size = size, usage = vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT |
                                                 vk.VK_BUFFER_USAGE_RAY_TRACING_BIT_NV)
end

function copyWithSize(this::BufferCreateInfo, size::Csize_t)
    handle = this.mHandleRef[]
    return BufferCreateInfo(
        pNext = handle.pNext, #::Ptr{Cvoid}
        flags = handle.flags, #::VkBufferCreateFlags
        size = size, #::VkDeviceSize
        usage = handle.usage, #::VkBufferUsageFlags
        sharingMode = handle.sharingMode, #::VkSharingMode
        queueFamilyIndexCount = handle.queueFamilyIndexCount, #::UInt32
        pQueueFamilyIndices = handle.pQueueFamilyIndices #::Ptr{UInt32}
    )
end

function copyWithUsage(this::BufferCreateInfo, usage::vk.VkBufferUsageFlags)
    handle = this.mHandleRef[]
    return BufferCreateInfo(
        pNext = handle.pNext, #::Ptr{Cvoid}
        flags = handle.flags, #::VkBufferCreateFlags
        size = handle.size, #::VkDeviceSize
        usage = usage, #::VkBufferUsageFlags
        sharingMode = handle.sharingMode, #::VkSharingMode
        queueFamilyIndexCount = handle.queueFamilyIndexCount, #::UInt32
        pQueueFamilyIndices = handle.pQueueFamilyIndices #::Ptr{UInt32}
    )
end


text:
# HZAM_Sym_Julia_V1.jl
# First public release of HZAM_Sym_Julia
# by Darren Irwin, begun in early July 2021
# Re-written from the R version of HZAM_Sym

# Please note that this version is not intended to be super user-friendly software
# for people who aren't familiar with the Julia language. Rather, I am making this available
# so that those who know the basics of Julia can run this and tinker if they desire,
# and so that the community can see exactly what code produced the results in the paper.

# If you use this code, please cite the paper below, and possibly my GitHub or Dryad repository
# where you obtained this code:

# Irwin, D., and Schluter, D. Hybridization and the coexistence of species. bioRxiv 2021.04.04.438369; doi: https://doi.org/10.1101/2021.04.04.438369

# You need to make sure these packages below are added to your Julia environment.
# to add, type "]" to get the "pkg>" prompt, then type e.g. "add Distributions";
# or do this:
# import Pkg;
# Pkg.add("Distributions")
# Pkg.add("Statistics")
# Pkg.add("JLD2")
# Pkg.add("CSV")
# Pkg.add("DataFrames")
# Pkg.add("Plots")
# Pkg.add("CategoricalArrays")
# Pkg.add("Colors")
# Pkg.add("ColorSchemes")
# Pkg.add("CategoricalArrays")

using Distributions # needed for "Poisson" function
using Statistics  # needed for "mean" function
using JLD2 # needed for saving / loading data in Julia format
using CSV # for saving in csv format
using DataFrames # for converting data to save as csv

# for plotting:
using Plots
gr()  # use GR backend for graphs
using CategoricalArrays
using Colors, ColorSchemes
import ColorSchemes.plasma
using Plots.PlotMeasures  # needed for plot margin adjustment

# to start Julia with multiple threads, type in terminal e.g.:
# julia --threads 4
# To check, type in Julia: Threads.nthreads()
# 4

# set up functions (should not define functions repeatedly in loop, as causes re-compilation, slows things)

function generate_genotype_array(N_pop0,N_pop1,loci)
    total_N = N_pop0 + N_pop1
    genotypes = Array{Int8, 3}(undef, 2, loci, total_N) # The "Int8" is the type (8-bit integer), and "undef" means an unitialized array, so values are meaningless
    genotypes[:,:,1:N_pop0] .= 0  # assigns genotypes of pop1
    genotypes[:,:,(N_pop0+1):total_N] .= 1  # assigns genotypes of pop
    return genotypes
end

function calc_traits_additive(genotypes)
    N = size(genotypes, 3)
    traits = Array{Float32, 1}(undef, N) # Float32 should be enough precision; memory saving compared to Float64
    for i in 1:N
        traits[i] = mean(genotypes[:,:,i])
    end
    return traits
end

function get_survival_fitnesses_epistasis(genotypes, w_hyb, beta=1)
    survival_HI = calc_traits_additive(genotypes)
    epistasis_fitnesses = 1 .- (1 - w_hyb) .* (4 .* survival_HI .* (1 .- survival_HI)).^beta
    return epistasis_fitnesses
end

function get_survival_fitnesses_hetdisadvantage(genotypes, w_hyb)
    N = size(genotypes, 3)
    num_loci = size(genotypes, 2)
    s_per_locus = 1 - w_hyb ^ (1/num_loci)  # loss in fitness due to each heterozygous locus
    num_hetloci = Array{Int16, 1}(undef, N)
    for ind in 1:N  # count number of het loci per individual
        num_hetloci[ind] = sum(genotypes[1,:,ind] .!= genotypes[2,:,ind])
    end
    hetdisadvantage_fitnesses = (1 - s_per_locus) .^ num_hetloci
    return hetdisadvantage_fitnesses
end

function run_HZAM(set_name::String, ecolDiff, intrinsic_R, replications;  # the semicolon makes the following optional keyword arguments
    K_total::Int = 1000, max_generations::Int = 1000,
    total_loci::Int = 6, female_mating_trait_loci = 1:3, male_mating_trait_loci = 1:3,
    competition_trait_loci = 1:3, hybrid_survival_loci = 1:3, neutral_loci = 4:6,
    survival_fitness_method::String = "epistasis", per_reject_cost = 0,
    starting_pop_ratio = 1.0)
    # ecolDiff should be from 0 to 1 (parameter "E" in the paper)
    # intrinsic_R is called "R" in the paper
    # replications should be somehting like "1:10" or just "1" for 1 replicate, or something like "2:5" to add replicates after 1 is done
    # K_total (default 1000) is the total carrying capacity, divided equally between K_alpha and K_beta
    # max_generations (default 1000) is the number of generations each simulation will run
    # total_loci (default 6) is the total number of loci, regardless of their role (with indices 1:total_loci, referred to below)
    # Need to specify indices (columns) of four types of functional loci (can be the same). At least one should begin with index 1:
    # female_mating_trait_loci (default 1:3) is indices of the loci that determine the female mating trait
    # male_mating_trait_loci (default 1:3) is indices of the loci that determine the male mating trait
    # competition_trait_loci (default 1:3) is indices of the loci that determine the ecological trait (used in fitness related to resource use)
    # hybrid_survival_loci (default 1:3) is indices of the loci that determine survival probability of offspring to adulthood (can be viewed as incompatibilities and/or fitness valley based on ecology)
    # Specify indices (columns) of neutral loci (which have no effect on anything, just along for the ride):
    # neutral_loci (default 4:6) is indices of neutral loci (used for neutral measure of hybrid index; not used in the HZAM-sym paper)
    # per_reject_cost (default 0, can take values of 0 to 1) is fitness loss of female per male rejected (due to search time, etc.)
    # starting_pop_ratio (default 1) is the starting ratio of pop A to pop B

    save_outcomes_JL = true
    save_outcomes_csv = true  # whether to save the whole outcome array as csv files (with each rep as separate file)
    save_each_sim = false  # whether to save detailed data for each simulation

    # the set of hybrid fitnesses (w_hyb) values that will be run
    w_hyb_set = [1, 0.98, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0] # for just one run, just put one number in this and next line
    # the set of assortative mating strengths (S_AM) that will be run
    S_AM_set = [1, 3, 10, 30, 100, 300, 1000, Inf]  # ratio of: probably of accepting homospecific vs. prob of accepting heterospecific

    # get the chosen survival fitness function
    if survival_fitness_method == "epistasis"
        get_survival_fitnesses = get_survival_fitnesses_epistasis
        short_survFitnessMethod = "Ep"
    elseif survival_fitness_method == "hetdisadvantage"
        get_survival_fitnesses = get_survival_fitnesses_hetdisadvantage
        short_survFitnessMethod = "Het"
    else
        println("ERROR--no survival fitness method chosen--should be either epistasis or hetdisadvantage")
    end

    total_functional_loci = max(maximum(female_mating_trait_loci), maximum(male_mating_trait_loci), maximum(competition_trait_loci), maximum(hybrid_survival_loci))
    functional_loci_range = 1:total_functional_loci
    num_neutral_loci = length(Vector(neutral_loci))
    if total_functional_loci + num_neutral_loci ≠ total_loci
        println("#### WARNING: Please examine your loci numbers and indices, as they don't all match up ####")
    end

    # intrinsic_R = 1.05  # Intrinsic growth rate, this is the average maximum expected number of offspring per individual, when pop size far below K
    K_A = K_total / 2  # EVEN NUMBER; carrying capacity (on resource alpha) of entire range (for two sexes combined), regardless of species
    K_B = K_total / 2   # EVEN NUMBER; carrying capacity (on resource beta) of entire range (for two sexes combined), regardless of species

    pop0_starting_N = K_A   # starting N of species 0
    pop0_starting_N_half = Int(pop0_starting_N/2)
    pop1_starting_N = Int(round(starting_pop_ratio * K_B))   # starting N of species 1, which can be lower if starting_pop_ratio is below 1)
    pop1_starting_N_half = Int(pop1_starting_N/2)

    beta = 1  # the epistasis parameter beta

    # specify ecological resource competitive abilities for two resources A and B
    # ecolDiff = 1.0 # this is "E" in the paper
    competAbility_useResourceA_species0 = (1 + ecolDiff)/2    # equals 1 when ecolDiff = 1
    competAbility_useResourceB_species0 = 1 - competAbility_useResourceA_species0
    competAbility_useResourceA_species1 = (1 - ecolDiff)/2   # equals 0 when ecolDiff = 1
    competAbility_useResourceB_species1 = 1 - competAbility_useResourceA_species1

    # set up array of strings to record outcomes
    outcome_array = Array{String, 3}(undef, length(w_hyb_set), length(S_AM_set), length(replications))

    for k in 1:length(replications)  # loop through the replicate runs
        replicate_ID = replications[k]

        run_set_name = string(set_name,"_rep", replicate_ID)

        # Loop through the different simulation sets
        Threads.@threads for i in 1:length(w_hyb_set)
            for j in 1:length(S_AM_set)
                w_hyb = w_hyb_set[i]
                S_AM = S_AM_set[j]
                println("w_hyb = ",w_hyb,"; S_AM = ",S_AM)

                run_name = string("HZAM_animation_run",run_set_name,"_surv",short_survFitnessMethod,"_ecolDiff",ecolDiff,"_growthrate",intrinsic_R,"_K",K_total,"_FL",total_functional_loci,"_NL",num_neutral_loci,"_gen",max_generations,"_SC",per_reject_cost,"_Whyb",w_hyb,"_SAM",S_AM)

                # convert S_AM to pref_ratio (for use in math below)
                if S_AM == 1
                    S_AM_for_math = 1 + 10^(-15)
                elseif S_AM == Inf
                    S_AM_for_math = 10^(15)
                else
                    S_AM_for_math = S_AM
                end
                pref_SD = sqrt( -1 / (2 * log(1/S_AM_for_math)))  # width of female acceptance curve for male trait

                # set up initial values for one simulation
                extinction = false
                outcome = []
                final_distribution = []
                # Generate genotype array for population of females:
                # this is a 3D array, where rows (D1) are alleles (row 1 from mother, row 2 from father),
                # columns (D2) are loci, and pages (D3) are individuals
                genotypes_F = generate_genotype_array(pop0_starting_N_half, pop1_starting_N_half, total_loci)
                genotypes_M = generate_genotype_array(pop0_starting_N_half, pop1_starting_N_half, total_loci)
                # functional loci are first, followed by neutral loci

                for generation in 1:max_generations

                    # Prepare for mating and reproduction
                    N_F = size(genotypes_F, 3)
                    N_M = size(genotypes_M, 3)

                    # println("generation: ", generation, "; individuals: ", N_F + N_M)

                    # calculate mating trait values (T) from genotypes
                    female_mating_traits = calc_traits_additive(genotypes_F[:,female_mating_trait_loci,:])
                    male_mating_traits = calc_traits_additive(genotypes_M[:,male_mating_trait_loci,:])

                    # calculate ecological competition trait values from genotypes
                    competition_traits_F = calc_traits_additive(genotypes_F[:,competition_trait_loci,:])
                    competition_traits_M = calc_traits_additive(genotypes_M[:,competition_trait_loci,:])

                    # calculate individual contributions to resource use, according to linear gradient between use of species 0 and species 1
                    ind_useResourceA_F = competAbility_useResourceA_species1 .+ ((1 .- competition_traits_F) .* (competAbility_useResourceA_species0 - competAbility_useResourceA_species1))
                    ind_useResourceB_F = competAbility_useResourceB_species0 .+ (competition_traits_F .* (competAbility_useResourceB_species1 - competAbility_useResourceB_species0))
                    ind_useResourceA_M = competAbility_useResourceA_species1 .+ ((1 .- competition_traits_M) .* (competAbility_useResourceA_species0 - competAbility_useResourceA_species1))
                    ind_useResourceB_M = competAbility_useResourceB_species0 .+ (competition_traits_M .* (competAbility_useResourceB_species1 - competAbility_useResourceB_species0))
                    # sum up the resource use over all individuals
                    total_useResourceA = sum(ind_useResourceA_F) + sum(ind_useResourceA_M)
                    total_useResourceB = sum(ind_useResourceB_F) + sum(ind_useResourceB_M)
                    # calculate growth rates due to each resource (according to discrete time logistic growth equation)
                    growth_rate_ResourceA = (intrinsic_R * K_A) / (K_A + ((total_useResourceA)*(intrinsic_R - 1)))
                    growth_rate_ResourceB = (intrinsic_R * K_B) / (K_B + ((total_useResourceB)*(intrinsic_R - 1)))

                    # Set up structure to record number of matings per male (and female, although almost always 1 for females),
                    # to determine sexual selection due to HI class:
                    matings_per_male = zeros(Int8, N_M)
                    matings_per_female = zeros(Int8, N_F)

                    # create empty data structures for keeping track of numbers offspring of parents
                    daughters_per_mother = zeros(Int16, N_F)
                    sons_per_mother = zeros(Int16, N_F)
                    daughters_per_father = zeros(Int16, N_M)
                    sons_per_father = zeros(Int16, N_M)

                    # make empty arrays for storing genotypes of daughters and sons
                    genotypes_daughters = Array{Int8, 3}(undef, 2, total_loci, 0)
                    genotypes_sons = Array{Int8, 3}(undef, 2, total_loci, 0)

                    # create structures for recording indices of mother and father (for error checking, and potentially for tracking genealogies)
                    # first column for mother index (3rd dim of genotypes_F) and second column for father index (3rd dim of genotypes_M)
                    daughter_parent_IDs = Array{Int}(undef, 0, 2)
                    son_parent_IDs = Array{Int}(undef, 0, 2)

                    # loop through mothers, mating and reproducing
                    for mother in 1:N_F
                        # initialize tracking variables
                        mate = false # becomes true when female is paired with male
                        rejects = 0 # will track number of rejected males (in cases there is cost--which there isn't in main HZAM-Sym paper)
                        father = [] # will contain the index of the male mate
                        # make vector of indices of eligible males
                        elig_M = Vector(1:N_M)
                        # determine male mate of female
                        while mate == false
                            # present female with random male, and remove him from list using "splice!"
                            focal_male = splice!(elig_M, rand(eachindex(elig_M)))
                            # compare male trait with female's trait (preference), and determine
                            # whether she accepts; note that match_strength is determined by a
                            # Gaussian, with a maximum of 1 and minimum of zero.
                            match_strength = (exp(1) ^ ((-(male_mating_traits[focal_male] - female_mating_traits[mother])^2) / (2 * (pref_SD ^2))))
                            if rand() < match_strength
                                # she accepts male, and mates
                                father = focal_male
                                matings_per_male[focal_male] += 1 # this adds 1 to the matings for that male
                                matings_per_female[mother] += 1
                                mate = true
                            else
                                # she rejects male
                                rejects += 1
                                if length(elig_M) == 0
                                    break
                                end
                            end
                        end
                        # determine fitness cost due to mate search (number of rejected males)
                        search_fitness = (1-per_reject_cost) ^ rejects    # (in most of HZAM-Sym paper, per_reject_cost = 0)

                        # determine fitness due to female use of available resources
                        growth_rate_of_focal_female = (ind_useResourceA_F[mother] * growth_rate_ResourceA) + (ind_useResourceB_F[mother] * growth_rate_ResourceB)
                        #combine for total fitness:
                        reproductive_fitness = 2 * growth_rate_of_focal_female * search_fitness  # the 2 is because only females, not males, produce offspring
                        # now draw the number of offspring from a poisson distribution with a mean of reproductive_fitness
                        if isempty(father)
                            offspring = 0  # if no mate (because all males were rejected), then no offspring
                        else
                            offspring = rand(Poisson(reproductive_fitness))
                        end
                        # if offspring, generate their genotypes and sexes
                        if offspring >= 1
                            for kid in 1:offspring
                                kid_info = Array{Int8, 2}(undef, 2, total_loci) # place to store genotype of one offspring
                                # generate genotypes; for each locus (column), first row for allele from mother, second row for allele from father
                                for locus in 1:total_loci
                                    kid_info[1,locus] = genotypes_F[rand([1 2]), locus, mother]  # for this locus, pick a random allele from the mother
                                    kid_info[2,locus] = genotypes_M[rand([1 2]), locus, father]
                                end
                                # determine sex of kid and add to table
                                if rand() > 0.5 # kid is daughter
                                    genotypes_daughters = cat(genotypes_daughters, kid_info, dims=3)
                                    daughter_parent_IDs = cat(daughter_parent_IDs, [mother father], dims=1)
                                    daughters_per_mother[mother] += 1
                                    daughters_per_father[father] += 1
                                else # kid is son
                                    genotypes_sons = cat(genotypes_sons, kid_info, dims=3)
                                    son_parent_IDs = cat(son_parent_IDs, [mother father], dims=1)
                                    sons_per_mother[mother] += 1
                                    sons_per_father[father] += 1
                                end
                            end
                        end
                    end # of loop through mothers

                    # check if either no daughters or no sons, and end the simulation if so
                    if (size(genotypes_daughters, 3) == 0) || (size(genotypes_sons, 3) == 0)
                        extinction = true # record an extinction of whole population (both "species")
                        break # break out of current loop (this simulation)
                    end

                    # For someday: add in here the option of tracking fitness?

                    # determine survival fitnesses of daughters due to epistasis

                    survival_fitness_daughters = get_survival_fitnesses(genotypes_daughters[:,hybrid_survival_loci,:], w_hyb)
                    daughters_survive = survival_fitness_daughters .> rand(length(survival_fitness_daughters))
                    # same for sons:
                    survival_fitness_sons = get_survival_fitnesses(genotypes_sons[:,hybrid_survival_loci,:], w_hyb)
                    sons_survive = survival_fitness_sons .> rand(length(survival_fitness_sons))

                    # check if either no surviving daughters or no surviving sons, and end the simulation if so
                    if (sum(daughters_survive) == 0) || (sum(sons_survive) == 0)
                        extinction = true # record an extinction of whole population (both "species")
                        break # break out of current loop (this simulation)
                    end

                    # assign surviving offspring to new adult population
                    genotypes_F = genotypes_daughters[:,:,daughters_survive]
                    genotypes_M = genotypes_sons[:,:,sons_survive]

                end # of loop through generations

                # Record results of the one simulation
                functional_HI_all_inds = []
                species0_proportion = []
                species1_proportion = []
                HI_NL_all_inds = []
                species0_proportion_NL = []
                species1_proportion_NL = []
                if extinction  # whole simulation went extinct
                    outcome = "extinction"
                    if save_each_sim
                        @save string("simulation_data.",run_name,".jld2") outcome
                    end
                else  # no complete extinction
                    # use trait loci to calculate HI of each individual
                    functional_HI_all_inds = [calc_traits_additive(genotypes_F[:,functional_loci_range,:]); calc_traits_additive(genotypes_M[:,functional_loci_range,:])]
                    # calculate proportion of all individuals who are species0 or species1 (defined as low and high 10% of HI distribution, respectively)
                    species0_proportion = sum(functional_HI_all_inds .< 0.1) / length(functional_HI_all_inds)
                    species1_proportion = sum(functional_HI_all_inds .> 0.9) / length(functional_HI_all_inds)
                    if species0_proportion >= 0.85 || species1_proportion >= 0.85
                        outcome = "one_species"
                    elseif (species0_proportion + species1_proportion >= 0.85) && (species0_proportion >= 0.15) && (species1_proportion >= 0.15)
                        outcome = "two_species"
                    else
                        outcome = "blended"
                    end
                    HI_NL_all_inds = [calc_traits_additive(genotypes_F[:,neutral_loci,:]); calc_traits_additive(genotypes_M[:,neutral_loci,:])]
                    species0_proportion_NL = sum(HI_NL_all_inds .== 0) / length(HI_NL_all_inds)
                    species1_proportion_NL = sum(HI_NL_all_inds .== 1) / length(HI_NL_all_inds)
                    if save_each_sim
                        @save string("HZAM_Sym_Julia_results_GitIgnore/simulation_data.",run_name,".jld2") outcome functional_HI_all_inds species0_proportion species1_proportion HI_NL_all_inds species0_proportion_NL species1_proportion_NL
                    end
                end
                println(run_name, "  outcome was: ", outcome)
                outcome_array[i,j,k] = outcome
            end # of S_AM loop
        end # of w_hyb loop
    end # of replicate loop

    if save_outcomes_JL
        filename = string("HZAM_Sym_Julia_results_GitIgnore/outcomeArray_set",set_name,"_surv",short_survFitnessMethod,"_ecolDiff",ecolDiff,"_growthrate",intrinsic_R,"_K",K_total,"_FL",total_functional_loci,"_NL",num_neutral_loci,"_gen",max_generations,"_SC",per_reject_cost,".jld2")
        save_object(filename, outcome_array)
    end

    if save_outcomes_csv
        for i in 1:size(outcome_array, 3)
            filename = string("HZAM_Sym_Julia_results_GitIgnore/outcomeArray_set",set_name,"_surv",short_survFitnessMethod,"_ecolDiff",ecolDiff,"_growthrate",intrinsic_R,"_K",K_total,"_FL",total_functional_loci,"_NL",num_neutral_loci,"_gen",max_generations,"_SC",per_reject_cost,"_rep",replications[i])
            CSV.write(filename, Tables.table(outcome_array[:,:,i]), writeheader=false)
        end
    end
    return outcome_array
end

#### functions for summarizing and plotting results

# convert outcome array (an array of strings) to categorical array
function convert_to_cat_array(outcome_array)
    cat_outcome_array = compress(CategoricalArray(outcome_array))
    levels!(cat_outcome_array, ["extinction", "blended", "one_species", "two_species"])
    return cat_outcome_array
end

# function for plotting grid of pie charts showing distribution of four outcomes, using categorical array as input
function plot_all_outcomes(cat_outcome_array)
    num_outcome_types = length(levels(cat_outcome_array))
    outcome_counts = Array{Int16, 2}(undef, num_outcome_types, (size(cat_outcome_array, 1)*size(cat_outcome_array, 2) ))
    for i in 1:size(cat_outcome_array, 1)
        for j in 1:size(cat_outcome_array, 2)
            for outcome_num in 1:num_outcome_types
                outcome_counts[outcome_num, j + (i-1)*size(cat_outcome_array, 2)] = sum(cat_outcome_array[i,j,:] .== levels(cat_outcome_array)[outcome_num])
            end
        end
    end
    colors_of_outcomes = [RGB(0,0,0), plasma[0.525], plasma[0.2], plasma[0.9]] # colors for 4 outcome categories
    pie(outcome_counts, layout = grid(size(cat_outcome_array, 1), size(cat_outcome_array, 2)), legend = false, palette = colors_of_outcomes, margin = -2.0mm)
    plot!(size=(800,1300))
end

function get_most_common_outcomes(cat_outcome_array)
    levels_of_outcomes = levels(cat_outcome_array)
    num_outcome_types = length(levels_of_outcomes)
    most_common_outcomes = CategoricalArray{String, 2}(undef, size(cat_outcome_array, 1), size(cat_outcome_array, 2))
    levels!(most_common_outcomes, ["extinction", "blended", "one_species", "two_species"])
    for i in 1:size(cat_outcome_array, 1)
        for j in 1:size(cat_outcome_array, 2)
            outcome_counts = Vector{Int}(undef, num_outcome_types)
            for outcome_num in 1:num_outcome_types
                outcome_counts[outcome_num] = sum(cat_outcome_array[i,j,:] .== levels(cat_outcome_array)[outcome_num])
            end
            outcomes_with_max_count = findall(outcome_counts .== maximum(outcome_counts))
            if length(outcomes_with_max_count) == 1
                most_common_outcomes[i,j] = levels_of_outcomes[outcomes_with_max_count[1]]
            elseif length(outcomes_with_max_count) >= 2
                # if tie in outcome count, choose one randomly
                most_common_outcomes[i,j] = levels_of_outcomes[sample(outcomes_with_max_count)]
            end
        end
    end
    return most_common_outcomes
end

function plot_common_outcomes(common_outcome_array)
    # make heat map of outcomes
    colors_of_outcomes = [RGB(0,0,0), plasma[0.525], plasma[0.2], plasma[0.9]] # colors for 4 outcome categories
    one_outcome_array = reverse(common_outcome_array, dims=1)
    x_midpoints = log10.([1, 3, 10, 30, 100, 300, 1000, 5000])  # the S_AM values, with Inf convert to 3000 for graphing
    w_hyb_set = [1, 0.98, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0] # for just one run, just put one number in this and next line
    y_midpoints = reverse(w_hyb_set)
    min_color = minimum(levelcode.(one_outcome_array)) # this and next line needed to choose the proper colors for the figure
    max_color = maximum(levelcode.(one_outcome_array))
    heatmap(x_midpoints, y_midpoints, one_outcome_array, c = colors_of_outcomes[min_color:max_color], yflip = false, tick_direction = :out, colorbar = false, size = (440,310), framestyle = :box)
    xaxis!("Strength of conspecific mate preference")
    xticklabels = ["1", "3", "10", "30", "100", "300", "1000", "complete"]
    plot!(xticks=(x_midpoints, xticklabels))
    yaxis!("Hybrid fitness")
    yticklabels = ["0.0", "0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "", "", "1.0"]
    plot!(yticks=(y_midpoints, yticklabels), tick_direction = :out)
    # add white lines
    plot!([xlims()[1], xlims()[2]], [0.05, 0.05], linecolor = :white, widen = false, legend = false, linewidth=3)
    x_for_line = mean(x_midpoints[[length(x_midpoints)-1 length(x_midpoints)]])
    plot!([x_for_line, x_for_line], [ylims()[1], ylims()[2]], linecolor = :white, widen = false, legend = false, linewidth=3)
end

function make_and_save_figs(ResultsFolder, RunName, RunOutcomes)
    cat_RunOutcomes = convert_to_cat_array(RunOutcomes)
    display(plot_all_outcomes(cat_RunOutcomes))
    savefig(string(ResultsFolder,"/",RunName,"_AllOutcomes.png"))
    savefig(string(ResultsFolder,"/",RunName,"_AllOutcomes.pdf"))
    most_common_outcomes = get_most_common_outcomes(cat_RunOutcomes)
    display(plot_common_outcomes(most_common_outcomes))
    savefig(string(ResultsFolder,"/",RunName,"_MostCommonOutcomes.png"))
    savefig(string(ResultsFolder,"/",RunName,"_MostCommonOutcomes.pdf"))
end

#### Run the actual simulation by calling the above functions:

# If wanting to run a test to make sure things are working (after above functions are loaded into memory):
# RunOutcomes = run_HZAM("TEST", 1.0, 1.05, 1:3, 100, 50)

ResultsFolder = "HZAM_Sym_Julia_results_GitIgnore" # or change to another folder where results will be saved

# To run simulations under the parameters in the paper, you can execute the following,
# under the heading for each Figure panel in the paper. Because these are stochastic simulations,
# results will differ a little bit (but not much) from those shown in the paper.
# Also, the "make_and_save_figs" function has saved the figure panels under each set of conditions run below,
# and the main script saved the Julia object where the results are stored--this starts with "outcomeArray_set"
# and will be available on Dryad.

# (Note that if the below have "Het" in their name, they were done using the heterozygote disadvantage (underdominance)
# option for survival fitness. If they do not have "Het" in their name, they were done using the epistasis option.)

# Figure 3A
RunName = "JL_fig3aHet"
RunOutcomes = run_HZAM(RunName, 0.0, 1.05, 1:25;
    survival_fitness_method = "hetdisadvantage")
# took about 5 hours to run, with 4 threads on a Mac desktop with 4 cores
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 3B
RunName = "JL_fig3bHet"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 4A
RunName = "JL_Fig4aHet"
RunOutcomes = run_HZAM(RunName, 1.0, 1.025, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 4B
RunName = "JL_Fig4bHet"
RunOutcomes = run_HZAM(RunName, 1.0, 1.2, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 4C
RunName = "JL_Fig4cHet"
RunOutcomes = run_HZAM(RunName, 1.0, 2.6, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 6 row 1 panel 2
RunName = "JL_Fig6(1,2)Het"
RunOutcomes = run_HZAM(RunName, 0.25, 1.05, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 6 row 1 panel 3
RunName = "JL_Fig6(1,3)Het"
RunOutcomes = run_HZAM(RunName, 0.5, 1.05, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 6 row 1 panel 4
RunName = "JL_Fig6(1,4)Het"
RunOutcomes = run_HZAM(RunName, 0.75, 1.05, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 6 row 2 panel 1
RunName = "JL_Fig6(2,1)Het"
RunOutcomes = run_HZAM(RunName, 0.0, 2.6, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 6 row 2 panel 2
RunName = "JL_Fig6(2,2)Het"
RunOutcomes = run_HZAM(RunName, 0.25, 2.6, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 6 row 2 panel 3
RunName = "JL_Fig6(2,3)Het"
RunOutcomes = run_HZAM(RunName, 0.5, 2.6, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure 6 row 2 panel 4
RunName = "JL_Fig6(2,4)Het"
RunOutcomes = run_HZAM(RunName, 0.75, 2.6, 1:25;
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S2A (epistasis):
RunName = "JL_fig3a"
RunOutcomes = run_HZAM(RunName, 0.0, 1.05, 1:25)
# this run (which does 2600 simulations) takes about 4 hours (with 4 threads on a somewhat powerful desktop Mac)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S2B (epistasis):
RunName = "JL_fig3b"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S7A
RunName = "JL_fig3bHet125gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    max_generations = 125, survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S7B
RunName = "JL_fig3bHet250gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    max_generations = 250, survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S7C
RunName = "JL_fig3bHet500gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    max_generations = 500, survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S7E
RunName = "JL_fig3bHet2000gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    max_generations = 2000, survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S7F
RunName = "JL_fig3bHet4000gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    max_generations = 4000, survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S8A
RunName = "JL_LikeFig3bHet_butFL1"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    K_total = 1000, max_generations = 1000,
    total_loci = 2, female_mating_trait_loci = 1:1, male_mating_trait_loci = 1:1,
    competition_trait_loci = 1:1, hybrid_survival_loci = 1:1, neutral_loci = 2:2,
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S8C
RunName = "JL_LikeFig3bHet_butFL9"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    K_total = 1000, max_generations = 1000,
    total_loci = 18, female_mating_trait_loci = 1:9, male_mating_trait_loci = 1:9,
    competition_trait_loci = 1:9, hybrid_survival_loci = 1:9, neutral_loci = 10:18,
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S8D
RunName = "JL_LikeFig3bHet_butFL27"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    K_total = 1000, max_generations = 1000,
    total_loci = 54, female_mating_trait_loci = 1:27, male_mating_trait_loci = 1:27,
    competition_trait_loci = 1:27, hybrid_survival_loci = 1:27, neutral_loci = 28:54,
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S9A
RunName = "JL_fig3bHetDiffLociTraitPref"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    K_total = 1000, max_generations = 1000,
    total_loci = 9, female_mating_trait_loci = 1:3, male_mating_trait_loci = 4:6,
    competition_trait_loci = 1:6, hybrid_survival_loci = 1:6, neutral_loci = 7:9,
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S9B
RunName = "JL_fig3bHetDiffLociTraitPrefEcol"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    K_total = 1000, max_generations = 1000,
    total_loci = 12, female_mating_trait_loci = 1:3, male_mating_trait_loci = 4:6,
    competition_trait_loci = 7:9, hybrid_survival_loci = 7:9, neutral_loci = 10:12,
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S10C
RunName = "JL_LikeFig3bHet_but_starting_pop_ratio_0.5"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    survival_fitness_method = "hetdisadvantage",
    starting_pop_ratio = 0.5)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S10D
RunName = "JL_LikeFig4bHet_but_starting_pop_ratio_0.5"
RunOutcomes = run_HZAM(RunName, 1.0, 1.2, 1:25;
    survival_fitness_method = "hetdisadvantage",
    starting_pop_ratio = 0.5)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S10E
RunName = "JL_LikeFig3bHet_but_starting_pop_ratio_0.1"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25;
    survival_fitness_method = "hetdisadvantage",
    starting_pop_ratio = 0.1)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Figure S10F
RunName = "JL_LikeFig4bHet_but_starting_pop_ratio_0.1"
RunOutcomes = run_HZAM(RunName, 1.0, 1.2, 1:25;
    survival_fitness_method = "hetdisadvantage",
    starting_pop_ratio = 0.1)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results):
RunName = "JL_LikeFig3bHet_butFL2"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    K_total = 1000, max_generations = 1000,
    total_loci = 4, female_mating_trait_loci = 1:2, male_mating_trait_loci = 1:2,
    competition_trait_loci = 1:2, hybrid_survival_loci = 1:2, neutral_loci = 3:4,
    survival_fitness_method = "hetdisadvantage")
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_LikeFig3b_butFL1"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000, 1000,
    2, 1:1, 1:1,
    1:1, 1:1, 2:2)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_LikeFig3b_butFL9"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000, 1000,
    18, 1:9, 1:9,
    1:9, 1:9, 10:18)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_LikeFig3b_butFL27"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000, 1000,
    54, 1:27, 1:27,
    1:27, 1:27, 28:54)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_LikeFig4c_butFL27"
RunOutcomes = run_HZAM(RunName, 1.0, 2.6, 1:25,
    1000, 1000,
    54, 1:27, 1:27,
    1:27, 1:27, 28:54)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig4a"
RunOutcomes = run_HZAM(RunName, 1.0, 1.025, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig4b"
RunOutcomes = run_HZAM(RunName, 1.0, 1.2, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig4c"
RunOutcomes = run_HZAM(RunName, 1.0, 2.6, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig6(1,2)"
RunOutcomes = run_HZAM(RunName, 0.25, 1.05, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig6(1,3)"
RunOutcomes = run_HZAM(RunName, 0.5, 1.05, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig6(1,4)"
RunOutcomes = run_HZAM(RunName, 0.75, 1.05, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig6(2,1)"
RunOutcomes = run_HZAM(RunName, 0.0, 2.6, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig6(2,2)"
RunOutcomes = run_HZAM(RunName, 0.25, 2.6, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig6(2,3)"
RunOutcomes = run_HZAM(RunName, 0.5, 2.6, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_Fig6(2,4)"
RunOutcomes = run_HZAM(RunName, 0.75, 2.6, 1:25)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_fig3b125gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000, 125)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_fig3b250gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000, 250)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_fig3b500gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000,500)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_fig3b2000gen"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000,2000)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_fig3bDiffLociTraitPref"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000, 1000,
    9, 1:3, 4:6,
    1:6, 1:6, 7:9)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# Not shown in paper (but available in saved results), with epistasis:
RunName = "JL_fig3bDiffLociTraitPrefEcol"
RunOutcomes = run_HZAM(RunName, 1.0, 1.05, 1:25,
    1000, 1000,
    12, 1:3, 4:6,
    7:9, 7:9, 10:12)
make_and_save_figs(ResultsFolder, RunName, RunOutcomes)

# make graph of epistasis fitness (for Figure 1Cii):
survival_HI = collect(0.0:0.01:1.0)
beta = 1
w_hyb = 0.8
epistasis_fitnesses = 1 .- (1 - w_hyb) .* (4 .* survival_HI .* (1 .- survival_HI)).^beta
using Plots
display(plot(survival_HI, epistasis_fitnesses, ylims = (0,1)))
savefig("plot.pdf")


text:
using Polynomials: Poly, polyder, polyint

struct LegendreBasisHierarchical{T, N, fT, pT, dT, iT} <: PolynomialBasis{T,N}
    x::Vector{T}
    factors::fT
    polys::pT
    derivs::dT
    ints::iT
end

function LegendreBasisHierarchical(T,N)

    if N < 2
        error("Number of Nodes of Hierarchical Legendre Basis must be at least 2 (Degree at least 1).")
    end

    x       = zeros(T,N)
    factors = Tuple{}()
    polys   = Tuple{}()
    derivs  = Tuple{}()
    ints    = Tuple{}()

    for p in 0:N-1
        tfac, tpol = legendre_polynomial_hierarchical(p, T)
        tder = polyder(tpol)
        tint = polyint(tpol)

        if p == 0
            x[p+1] = zero(T)
        elseif p == 1
            x[p+1] = one(T)
        else
            x[p+1] = T(0.5)
        end

        factors = (factors..., tfac)
        polys   = (polys..., tpol)
        derivs  = (derivs..., tder)
        ints    = (ints..., tint)
    end

    LegendreBasisHierarchical{T, N, typeof(factors), typeof(polys), typeof(derivs), typeof(ints)}(x, factors, polys, derivs, ints)
end

nodes(b::LegendreBasisHierarchical{T,N})  where {T,N} = b.x

Base.hash(b::LegendreBasisHierarchical{T,N}, h::UInt) where {T,N} = hash(T, hash(N, h))

Base.:(==)(b1::LegendreBasisHierarchical{T1,N1}, b2::LegendreBasisHierarchical{T2,N2}) where {T1,N1,T2,N2} = (T1 == T2 && N1 == N2)

Base.isequal(b1::LegendreBasisHierarchical{T1,N1}, b2::LegendreBasisHierarchical{T2,N2}) where {T1,N1,T2,N2} = (T1 == T2 && N1 == N2)

function legendre_polynomial_hierarchical(p, T=Float64)
    if p == 0
        (one(T), Poly([+1, -1]))
    elseif p == 1
        (one(T), Poly([ 0, +1]))
    else
        l1 = legendre_polynomial(p,   T)
        l2 = legendre_polynomial(p-2, T)

        l1p = l1[2]
        l2p = l2[2]

        (1/sqrt(2 * (2*p-1)), l1p - l2p)
    end
end

function eval_basis(b::LegendreBasisHierarchical{T,N}, j::Int, x::T) where {T,N}
    return b.factors[j] * b.polys[j](x)
end

function deriv_basis(b::LegendreBasisHierarchical{T,N}, j::Int, x::T) where {T,N}
    return b.factors[j] * b.derivs[j](x)
end

function int_basis(b::LegendreBasisHierarchical{T,N}, j::Int, x::T) where {T,N}
    return b.factors[j] * b.ints[j](x)
end


text:
using JuMP, AmplNLWriter, MathOptInterfaceToMathProgBase, Test, Ipopt
using Random, Statistics, Test

const MOI = JuMP.MOI

model = JuMP.Model(
    with_optimizer(
        MathOptInterfaceToMathProgBase.Optimizer,
        solver = AmplNLWriter.AmplNLSolver(Ipopt.amplexe)
    )
)

N = 1000
h = 1/N
alpha = 350

@variables(model, begin
    -1 <= t[1:(N + 1)] <= 1
    -0.05 <= x[1:(N + 1)] <= 0.05
    u[1:(N + 1)]
end)
@NLobjective(model, Min, sum(0.5 * h * (u[i + 1]^2 + u[i]^2) +
    0.5 * alpha * h * (cos(t[i + 1]) + cos(t[i])) for i in 1:N)
)
@NLconstraint(model, [i = 1:N],
    x[i + 1] - x[i] - 0.5 * h * (sin(t[i + 1]) + sin(t[i])) == 0
)
@constraint(model, [i = 1:N],
    t[i + 1] - t[i] - 0.5 * h * u[i + 1] - 0.5 * h * u[i] == 0
)
JuMP.optimize!(model)

@test JuMP.termination_status(model) == MOI.OPTIMAL
@test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
@test JuMP.objective_value(model) ≈ 350.0

function example_mle(; verbose = true)
    n = 1_000
    Random.seed!(1234)
    data = randn(n)

    model = JuMP.Model(
        with_optimizer(
            MathOptInterfaceToMathProgBase.Optimizer,
            solver = AmplNLWriter.AmplNLSolver(Ipopt.amplexe)
        )
    )

    @variable(model, μ, start = 0.0)
    @variable(model, σ >= 0.0, start = 1.0)
    @NLobjective(model, Max, n / 2 * log(1 / (2 * π * σ^2)) -
        sum((data[i] - μ)^2 for i in 1:n) / (2 * σ^2)
    )
    JuMP.optimize!(model)
    if verbose
        println("μ = ", JuMP.value(μ))
        println("mean(data) = ", mean(data))
        println("σ^2 = ", JuMP.value(σ)^2)
        println("var(data) = ", var(data))
        println("MLE objective: ", JuMP.objective_value(m))
    end
    @test JuMP.value(μ) ≈ mean(data) atol = 1e-3
    @test JuMP.value(σ)^2 ≈ var(data) atol = 1e-2

    # constrained MLE?
    @NLconstraint(model, μ == σ^2)

    JuMP.optimize!(model)
    if verbose
        println("\nWith constraint μ == σ^2:")
        println("μ = ", JuMP.value(μ))
        println("σ^2 = ", JuMP.value(σ)^2)
        println("Constrained MLE objective: ", JuMP.objective_value(model))
    end
    @test JuMP.value(μ) ≈ JuMP.value(σ)^2
end

example_mle(verbose = false)


text:
export read_ecfp, read_sparse, read_rowcol
export read_binary_int32, filter_rare, write_binary_int32
export write_binary_matrix
export read_binary_float32
export read_sparse_float32, write_sparse_float32
export read_sparse_float64, write_sparse_float64
export read_sparse_binary_matrix, write_sparse_binary_matrix
export read_matrix_market

function read_ecfp(filename)
    i = 0
    next_fp = 1
    fp = Dict{Int32, Int32}()
    rows = Int32[]
    cols = Int32[]
    open(filename) do f
        for line in eachline(f)
            i += 1
            a = split(line, ",")
            for j = 2:length(a)
                fp_raw = parse(Int, a[j])
                local fp_id::Int32
                ## fetch fingerprint id range [1,...,max_fp]
                if haskey(fp, fp_raw)
                    fp_id = fp[fp_raw]
                else
                    fp_id      = next_fp
                    fp[fp_raw] = fp_id
                    next_fp += 1
                end
                push!(rows, i)
                push!(cols, fp_id)
            end
        end
    end
    println("Number of lines: $i")
    return rows, cols, fp
end

function read_rowcol(filename)
    rows = Int32[]
    cols = Int32[]
    open(filename) do f
        for line in eachline(f)
            a = split(line, ",")
            push!( rows, parse(Int32, a[1]) )
            push!( cols, parse(Int32, a[2]) )
        end
    end
    return rows, cols
end

function read_binary_int32(filename)
    open(filename) do f
        nrows = read(f, Int64)
        ncols = read(f, Int64)
        return read(f, Int32, (nrows, ncols))
    end
end

function read_binary_float32(filename)
    open(filename) do f
        nrows = read(f, Int64)
        ncols = read(f, Int64)
        return read(f, Float32, (nrows, ncols))
    end
end

function read_sparse_float32(filename)
  open(filename) do f
    nnz = read(f, Int64)
    rows = read(f, Int32, nnz)
    cols = read(f, Int32, nnz)
    vals = read(f, Float32, nnz)
    return rows, cols, vals
  end
end

function read_sparse(filename)
    rc = read_rowcol(filename)
    return sparse(rc[1], rc[2], 1f0)
end

function filter_rare(X::SparseMatrixCSC, nmin)
    featn = vec(sum(X, 1))
    return X[:, featn .>= nmin]
end

function write_binary_int32(filename, X::Matrix{Int32})
    write_binary_matrix(filename, X)
end

function write_binary_matrix(filename, X)
  open(filename, "w") do f
    write(f, size(X, 1))
    write(f, size(X, 2))
    write(f, X)
  end
end

function write_sparse_float32(filename, X::SparseMatrixCSC)
  I, J, V = findnz(X)
  write_sparse_float32(
    filename,
    convert(Vector{Int32}, I),
    convert(Vector{Int32}, J),
    convert(Vector{Float32}, V))
  nothing
end

function write_sparse_float32(filename, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float32})
  open(filename, "w") do f
    write(f, length(rows))
    write(f, rows)
    write(f, cols)
    write(f, values)
  end
  nothing
end

## writes the non-zero coordinates of X
function write_sparse_binary_matrix(filename, X::SparseMatrixCSC)
  rows, cols, V = findnz(X)
  open(filename, "w") do f
    write(f, size(X, 1))  ## nrows
    write(f, size(X, 2))  ## ncols
    write(f, nnz(X))      ## nnz
    write(f, convert(Vector{Int32}, rows))
    write(f, convert(Vector{Int32}, cols))
  end
  nothing
end

function read_sparse_binary_matrix(filename)
  open(filename) do f
    nrows = read(f, Int64)
    ncols = read(f, Int64)
    nnz   = read(f, Int64)
    rows  = read(f, Int32, nnz)
    cols  = read(f, Int32, nnz)
    return sparse(rows, cols, 1, nrows, ncols)
  end
end

function read_matrix_market(filename)
  nrows = 0
  ncols = 0
  nnz   = 0
  local rows::Vector{Int32}
  local cols::Vector{Int32}
  local vals::Vector{Float64}
  open(filename) do f
    ## reading the first line
    while !eof(f)
      ln = readline(f)
      ln[1] == '%' && continue
      arr = split(ln, [' ','\t'])
      nrows = parse(Int, arr[1])
      ncols = parse(Int, arr[2])
      nnz   = parse(Int, arr[3])
      break
    end
    rows = zeros(Int32, nnz)
    cols = zeros(Int32, nnz)
    vals = zeros(nnz)
    i = 1
    ## reading the rest
    while !eof(f)
      ln = readline(f)
      ln[1] == '%' && continue
      arr = split(ln, [' ','\t'])

      rows[i] = parse(Int32, arr[1])
      cols[i] = parse(Int32, arr[2])
      vals[i] = parse(Float64, arr[3])
      i += 1
    end
  end
  return sparse(rows, cols, vals, nrows, ncols)
end

function write_matrix_market(filename, X::DataFrame)
  nrows = maximum(X[:, 1])
  ncols = maximum(X[:, 2])
  nnz   = size(X, 1)
  ## writing the header line
  open(filename, "w") do f
    write(f, "%%MatrixMarket matrix coordinate real general\n")
    write(f, @sprintf("%d\t%d\t%d\n", nrows, ncols, nnz) )
    writedlm(f, convert(Matrix, X[:,1:3]))
  end
  nothing
end

function write_sparse_float64(filename, X)
  open(filename, "w") do f
    nz = findnz(X)
    write(f, size(X, 1))  ## nrow
    write(f, size(X, 2))  ## ncol
    write(f, length(nz[1]))  ## nnz
    write(f, convert(Vector{Int32}, nz[1])) ## row_idx
    write(f, convert(Vector{Int32}, nz[2])) ## col_idx
    write(f, convert(Vector{Float64}, nz[3])) ## values
  end
  nothing
end

function read_sparse_float64(filename)
  open(filename) do f
    nrow = read(f, Int64)
    ncol = read(f, Int64)
    nnz  = read(f, Int64)
    rows = read(f, Int32, nnz)
    cols = read(f, Int32, nnz)
    vals = read(f, Float64, nnz)
    return sparse(rows, cols, vals, nrow, ncol)
  end
end


text:
using NMarkov
using Test
using Printf
using Distributions
using SparseMatrix
using SparseArrays

import NMarkov: rightbound, poipmf, cpoipmf, convunifstep!

@testset "PoissonRight" begin
    for q in [1.0e-6, 1.0e-7, 1.0e-8, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            @test 0 <= rightbound(lambda, q) - cquantile(Poisson(lambda), q) <= 6
        end
    end
    @time for q in [1.0e-6, 1.0e-7, 1.0e-8, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            rightbound(lambda, q)
        end
    end
    @time for q in [1.0e-6, 1.0e-7, 1.0e-8, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            cquantile(Poisson(lambda), q)
        end
    end
end

@testset "PoissonPMF" begin
    for q in [1.0e-8, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            right = rightbound(lambda, q)
            weight,p = poipmf(lambda, right)
            v1 = p / weight
            v2 = pdf.(Poisson(lambda), 0:right)
            # println(v1)
            # println(v2)
            @test v1 ≈ v2
        end
    end
    @time for q in [1.0e-8, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            right = rightbound(lambda, q)
            weight,p = poipmf(lambda, right)
        end
    end
    @time for q in [1.0e-8, 1.0e-9, 1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            right = rightbound(lambda, q)
            v2 = pdf.(Poisson(lambda), 0:right)
        end
    end
end

@testset "PoissonCCDF" begin
    for q in [1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            right = rightbound(lambda, q)
            weight,p,cp = cpoipmf(lambda, right)
            v1 = cp / weight
            v2 = ccdf.(Poisson(lambda), 0:right)
            # println(v1)
            # println(v2)
            @test v1 ≈ v2
        end
    end
    @time for q in [1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            right = rightbound(lambda, q)
            weight,p,cp = cpoipmf(lambda, right)
            v1 = cp / weight
        end
    end
    @time for q in [1.0e-10, 1.0e-11, 1.0e-12]
        for lambda in [0.01, 0.1, 1.0, 10.0, 100.0, 10000.0, 100000.0, 1000000.0]
            right = rightbound(lambda, q)
            v2 = ccdf.(Poisson(lambda), 0:right)
        end
    end
end

@testset "GTH" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    @time x = gth(Q)
    Q1 = copy(Q)
    Q1[:,1] .= 1
    b = zeros(size(Q1)[1])
    b[1] = 1.0
    @time x0 = Q1' \ b
    @test x0 ≈ x
    @time x = gth(Q, [2,3,1]) # permutation
    @test x0 ≈ x
end

@testset "GS1" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    @time x, = stgs(SparseCSC(Q), x0 = [0.3, 0.4, 0.3])
    Q1 = copy(Q)
    Q1[:,1] .= 1
    b = zeros(size(Q1)[1])
    b[1] = 1.0
    @time x0 = Q1' \ b
    @test x0 ≈ x
    @time x, = stgs(SparseCSC(Q))
    @test x0 ≈ x
end

@testset "GS2" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    @time x, = stgs(SparseMatrixCSC(Q), x0 = [0.3, 0.4, 0.3])
    Q1 = copy(Q)
    Q1[:,1] .= 1
    b = zeros(size(Q1)[1])
    b[1] = 1.0
    @time x0 = Q1' \ b
    @test x0 ≈ x
    @time x, = stgs(SparseCSC(Q))
    @test x0 ≈ x
end

@testset "Power" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    @time x, = stpower(unif(Q)[1], x0 = [0.3, 0.4, 0.3])
    Q1 = copy(Q)
    Q1[:,1] .= 1
    b = zeros(size(Q1)[1])
    b[1] = 1.0
    @time x0 = Q1' \ b
    @test x0 ≈ x
    @time x, = stpower(unif(Q)[1])
    @test x0 ≈ x
end

@testset "STSENGS1" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    Qdash = [
        -1.0 1.0 0.0;
        0.0 0.0 0.0;
        0.0 0.0 0.0
    ]
    pis = gth(Q)
    b = Qdash' * pis
    y0 = stsengs(SparseCSC(Q), pis, b)

    Q1 = copy(Q)
    Q1[:,1] .= 1
    b[1] = 0.0
    x0 = Q1' \ (-b)

    @test x0 ≈ y0[1]
end

@testset "STSENGS2" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    Qdash = [
        -1.0 1.0 0.0;
        0.0 0.0 0.0;
        0.0 0.0 0.0
    ]
    pis = gth(Q)
    b = Qdash' * pis
    y0 = stsengs(SparseMatrixCSC(Q), pis, b)

    Q1 = copy(Q)
    Q1[:,1] .= 1
    b[1] = 0.0
    x0 = Q1' \ (-b)

    @test x0 ≈ y0[1]
end

@testset "STSEN1" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    Qdash = [
        -1.0 1.0 0.0;
        0.0 0.0 0.0;
        0.0 0.0 0.0
    ]
    pis = gth(Q)
    b = Qdash' * pis
    y0 = stsen(Q, pis, b)

    Q1 = copy(Q)
    Q1[:,1] .= 1
    b[1] = 0.0
    x0 = Q1' \ (-b)

    @test x0 ≈ y0
end

@testset "STSEN12" begin
    Q = [
        -2.0 2.0 0.0;
        1.0 -2.0 1.0;
        0.0 1.0 -1.0
    ]
    Qdash = [
        -1.0 1.0 0.0;
        0.0 0.0 0.0;
        0.0 0.0 0.0
    ]
    pis = gth(Q)
    b = Qdash' * pis
    y0 = stsen(Q, pis, b)

    Q1 = copy(Q)
    Q1[:,1] .= 1
    b[1] = 0.0
    x0 = Q1' \ (-b)

    @test x0 ≈ y0
end

@testset "STSENPower" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    Qdash = [
        -1.0 1.0 0.0;
        0.0 0.0 0.0;
        0.0 0.0 0.0
    ]
    P, qv = unif(SparseCSC(Q))
    Pdash = Qdash / qv
    pis = gth(Q)
    y0 = stsenpower(P, pis, Pdash' * pis)

    Q1 = copy(Q)
    Q1[:,1] .= 1
    b = Qdash' * pis
    b[1] = 0.0
    x0 = Q1' \ (-b)

    @test x0 ≈ y0[1]
end

@testset "QSTGS1" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    xi = [1.0, 0.0, 0.0]
    @time x = qstgs(SparseCSC(Q), xi)

    P, qv = unif(Q)
    pxi = xi / qv
    @time x2 = qstpower(P, pxi)

    @test x[1] ≈ x2[1]
    @test x[2] ≈ x2[2]*qv
end

@testset "QSTGS2" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    xi = [1.0, 0.0, 0.0]
    @time x = qstgs(SparseMatrixCSC(Q), xi)

    P, qv = unif(Q)
    pxi = xi / qv
    @time x2 = qstpower(P, pxi)

    @test x[1] ≈ x2[1]
    @test x[2] ≈ x2[2]*qv
end

@testset "mexp1" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        0.5 0.0;
        0.5 0.5;
        0.0 0.5
    ]
    t = 10.0
    res0 = exp(Q*t) * x0
    res1 = mexp(Q, x0, t)
    res2 = mexp(SparseCSR(Q), x0, t)
    res3 = mexp(SparseCSC(Q), x0, t)
    res4 = mexp(SparseCOO(Q), x0, t)
    res5 = mexp(SparseMatrixCSC(Q), x0, t)
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
    res0 = exp(Q'*t) * x0
    res1 = mexp(Q, x0, t, transpose=:T)
    res2 = mexp(SparseCSR(Q), x0, t, transpose=:T)
    res3 = mexp(SparseCSC(Q), x0, t, transpose=:T)
    res4 = mexp(SparseCOO(Q), x0, t, transpose=:T)
    res5 = mexp(SparseMatrixCSC(Q), x0, t, transpose=:T)
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
end

@testset "mexpc1" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    I = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    x0 = [
        0.5 0.0;
        0.5 0.5;
        0.0 0.5
    ]
    t = 10.0

    Qdash = [
        Q zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    res0 = exp(Qdash*t) * x0dash
    res1 = mexpc(Q, x0, t)
    res2 = mexpc(SparseCSR(Q), x0, t)
    res3 = mexpc(SparseCSC(Q), x0, t)
    res4 = mexpc(SparseCOO(Q), x0, t)
    res5 = mexpc(SparseMatrixCSC(Q), x0, t)
    @test res0[1:3,1:2] ≈ res1[1]
    @test res0[1:3,1:2] ≈ res2[1]
    @test res0[1:3,1:2] ≈ res3[1]
    @test res0[1:3,1:2] ≈ res4[1]
    @test res0[1:3,1:2] ≈ res5[1]
    @test res0[4:6,1:2] ≈ res1[2]
    @test res0[4:6,1:2] ≈ res2[2]
    @test res0[4:6,1:2] ≈ res3[2]
    @test res0[4:6,1:2] ≈ res4[2]
    @test res0[4:6,1:2] ≈ res5[2]

    Qdash = [
        Q' zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    res0 = exp(Qdash*t) * x0dash
    res1 = mexpc(Q, x0, t, transpose=:T)
    res2 = mexpc(SparseCSR(Q), x0, t, transpose=:T)
    res3 = mexpc(SparseCSC(Q), x0, t, transpose=:T)
    res4 = mexpc(SparseCOO(Q), x0, t, transpose=:T)
    res5 = mexpc(SparseMatrixCSC(Q), x0, t, transpose=:T)
    @test res0[1:3,1:2] ≈ res1[1]
    @test res0[1:3,1:2] ≈ res2[1]
    @test res0[1:3,1:2] ≈ res3[1]
    @test res0[1:3,1:2] ≈ res4[1]
    @test res0[1:3,1:2] ≈ res5[1]
    @test res0[4:6,1:2] ≈ res1[2]
    @test res0[4:6,1:2] ≈ res2[2]
    @test res0[4:6,1:2] ≈ res3[2]
    @test res0[4:6,1:2] ≈ res4[2]
    @test res0[4:6,1:2] ≈ res5[2]
end

@testset "mexp2" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        0.5 0.0;
        0.5 0.5;
        0.0 0.5
    ]
    ts = LinRange(0.0, 10.0, 10)
    res0 = [exp(Q*t) * x0 for t = ts]
    res1 = mexp(Q, x0, ts)
    res2 = mexp(SparseCSR(Q), x0, ts)
    res3 = mexp(SparseCSC(Q), x0, ts)
    res4 = mexp(SparseCOO(Q), x0, ts)
    res5 = mexp(SparseMatrixCSC(Q), x0, ts)
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
    @time res0 = [exp(Q'*t) * x0 for t = ts]
    @time res0 = [exp(Q'*t) * x0 for t = ts]
    @time res1 = mexp(Q, x0, ts, transpose=:T)
    @time res1 = mexp(Q, x0, ts, transpose=:T)
    @time res2 = mexp(SparseCSR(Q), x0, ts, transpose=:T)
    @time res2 = mexp(SparseCSR(Q), x0, ts, transpose=:T)
    @time res3 = mexp(SparseCSC(Q), x0, ts, transpose=:T)
    @time res3 = mexp(SparseCSC(Q), x0, ts, transpose=:T)
    @time res4 = mexp(SparseCOO(Q), x0, ts, transpose=:T)
    @time res4 = mexp(SparseCOO(Q), x0, ts, transpose=:T)
    @time res5 = mexp(SparseMatrixCSC(Q), x0, ts, transpose=:T)
    @time res5 = mexp(SparseMatrixCSC(Q), x0, ts, transpose=:T)
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
end

@testset "mexpc2" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    I = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    x0 = [
        0.5 0.0;
        0.5 0.5;
        0.0 0.5
    ]
    ts = LinRange(0.0, 10.0, 10)

    Qdash = [
        Q zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res1 = mexpc(Q, x0, ts)
    @time res1 = mexpc(Q, x0, ts)
    @time res2 = mexpc(SparseCSR(Q), x0, ts)
    @time res2 = mexpc(SparseCSR(Q), x0, ts)
    @time res3 = mexpc(SparseCSC(Q), x0, ts)
    @time res3 = mexpc(SparseCSC(Q), x0, ts)
    @time res4 = mexpc(SparseCOO(Q), x0, ts)
    @time res4 = mexpc(SparseCOO(Q), x0, ts)
    @time res5 = mexpc(SparseMatrixCSC(Q), x0, ts)
    @time res5 = mexpc(SparseMatrixCSC(Q), x0, ts)
    @test res0[1] ≈ res1[1]
    @test res0[1] ≈ res2[1]
    @test res0[1] ≈ res3[1]
    @test res0[1] ≈ res4[1]
    @test res0[1] ≈ res5[1]
    @test res0[2] ≈ res1[2]
    @test res0[2] ≈ res2[2]
    @test res0[2] ≈ res3[2]
    @test res0[2] ≈ res4[2]
    @test res0[2] ≈ res5[2]

    Qdash = [
        Q' zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res1 = mexpc(Q, x0, ts, transpose=:T)
    @time res1 = mexpc(Q, x0, ts, transpose=:T)
    @time res2 = mexpc(SparseCSR(Q), x0, ts, transpose=:T)
    @time res2 = mexpc(SparseCSR(Q), x0, ts, transpose=:T)
    @time res3 = mexpc(SparseCSC(Q), x0, ts, transpose=:T)
    @time res3 = mexpc(SparseCSC(Q), x0, ts, transpose=:T)
    @time res4 = mexpc(SparseCOO(Q), x0, ts, transpose=:T)
    @time res4 = mexpc(SparseCOO(Q), x0, ts, transpose=:T)
    @time res5 = mexpc(SparseMatrixCSC(Q), x0, ts, transpose=:T)
    @time res5 = mexpc(SparseMatrixCSC(Q), x0, ts, transpose=:T)
    @test res0[1] ≈ res1[1]
    @test res0[1] ≈ res2[1]
    @test res0[1] ≈ res3[1]
    @test res0[1] ≈ res4[1]
    @test res0[1] ≈ res5[1]
    @test res0[2] ≈ res1[2]
    @test res0[2] ≈ res2[2]
    @test res0[2] ≈ res3[2]
    @test res0[2] ≈ res4[2]
    @test res0[2] ≈ res5[2]
end

@testset "mixexp1" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    ans = [
        0.2510060736298584394000 0.1554210097708885740531 0.2424928563850387575052;
        0.1989569330779636657791 0.2168314920514892163439 0.4320884477342963880808;
        0.1989569330779636657791 0.1686453260298337586409 0.4802746137559519845617
    ]
    @time res1 = mexp(Q, x0, Weibull(2.0, 1.0))
    @time res1 = mexp(Q, x0, Weibull(2.0, 1.0))
    for i = eachindex(ans)
        @test abs(ans[i] - res1[i]) < 1.0e-8
    end
end

@testset "mixexpc1" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    ans_y = [
        0.3209844251001398029999 0.2157169520527697503809 0.4632986228470717948724;
        0.2263385249666138854252 0.2406732432647746833254 0.5329882317685930015472;
        0.2263385249666138854252 0.1924870824122626444819 0.5811743926211048183461
    ]
    ans_cy = [
        0.3913106241978086541344 0.1769700990060719680841 0.3179461986875124490659;
        0.1649720992311948519760 0.3017542963096322417016 0.4195005263505657833178;
        0.1649720992311948519760 0.1431186564517208936742 0.5781361662084771868564
    ]
    # println(Q)
    # println(x0)
    @time y, cy = mexpc(Q, x0, Weibull(2.0, 1.0))
    # println(Q)
    # println(x0)
    @time y, cy = mexpc(Q, x0, Weibull(2.0, 1.0))
    for i = eachindex(ans_y)
        @test abs(ans_y[i] - y[i]) < 1.0e-8
    end
    for i = eachindex(ans_cy)
        @test abs(ans_cy[i] - cy[i]) < 1.0e-8
    end
end

@testset "conv" begin
    @testset "conv notrans notrans" begin
        Q = [
            -3.0 2.0 0.0;
            1.0 -5.0 4.0;
            0.0 1.0 -2.0
        ]
        x = rand(3)
        y = rand(3)
        A = [
            Q reshape(x, (3,1)) * reshape(y, (1,3));
            zeros(3,3) Q
        ]

        tau = 0.5
        X = exp(A*tau)
        z0 = X[1:3,1:3] * x

        P, qv = unif(Q)
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(Q)
        convunifstep!(:N, :N,
            P, poi, (0, right), weight, qv * weight, copy(x), y, z, H)
        @test z ≈ z0
        @test X[1:3,4:6] ≈ H

        P, qv = unif(SparseCSR(Q))
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(P)
        convunifstep!(:N, :N,
            P, poi, (0, right), weight, qv * weight, copy(x), y, z, H)
        @test z ≈ z0
        for i = 1:3
            for z = H.rowptr[i]:H.rowptr[i+1]-1
                j = H.colind[z]
                @test X[i, j+3] ≈ H.val[z]
            end
        end

        P, qv = unif(SparseCSC(Q))
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(P)
        convunifstep!(:N, :N,
            P, poi, (0, right), weight, qv * weight, copy(x), y, z, H)
        @test z ≈ z0
        for j = 1:3
            for z = H.colptr[j]:H.colptr[j+1]-1
                i = H.rowind[z]
                @test X[i, j+3] ≈ H.val[z]
            end
        end

        P, qv = unif(SparseCOO(Q))
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(P)
        convunifstep!(:N, :N,
            P, poi, (0, right), weight, qv * weight, copy(x), y, z, H)
        @test z ≈ z0
        for z = 1:nnz(H)
            i = H.rowind[z]
            j = H.colind[z]
            @test X[i, j+3] ≈ H.val[z]
        end

        P, qv = unif(SparseMatrixCSC(Q))
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(P)
        convunifstep!(:N, :N,
            P, poi, (0, right), weight, qv * weight, copy(x), y, z, H)
        @test z ≈ z0
        for j = 1:3
            for z = H.colptr[j]:H.colptr[j+1]-1
                i = H.rowval[z]
                @test X[i, j+3] ≈ H.nzval[z]
            end
        end
    end

    @testset "conv notrans trans" begin
        Q = [
            -3.0 2.0 0.0;
            1.0 -5.0 4.0;
            0.0 1.0 -2.0
        ]
        x = rand(3)
        y = rand(3)
        A = [
            Q reshape(x, (3,1)) * reshape(y, (1,3));
            zeros(3,3) Q
        ]

        tau = 0.5
        X = exp(A*tau)
        z0 = X[1:3,1:3] * x

        P, qv = unif(Q)
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(Q)
        convunifstep!(:N, :T,
            P, poi, (0, right), weight, qv * weight, x, y, z, H)
        @test z ≈ z0
        @test X[1:3,4:6]' ≈ H
    end

    @testset "conv trans notrans" begin
        Q = [
            -3.0 2.0 0.0;
            1.0 -5.0 4.0;
            0.0 1.0 -2.0
        ]
        x = rand(3)
        y = rand(3)
        A = [
            Q' reshape(x, (3,1)) * reshape(y, (1,3));
            zeros(3,3) Q'
        ]

        tau = 0.5
        X = exp(A*tau)
        z0 = X[1:3,1:3] * x

        P, qv = unif(Q)
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(Q)
        convunifstep!(:T, :N,
            P, poi, (0, right), weight, qv * weight, x, y, z, H)
        @test z ≈ z0
        @test X[1:3,4:6] ≈ H
    end

    @testset "conv trans trans" begin
        Q = [
            -3.0 2.0 0.0;
            1.0 -5.0 4.0;
            0.0 1.0 -2.0
        ]
        x = rand(3)
        y = rand(3)
        A = [
            Q' reshape(x, (3,1)) * reshape(y, (1,3));
            zeros(3,3) Q'
        ]

        tau = 0.5
        X = exp(A*tau)
        z0 = X[1:3,1:3] * x

        P, qv = unif(Q)
        right = rightbound(qv*tau, 1.0e-8)
        weight, poi = poipmf(qv*tau, right, left = 0)
        z = zeros(3)
        H = zero(Q)
        convunifstep!(:T, :T,
            P, poi, (0, right), weight, qv * weight, x, y, z, H)
        @test z ≈ z0
        @test X[1:3,4:6]' ≈ H
    end
end

@testset "tran 1" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x = Float64[1, 0, 0]
    r = Float64[1, 1, 0]

    ts = LinRange(0.0, 10.0, 10)
    @time irwd, crwd, y, cy = tran(Q, x, r, ts, forward=:T)

    res,cres = mexpc(Q, x, ts, transpose=:T)
    @test isapprox(irwd, [sum(x .* r) for x = res])
    @test isapprox(crwd, [sum(x .* r) for x = cres])
end

@testset "tran 2" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x = Float64[1, 0, 0]
    r = Float64[1, 1, 0]

    ts = LinRange(0.0, 10.0, 10)
    @time irwd, crwd, y, cy = tran(Q, x, r, ts, forward=:N)

    res,cres = mexpc(Q, x, ts, transpose=:T)
    @test isapprox(irwd, [sum(x .* r) for x = res])
    @test isapprox(crwd, [sum(x .* r) for x = cres])
end

@testset "tran 3" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x = Float64[1 0 0; 0 1 0]
    r = Float64[1, 1, 0]

    ts = LinRange(0.0, 10.0, 10)
    @time irwd, crwd, y, cy = tran(Q, x, r, ts, forward=:T)

    res,cres = mexpc(Q, r, ts, transpose=:N)
    @test isapprox(irwd, [x* r for r = res])
    @test isapprox(crwd, [x* r for r = cres])
end

@testset "tran 4" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x = Float64[1 0 0; 0 1 0]
    r = Float64[1, 1, 0]

    ts = LinRange(0.0, 10.0, 10)
    @time irwd, crwd, y, cy = tran(Q, x, r, ts, forward=:N)

    res,cres = mexpc(Q, r, ts, transpose=:N)
    @test isapprox(irwd, [x* r for r = res])
    @test isapprox(crwd, [x* r for r = cres])
end

@testset "tran 5" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x = Float64[1 0 0; 0 1 0]
    r = Float64[1 1 0; 0 1 0]

    ts = LinRange(0.0, 10.0, 10)
    @time irwd, crwd, y, cy = tran(Q, x, r', ts, forward=:N)

    res,cres = mexpc(Q, r', ts, transpose=:N)
    @test isapprox(irwd, [x * r for r = res])
    @test isapprox(crwd, [x * r for r = cres])
end

@testset "tran 6" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x = Float64[1 0 0; 0 1 0]
    r = Float64[1 1 0; 0 1 0]

    ts = LinRange(0.0, 10.0, 10)
    @time irwd, crwd, y, cy = tran(Q, x, r', ts, forward=:T)

    res,cres = mexpc(Q, r', ts, transpose=:N)
    @test isapprox(irwd, [x * r for r = res])
    @test isapprox(crwd, [x * r for r = cres])
end

@testset "mexp3" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        0.5 0.5 0.0;
        0.5 0.5 0.5
    ]
    t = 10.0
    res0 = exp(Q*t) * x0'
    res1 = mexp(Q, x0', t)
    res2 = mexp(SparseCSR(Q), x0', t)
    res3 = mexp(SparseCSC(Q), x0', t)
    res4 = mexp(SparseCOO(Q), x0', t)
    res5 = mexp(SparseMatrixCSC(Q), x0', t)
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
    res0 = exp(Q'*t) * x0'
    res1 = mexp(Q, x0', t, transpose=:T)
    res2 = mexp(SparseCSR(Q), x0', t, transpose=:T)
    res3 = mexp(SparseCSC(Q), x0', t, transpose=:T)
    res4 = mexp(SparseCOO(Q), x0', t, transpose=:T)
    res5 = mexp(SparseMatrixCSC(Q), x0', t, transpose=:T)
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
end

@testset "mexp4" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        0.5 0.5 0.0;
        0.5 0.5 0.5
    ]
    ts = LinRange(0.0, 10.0, 10)
    res0 = [exp(Q*t) * x0' for t = ts]
    res1 = mexp(Q, x0', ts)
    res2 = mexp(SparseCSR(Q), x0', ts)
    res3 = mexp(SparseCSC(Q), x0', ts)
    res4 = mexp(SparseCOO(Q), x0', ts)
    res5 = mexp(SparseMatrixCSC(Q), x0', ts)
    # println(typeof(res1))
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
    @time res0 = [exp(Q'*t) * x0' for t = ts]
    @time res0 = [exp(Q'*t) * x0' for t = ts]
    @time res1 = mexp(Q, x0', ts, transpose=:T)
    @time res1 = mexp(Q, x0', ts, transpose=:T)
    @time res2 = mexp(SparseCSR(Q), x0', ts, transpose=:T)
    @time res2 = mexp(SparseCSR(Q), x0', ts, transpose=:T)
    @time res3 = mexp(SparseCSC(Q), x0', ts, transpose=:T)
    @time res3 = mexp(SparseCSC(Q), x0', ts, transpose=:T)
    @time res4 = mexp(SparseCOO(Q), x0', ts, transpose=:T)
    @time res4 = mexp(SparseCOO(Q), x0', ts, transpose=:T)
    @time res5 = mexp(SparseMatrixCSC(Q), x0', ts, transpose=:T)
    @time res5 = mexp(SparseMatrixCSC(Q), x0', ts, transpose=:T)
    @test res0 ≈ res1
    @test res0 ≈ res2
    @test res0 ≈ res3
    @test res0 ≈ res4
    @test res0 ≈ res5
end

@testset "mexpc3" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    I = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    x0 = [
        0.5 0.0;
        0.5 0.5;
        0.0 0.5
    ]
    x1 = copy(x0')
    t = 10.0

    Qdash = [
        Q zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    res0 = exp(Qdash*t) * x0dash
    res1 = mexpc(Q, x1', t)
    res2 = mexpc(SparseCSR(Q), x1', t)
    res3 = mexpc(SparseCSC(Q), x1', t)
    res4 = mexpc(SparseCOO(Q), x1', t)
    res5 = mexpc(SparseMatrixCSC(Q), x1', t)
    @test res0[1:3,1:2] ≈ res1[1]
    @test res0[1:3,1:2] ≈ res2[1]
    @test res0[1:3,1:2] ≈ res3[1]
    @test res0[1:3,1:2] ≈ res4[1]
    @test res0[1:3,1:2] ≈ res5[1]
    @test res0[4:6,1:2] ≈ res1[2]
    @test res0[4:6,1:2] ≈ res2[2]
    @test res0[4:6,1:2] ≈ res3[2]
    @test res0[4:6,1:2] ≈ res4[2]
    @test res0[4:6,1:2] ≈ res5[2]

    Qdash = [
        Q' zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    res0 = exp(Qdash*t) * x0dash
    res1 = mexpc(Q, x1', t, transpose=:T)
    res2 = mexpc(SparseCSR(Q), x1', t, transpose=:T)
    res3 = mexpc(SparseCSC(Q), x1', t, transpose=:T)
    res4 = mexpc(SparseCOO(Q), x1', t, transpose=:T)
    res5 = mexpc(SparseMatrixCSC(Q), x1', t, transpose=:T)
    @test res0[1:3,1:2] ≈ res1[1]
    @test res0[1:3,1:2] ≈ res2[1]
    @test res0[1:3,1:2] ≈ res3[1]
    @test res0[1:3,1:2] ≈ res4[1]
    @test res0[1:3,1:2] ≈ res5[1]
    @test res0[4:6,1:2] ≈ res1[2]
    @test res0[4:6,1:2] ≈ res2[2]
    @test res0[4:6,1:2] ≈ res3[2]
    @test res0[4:6,1:2] ≈ res4[2]
    @test res0[4:6,1:2] ≈ res5[2]
end

@testset "mexpc4" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    I = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    x0 = [
        0.5 0.0;
        0.5 0.5;
        0.0 0.5
    ]
    x1 = copy(x0')
    ts = LinRange(0.0, 10.0, 10)

    Qdash = [
        Q zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res1 = mexpc(Q, x1', ts)
    @time res1 = mexpc(Q, x1', ts)
    @time res2 = mexpc(SparseCSR(Q), x1', ts)
    @time res2 = mexpc(SparseCSR(Q), x1', ts)
    @time res3 = mexpc(SparseCSC(Q), x1', ts)
    @time res3 = mexpc(SparseCSC(Q), x1', ts)
    @time res4 = mexpc(SparseCOO(Q), x1', ts)
    @time res4 = mexpc(SparseCOO(Q), x1', ts)
    @time res5 = mexpc(SparseMatrixCSC(Q), x1', ts)
    @time res5 = mexpc(SparseMatrixCSC(Q), x1', ts)
    @test res0[1] ≈ res1[1]
    @test res0[1] ≈ res2[1]
    @test res0[1] ≈ res3[1]
    @test res0[1] ≈ res4[1]
    @test res0[1] ≈ res5[1]
    @test res0[2] ≈ res1[2]
    @test res0[2] ≈ res2[2]
    @test res0[2] ≈ res3[2]
    @test res0[2] ≈ res4[2]
    @test res0[2] ≈ res5[2]

    Qdash = [
        Q' zeros(3,3);
        I zeros(3,3)
    ]
    x0dash = [
        x0;
        zeros(3,2)
    ]
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res0 = ([(exp(Qdash*t) * x0dash)[1:3,1:2] for t = ts], [(exp(Qdash*t) * x0dash)[4:6,1:2] for t = ts])
    @time res1 = mexpc(Q, x1', ts, transpose=:T)
    @time res1 = mexpc(Q, x1', ts, transpose=:T)
    @time res2 = mexpc(SparseCSR(Q), x1', ts, transpose=:T)
    @time res2 = mexpc(SparseCSR(Q), x1', ts, transpose=:T)
    @time res3 = mexpc(SparseCSC(Q), x1', ts, transpose=:T)
    @time res3 = mexpc(SparseCSC(Q), x1', ts, transpose=:T)
    @time res4 = mexpc(SparseCOO(Q), x1', ts, transpose=:T)
    @time res4 = mexpc(SparseCOO(Q), x1', ts, transpose=:T)
    @time res5 = mexpc(SparseMatrixCSC(Q), x1', ts, transpose=:T)
    @time res5 = mexpc(SparseMatrixCSC(Q), x1', ts, transpose=:T)
    @test res0[1] ≈ res1[1]
    @test res0[1] ≈ res2[1]
    @test res0[1] ≈ res3[1]
    @test res0[1] ≈ res4[1]
    @test res0[1] ≈ res5[1]
    @test res0[2] ≈ res1[2]
    @test res0[2] ≈ res2[2]
    @test res0[2] ≈ res3[2]
    @test res0[2] ≈ res4[2]
    @test res0[2] ≈ res5[2]
end

@testset "mixexp2" begin
    Q = [
        -3.0 2.0 0.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    x1 = copy(x0')
    ans = [
        0.2510060736298584394000 0.1554210097708885740531 0.2424928563850387575052;
        0.1989569330779636657791 0.2168314920514892163439 0.4320884477342963880808;
        0.1989569330779636657791 0.1686453260298337586409 0.4802746137559519845617
    ]
    @time res1 = mexp(Q, x1', Weibull(2.0, 1.0))
    @time res1 = mexp(Q, x1', Weibull(2.0, 1.0))
    for i = eachindex(ans)
        @test abs(ans[i] - res1[i]) < 1.0e-8
    end
end

@testset "mixexpc2" begin
    Q = [
        -3.0 2.0 1.0;
        1.0 -5.0 4.0;
        1.0 1.0 -2.0
    ]
    x0 = [
        1.0 0.0 0.0;
        0.0 1.0 0.0;
        0.0 0.0 1.0
    ]
    x1 = copy(x0')
    ans_y = [
        0.3209844251001398029999 0.2157169520527697503809 0.4632986228470717948724;
        0.2263385249666138854252 0.2406732432647746833254 0.5329882317685930015472;
        0.2263385249666138854252 0.1924870824122626444819 0.5811743926211048183461
    ]
    ans_cy = [
        0.3913106241978086541344 0.1769700990060719680841 0.3179461986875124490659;
        0.1649720992311948519760 0.3017542963096322417016 0.4195005263505657833178;
        0.1649720992311948519760 0.1431186564517208936742 0.5781361662084771868564
    ]
    @time y, cy = mexpc(Q, x1', Weibull(2.0, 1.0))
    @time y, cy = mexpc(Q, x1', Weibull(2.0, 1.0))
    for i = eachindex(ans_y)
        @test abs(ans_y[i] - y[i]) < 1.0e-8
    end
    for i = eachindex(ans_cy)
        @test abs(ans_cy[i] - cy[i]) < 1.0e-8
    end
end

@testset "zeromat" begin
    Q = [0.0][:,:]
    x0 = eye(1)
    y = mexp(Q, x0, 1.0)
    @test y == [1.0][:,:]
    y, cy = mexpc(Q, x0, 10.0)
    @test y == [1.0][:,:]
    @test isapprox(cy, [10.0][:,:])
end

@testset "zeromat_mix" begin
    Q = [0.0][:,:]
    x0 = eye(1)
    y = mexpmix(Q, x0) do t
        exp(-t)
    end
    @test isapprox(y, [1.0][:,:])
    y, cy = mexpcmix(Q, x0, bounds=(0.0, 10.0)) do t
        exp(-t)
    end
    println(cy)
    @test isapprox(y[1,1], 1-exp(-10.0))
    @test isapprox(cy[1,1], (1-(1+10.0*1.0)*exp(-10.0*1.0))/1.0)
end