text:
"""
    run_omniscape("path/to/config.ini")

Run the Omniscape algorithm using the files and options specified in a .ini file

Visit https://circuitscape.github.io/Omniscape.jl/stable for detailed
documentation
"""
function run_omniscape(path::String)
    start_time = time()

    cfg_user = parse_cfg(path)
    cfg = init_cfg()
    update_cfg!(cfg, cfg_user)

    ## Parse commonly called integer arguments
    int_arguments = Dict{String, Int64}()

    int_arguments["block_size"] = Int64(round(parse(Float64,
                                                    cfg["block_size"])))

    if iseven(int_arguments["block_size"])
        @warn "block_size is even, but must be odd. Using block_size + 1."
        int_arguments["block_size"] = int_arguments["block_size"] + 1
    end

    int_arguments["block_radius"] = Int64((int_arguments["block_size"] - 1) / 2)
    int_arguments["radius"] = Int64(round(parse(Float64, cfg["radius"])))
    int_arguments["buffer"] = Int64(round(parse(Float64, cfg["buffer"])))

    ## Parse other arguments
    # flags
    calc_flow_potential = lowercase(cfg["calc_flow_potential"]) == "true"
    write_flow_potential = lowercase(cfg["write_flow_potential"]) == "true"
    write_normalized_currmap = lowercase(cfg["write_normalized_currmap"]) == "true"
    write_raw_currmap = lowercase(cfg["write_raw_currmap"]) == "true"
    parallelize = lowercase(cfg["parallelize"]) == "true"
    correct_artifacts = lowercase(cfg["correct_artifacts"]) == "true"
    source_from_resistance = lowercase(cfg["source_from_resistance"]) == "true"
    if int_arguments["block_size"] == 1
        correct_artifacts = false
    end

    # other
    source_threshold = parse(Float64, cfg["source_threshold"])
    project_name = cfg["project_name"]
    n_workers = parse(Int64, cfg["max_parallel"])
    r_cutoff = parse(Float64, cfg["r_cutoff"])

    ## Store ascii header
    final_header = parse_ascii_header("$(cfg["resistance_file"])")

    ## Import sources and resistances
    resistance_raw = float(read_ascii("$(cfg["resistance_file"])"))

    # Adjust nodata value if not already -9999
    if parse(Float64, final_header["nodata_value"]) != -9999
        no_data_val = parse(Float64, final_header["nodata_value"])
        resistance_raw[resistance_raw .== no_data_val] .= -9999
        final_header["nodata_value"] = "-9999"
    end

    if source_from_resistance
        sources_raw = deepcopy(resistance_raw)
        sources_raw = 1.0 ./ sources_raw
        sources_raw[resistance_raw .> r_cutoff] .= 0.0
        sources_raw[resistance_raw .== -9999] .= 0.0
    else
        sources_raw = float(read_ascii("$(cfg["source_file"])"))
    end

    int_arguments["nrows"] = size(sources_raw, 1)
    int_arguments["ncols"] = size(sources_raw, 2)

    ## Setup Circuitscape configuration
    cs_cfg_dict = init_csdict(cfg)
    cs_cfg = Circuitscape.init_config()
    Circuitscape.update!(cs_cfg, cs_cfg_dict)

    ## Calculate targets
    targets = get_targets(sources_raw,
                          int_arguments,
                          threshold = source_threshold)

    ## Initialize temporary ascii header for CS advanced mode
    temp_header = init_ascii_header()

    ## Circuitscape calls in loop over targets
    n_targets = size(targets, 1)

    ## Define parameters for cs
    # Get flags
    o = Circuitscape.OutputFlags(false, false,
                                 false, false,
                                 false, false,
                                 false, false)

    ## Add parallel workers
    if parallelize
        println("Starting up Omniscape to use $(n_workers) processes in parallel")
        myaddprocs(n_workers)

        for i in workers()
            remotecall(copyvars, i, sources_raw)
        end

        for i in workers()
            @spawnat i eval(:(cum_currmap = fill(0.,
                                                 nrows,
                                                 ncols)))
        end

        if calc_flow_potential
            for i in workers()
                @spawnat i eval(:(fp_cum_currmap = fill(0.,
                                                        nrows,
                                                        ncols)))
            end
        end
    else
        cum_currmap = fill(0., int_arguments["nrows"], int_arguments["ncols"])
        if calc_flow_potential
            fp_cum_currmap = fill(0.,
                                  int_arguments["nrows"],
                                  int_arguments["ncols"])
        else
            fp_cum_currmap = Array{Float64, 2}(undef, 1, 1)
        end
    end

    if correct_artifacts
        art_start = time()
        correction_array = calc_correction(int_arguments,
                                           cs_cfg,
                                           o)
        println("time taken to calculate artifact correction array: $(time()- art_start) seconds")
    else
        correction_array = Array{Float64, 2}(undef, 1, 1)
    end

    ## Calculate and accumulate currents on each worker
    println("Solving targets")
    if parallelize
    pmap(x -> solve_target!(x, n_targets, int_arguments, targets,
                            sources_raw, resistance_raw, cs_cfg, o,
                            calc_flow_potential, correct_artifacts,
                            correction_array),
         1:n_targets)
    else
        for i in 1:n_targets
            solve_target!(i,
                          n_targets,
                          int_arguments,
                          targets,
                          sources_raw,
                          resistance_raw,
                          cs_cfg,
                          o,
                          calc_flow_potential,
                          correct_artifacts,
                          correction_array,
                          cum_currmap,
                          fp_cum_currmap)
        end
    end

    ## Add together remote cumulative maps
    if parallelize
        println("Combining maps across workers")

        cum_currmap = sum_currmaps(int_arguments)

        if calc_flow_potential
            fp_cum_currmap = sum_fpmaps(int_arguments)
        end
    end

    if calc_flow_potential == true
        normalized_cum_currmap = cum_currmap ./ fp_cum_currmap
    end

    ## Make output directory
    mkdir("$(project_name)_output")

    ## Write outputs
    if write_raw_currmap == true
        write_ascii(cum_currmap,
                    "$(project_name)_output/cum_currmap.asc",
                    final_header)
    end

    if calc_flow_potential == true
        if write_flow_potential == true
            write_ascii(fp_cum_currmap,
                        "$(project_name)_output/flow_potential.asc",
                        final_header)
        end

        if write_normalized_currmap == true
            write_ascii(normalized_cum_currmap,
                        "$(project_name)_output/normalized_cum_currmap.asc",
                        final_header)
        end
    end

    if parallelize
        rmprocs(workers())
    end

    println("Done")
    println("Time taken to complete job: $(round(time() - start_time; digits = 1)) seconds")

    ## Return outputs
    if calc_flow_potential == true
        return normalized_cum_currmap, cum_currmap, fp_cum_currmap
    else
        return cum_currmap
    end
end


text:
using Test
using Printf
using Random
using Statistics
using LinearAlgebra
using Logging

using CUDA
using MPI
using JLD2
using FFTW
using OffsetArrays
using SeawaterPolynomials

using Oceananigans
using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.Advection
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.AbstractOperations
using Oceananigans.Coriolis
using Oceananigans.BuoyancyModels
using Oceananigans.Forcings
using Oceananigans.Solvers
using Oceananigans.Models
using Oceananigans.Simulations
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.TurbulenceClosures
using Oceananigans.Distributed
using Oceananigans.Logger
using Oceananigans.Units
using Oceananigans.Utils
using Oceananigans.MultiRegion
using Oceananigans.Architectures: device, array_type # to resolve conflict with CUDA.device

using Oceananigans: Clock
using Dates: DateTime, Nanosecond
using TimesDates: TimeDate
using Statistics: mean
using LinearAlgebra: norm
using NCDatasets: Dataset
using KernelAbstractions: @kernel, @index, Event

import Oceananigans.Fields: interior
import Oceananigans.Utils: launch!, datatuple

Logging.global_logger(OceananigansLogger())

#####
##### Testing parameters
#####

float_types = (Float32, Float64)

closures = (
    :ScalarDiffusivity,
    :ScalarBiharmonicDiffusivity,
    :TwoDimensionalLeith,
    :SmagorinskyLilly,
    :AnisotropicMinimumDissipation,
    :ConvectiveAdjustmentVerticalDiffusivity
)

#####
##### Run tests!
#####

CUDA.allowscalar(true)

float_types = (Float32, Float64)

include("data_dependencies.jl")
include("utils_for_runtests.jl")

archs = test_architectures()

group     = get(ENV, "TEST_GROUP", :all) |> Symbol
test_file = get(ENV, "TEST_FILE", :none) |> Symbol

ext: jl


text:
#Running this code will take a while, it is included here for illustration purposed. The resulting circuit library CSV file (Circuitlibrary.csv) can already conveniently
#be found in the root folder.

initialize_circuitlibrary("R1-[C2,R3-[C4,R5]]",[20,4e-9,3400,4e-6,2500],"Nonbiological","https://www.gamry.com/application-notes/EIS/basics-of-electrochemical-impedance-spectroscopy/")
add_to_circuitlibrary("[R1,[C2,R3]-R4]",[4310,3.5e-9,2.9e6,220],"Biological","Electrical Impedance Studies on Potato and Alfalfa Tissue, Hayden et al")
add_to_circuitlibrary("[R1,C2-[R3,C4-R5]]",[7355,5.11e-9,3274,5e-9,346],"Biological","Electrical Impedance Analysis in Plant Tissues: A Double Shell Model, Zhang et al")
add_to_circuitlibrary("R1-[C2,R3]-[C4,R5]","Biological","Bioimpedance analysis of vascular tissue and fluid flow in human and plant body: A review, Prasad et al (plant cross-section)")
add_to_circuitlibrary("[R1-C2,R3]","Biological","Bioimpedance analysis of vascular tissue and fluid flow in human and plant body: A review, Prasad et al (body segment tissue model)")
add_to_circuitlibrary("R1-L2-[C3,R4]","Biological","Bioimpedance analysis of vascular tissue and fluid flow in human and plant body: A review, Prasad et al (arterial system RLCR model)")
add_to_circuitlibrary("P1-R2-[C3,R4]","Biological","Electrochemical impedance spectroscopy as a versatile tool for the characterization of neural tissue: A mini review, Krukiewicz et al")
add_to_circuitlibrary("[R1,C2-R3-[C4,R5]]",[2.6e3,1e-6,470,10,10e3],"Biological","Impedance spectroscopy on living tissue for determination of the state of organs, Gersing et al (porcine liver)")
add_to_circuitlibrary("[R1-[P2,R3],C4]","Biological","A suggested circuit for bacterial biofilms (Erkuden)")
add_to_circuitlibrary("R1-[R2,C3]",[20,250,40e-6],"Nonbiological","Randles circuit, https://www.gamry.com/application-notes/EIS/basics-of-electrochemical-impedance-spectroscopy/")
add_to_circuitlibrary("R1-[C2,R3-P4]",[20,40e-6,250,[150,0.5]],"Nonbiological"," Mixed Kinetic and Charge-Transfer Control Randles cell, https://www.gamry.com/application-notes/EIS/basics-of-electrochemical-impedance-spectroscopy/")
add_to_circuitlibrary("R1-C2-[R3,P4]","Biological"," Rapid detection of bacterial proliferation in food samples,Puttaswamy et al")
add_to_circuitlibrary("R1-[P2,R3]",[234,[2e8,0.56],15789],"Biological"," Electrical impedance spectroscopic study of mandarin orange during ripening, Chowdhury et al")
add_encoding_to_circuitlibrary("+++--+PPRPRLR","Nonbiological","Characterization of high-power lithium-ion batteries by electrochemical impedance spectroscopy, Andre et al")
add_to_circuitlibrary("[C1,R2-C3]",[2.1e-12,1.6e3,2.2e-12],"Nonbiological","Impedance spectroscopy study of hardened Portland cement paste, cabeza et al")


text:
using Iris
using Test

@testset "Linear Eigenvalue Problem" begin
    @testset "Inference" begin
        ω = 10
        sim = Iris.Examples.slab_pml(ω=ω)
        lep = @test_nowarn HelmholtzLEP(sim)
        cf = @test_nowarn HelmholtzCF(sim)
        nep = @test_nowarn HelmholtzNEP(sim)
        f(x) = getproperty(x,:helmholtz)
        @test_nowarn @inferred f(lep)
        @test_nowarn @inferred f(cf)
        @test_nowarn @inferred f(nep)
        f(x) = getproperty(x,:saturated)
        @test_nowarn @inferred f(lep)
        @test_nowarn @inferred f(cf)
        @test_nowarn @inferred f(nep)
        f(x) = getproperty(x,:αεpFχ)
        @test_nowarn @inferred f(lep)
        f(x) = getproperty(x,:F)
        @test_nowarn @inferred f(cf)
        f(x) = getproperty(x,:Fs)
        @test_nowarn @inferred f(nep)
        f(x) = getproperty(x,:nep)
        @test_nowarn @inferred f(nep)
        f(x) = getproperty(x,:simulation)
        @test_nowarn @inferred f(lep)
        @test_nowarn @inferred f(cf)
        @test_nowarn @inferred f(nep)
        @test_nowarn @inferred ScalarField(lep)
        @test_nowarn @inferred ScalarField(cf)
        @test_nowarn @inferred ScalarField(nep)
        @test_nowarn @inferred VectorField(lep)
        @test_nowarn @inferred VectorField(cf)
        @test_nowarn @inferred VectorField(nep)
        @test_nowarn @inferred ScalarField(lep,5)
        @test_nowarn @inferred ScalarField(cf,5)
        @test_nowarn @inferred ScalarField(nep,5)
        @test_nowarn @inferred VectorField(5,lep)
        @test_nowarn @inferred VectorField(5,cf)
        @test_nowarn @inferred VectorField(5,nep)
        @test_nowarn @inferred ScalarField(lep,rand(length(sim.x),12))
        @test_nowarn @inferred ScalarField(cf,rand(length(sim.x),2))
        @test_nowarn @inferred ScalarField(nep,rand(length(sim.x),3))
        @test_nowarn @inferred VectorField(lep,rand(length(sim.x),12))
        @test_nowarn @inferred VectorField(cf,rand(length(sim.x),2))
        @test_nowarn @inferred VectorField(nep,rand(length(sim.x),3))
        @test_nowarn @inferred ScalarField(lep,rand(length(sim.x),12))
        @test_nowarn @inferred ScalarField(rand(length(sim.x),2), cf)
        @test_nowarn @inferred ScalarField(rand(length(sim.x),3), nep)
        @test_nowarn @inferred VectorField(rand(length(sim.x),12), lep)
        @test_nowarn @inferred VectorField(rand(length(sim.x),2), cf)
        @test_nowarn @inferred VectorField(rand(length(sim.x),3), nep)
        @test_nowarn @inferred lep(ω)
        @test_nowarn @inferred lep(ω,[ω],ScalarField(lep))
        @test_nowarn @inferred cf(ω)
        @test_nowarn @inferred cf(ω,[ω],ScalarField(cf))
        @test_nowarn @inferred nep(ω)
        @test_nowarn @inferred nep(ω,[ω],ScalarField(nep))
    end
    @testset "Helmholtz 1D" begin
        @testset "PML" begin
            ω0 = 10
            sim = Iris.Examples.slab_pml(ω=ω0)
            lep = @test_nowarn HelmholtzLEP(sim)
            ω, ψ = @test_nowarn helmholtzeigen(lep,ω0;nev=4)
            @test_nowarn @inferred helmholtzeigen(lep,ω0;nev=4)
            sort!(ω; by=real)
            @test isapprox(ω[2] - ω[1], π/3; atol=2e-2)
            @test isapprox(ω[3] - ω[2], π/3; atol=2e-2)
            @test isapprox(ω[4] - ω[3], π/3; atol=2e-2)
            @test isapprox(imag(ω[1]), -.23; atol=2e-2)
            @test isapprox(imag(ω[2]), -.23; atol=2e-2)
            @test isapprox(imag(ω[3]), -.23; atol=2e-2)
            @test isapprox(imag(ω[4]), -.23; atol=2e-2)
        end
        @testset "cPML" begin
            ω0 = 10
            sim = Iris.Examples.slab_cpml(ω=ω0)
            lep = @test_nowarn HelmholtzLEP(sim)
            ω, ψ = @test_nowarn helmholtzeigen(lep,ω0;nev=4)
            @test_nowarn @inferred helmholtzeigen(lep,ω0;nev=4)
            sort!(ω; by=real)
            @test isapprox(ω[2] - ω[1], π/3; atol=2e-2)
            @test isapprox(ω[3] - ω[2], π/3; atol=2e-2)
            @test isapprox(ω[4] - ω[3], π/3; atol=2e-2)
            @test isapprox(imag(ω[1]), +.23; atol=2e-2)
            @test isapprox(imag(ω[2]), +.23; atol=2e-2)
            @test isapprox(imag(ω[3]), +.23; atol=2e-2)
            @test isapprox(imag(ω[4]), +.23; atol=2e-2)
        end
    end
    @testset "Helmholtz 2D" begin
        @testset "Dirichlet" begin
            bnd = Boundary(Square(π,0,0),DirichletBC)
            lat = Cartesian(.004,.004)
            latdom = LatticeDomain(bnd,lat)
            sim = Simulation(1,latdom)
            lep = @test_nowarn HelmholtzLEP(sim)
            ω, ψ = @test_nowarn helmholtzeigen(lep,1;nev=10)
            @test_nowarn @inferred helmholtzeigen(lep,1;nev=1)
            sort!(ω; by=real)
            @test isapprox(ω[1]^2, 1^2 + 1^2; atol = 3e-2)
            @test isapprox(ω[2]^2, 2^2 + 1^2; atol = 3e-2)
            @test isapprox(ω[3]^2, 1^2 + 2^2; atol = 3e-2)
            @test isapprox(ω[4]^2, 2^2 + 2^2; atol = 3e-2)
            @test isapprox(ω[5]^2, 1^2 + 3^2; atol = 3e-2)
            @test isapprox(ω[6]^2, 3^2 + 1^2; atol = 3e-2)
            @test isapprox(ω[7]^2, 2^2 + 3^2; atol = 3e-2)
            @test isapprox(ω[8]^2, 3^2 + 2^2; atol = 3e-2)
            @test isapprox(ω[9]^2, 1^2 + 4^2; atol = 3e-2)
            @test isapprox(ω[10]^2, 4^2 + 1^2; atol = 3e-2)
        end
    end
end
@testset "CF Eigenvalue Problem" begin
    @testset "Helmholtz 1D" begin
        @testset "PML" begin
            ω = 10.47
            sim = Iris.Examples.slab_pml(ω=ω)
            cf = @test_nowarn HelmholtzCF(sim)
            η, u = @test_nowarn helmholtzeigen(cf,ω; nev=4)
            @test_nowarn @inferred helmholtzeigen(cf,ω;nev=2)
            @test all(imag(η).<0)
            @test isapprox(η[1], 0-.4im; atol = 3e-2)
            @test isapprox(imag(η[2]), -.4; atol = 3e-2)
        end
        @testset "cPML" begin
            cavity = Interval(-.5,.5)
            pump = Interval(-.5,0)
            bnd = Boundary(Interval(-1,1),cPML)
            dom = NondispersiveDomain(cavity,3)
            pump = DispersiveDomain(Interval(-.8,.8),TwoLevelSystem())
            lat = Cartesian(.001)
            latdom = LatticeDomain(bnd,lat)
            sim = Simulation(10,dom,latdom,pump)
            cf = @test_nowarn HelmholtzCF(sim)
            η , u = @test_nowarn helmholtzeigen(cf,10;nev=4)
            @test_nowarn @inferred helmholtzeigen(cf,10;nev=2)
        end
    end
end
@testset "Nonlinear Eigenvalue Problem" begin
    @testset "Helmholtz 1D" begin
        @testset "PML" begin
            cavity = Interval(-.5,.5)
            pump = Interval(-.5,0)
            bnd = Boundary(Interval(-1,1),PML)
            dom = NondispersiveDomain(cavity,3)
            fdom = DispersiveDomain(pump,TwoLevelSystem(10,.1,2))
            lat = Cartesian(.001)
            latdom = LatticeDomain(bnd,lat)
            sim = Simulation(10,dom,latdom,fdom)
            nep = @test_nowarn HelmholtzNEP(sim)
            @test_nowarn nep(10)
            @test_nowarn @inferred nep(10)
            ω, ψ = @test_nowarn helmholtzeigen(nep,10;neigs=4)
            @test_nowarn @inferred helmholtzeigen(nep,10;neigs=2)
        end
        @testset "cPML" begin
            cavity = Interval(-.5,.5)
            pump = Interval(-.5,0)
            bnd = Boundary(Interval(-1,1),cPML)
            dom = NondispersiveDomain(cavity,3)
            fdom = DispersiveDomain(pump,TwoLevelSystem(10,.1,2))
            lat = Cartesian(.001)
            latdom = LatticeDomain(bnd,lat)
            sim = Simulation(10,dom,latdom,fdom)
            nep = @test_nowarn HelmholtzNEP(sim)
            @test_nowarn nep(10)
            @test_nowarn @inferred nep(10)
            ω, ψ = @test_nowarn helmholtzeigen(nep,10;neigs=4)
            @test_nowarn @inferred helmholtzeigen(nep,10;neigs=2)
        end
        @testset "MathcedBC sans dispersion" begin
            cavity = Interval(-.5,.5)
            pump = Interval(-.5,0)
            bnd = Boundary(Interval(-1,1),MatchedBC{1}(out=1),MatchedBC{2}(out=1))
            dom = NondispersiveDomain(cavity,3)
            lat = Cartesian(.001)
            latdom = LatticeDomain(bnd,lat)
            sim = Simulation(10,dom,latdom)
            nep = @test_nowarn HelmholtzNEP(sim)
            @test_nowarn @inferred HelmholtzNEP(sim)
            @test_nowarn nep(10)
            @test_nowarn @inferred nep(10)
            ω, ψ = @test_nowarn helmholtzeigen(nep,10;neigs=4,maxit=50)
            @test_nowarn @inferred helmholtzeigen(nep,10;neigs=1,maxit=50)
            sort!(ω; by=real)
            @test isapprox(ω[2] - ω[1], π/3; atol=1e-3)
            @test isapprox(ω[3] - ω[2], π/3; atol=1e-3)
            @test isapprox(ω[4] - ω[3], π/3; atol=1e-3)
        end
        @testset "MathcedBC with complex dispersion" begin
            cavity = Interval(-.5,.5)
            bnd = Boundary(Interval(-1,1),MatchedBC{1}(out=1),MatchedBC{2}(out=1))
            dom = NondispersiveDomain(cavity,3)
            fdom1 = DispersiveDomain(Interval(-.5,-.1),TwoLevelSystem(42,-.1,1))
            fdom2 = DispersiveDomain(Interval(.1,.5),TwoLevelSystem(45,.3,2))
            lat = Cartesian(.001)
            latdom = LatticeDomain(bnd,lat)
            sim = Simulation(40,dom,latdom,fdom1,fdom2)
            nep = @test_nowarn HelmholtzNEP(sim)
            @test_nowarn @inferred HelmholtzNEP(sim)
            @test_nowarn nep(40)
            @test_nowarn @inferred nep(40)
            ω, ψ = @test_nowarn helmholtzeigen(nep,43;neigs=4,maxit=50)
            @test_nowarn @inferred helmholtzeigen(nep,43;neigs=1,maxit=50)
        end
    end
end


text:
"""
    fit_single_peak_histogram(h::StatsBase.Histogram, fitrange::Union{Tuple{T,T}, Missing} = missing; fit_function = :Gauss, weights = missing )::FitFunction{T, 1, NP} where {T <: AbstractFloat, NP}

This method is meant for quick and stable fitting. The estimation of initial parameters is done automatically, this allows and easy fit with just one line of code. The user is expected to pass a histogram containing only one peak or specify corresponding fitrange.
The method performs a least square fit with the model specified under  `fit_function = ...`. Currently supported are :Gauss (default), :Gauss_pol1, Cauchy.
"""
function fit_single_peak_histogram(h::StatsBase.Histogram, fitrange::Union{Tuple{T,T}, Missing} = missing; fit_function = :Gauss, weights = missing ) where T<:Real
    first_bin, last_bin = ismissing(fitrange) ? (1, length(h.weights)) : (StatsBase.binindex(h, fitrange[1]), StatsBase.binindex(h, fitrange[2]))
    if first_bin < 1 first_bin = 1 end
    if (last_bin > length(h.weights)) last_bin = length(h.weights) end
    hist_weights = h.weights[first_bin:last_bin]
    bin_centers = StatsBase.midpoints(h.edges[1])[first_bin:last_bin]
    bin_widths = map(i -> StatsBase.binvolume(h, i), first_bin:last_bin)

    p0_mean, p0_sigma = mean_and_std(bin_centers, FrequencyWeights(hist_weights), corrected=true)
    p0_scale = sum(h.weights[StatsBase.binindex(h, p0_mean-p0_sigma):StatsBase.binindex(h, p0_mean+p0_sigma)]) / 0.68

    p0, fit_obj = if fit_function in [:Gauss, :Cauchy]
            [p0_scale, p0_sigma, p0_mean], FitFunction(fit_function)
        elseif fit_function in [:Gauss_pol1]
            p0_bg_offset = (hist_weights[1]/bin_widths[1] + hist_weights[end]/bin_widths[end]) / 2
            p0_bg_slope = (hist_weights[1]/bin_widths[1] - hist_weights[end]/bin_widths[end]) / (fitrange[2] - fitrange[1])
            [p0_scale, p0_sigma, p0_mean, p0_bg_offset, p0_bg_slope], FitFunction(:GaussPlusLinearBackground)
    end
    set_fitranges!(fit_obj,([bin_centers[1], bin_centers[end]],))
    set_initial_parameters!(fit_obj, p0)
    lsqfit!(fit_obj, h, weights = weights)
    fit_obj
end

function fit_single_peak_histogram_w_gaussian(h::StatsBase.Histogram, fitrange = missing; weights = missing) where T<:Real
    fit_single_peak_histogram(h, fitrange, fit_function = :Gauss, weights = weights)
end

function fit_single_peak_histogram_w_cauchy(h::StatsBase.Histogram, fitrange = missing; weights = missing) where T<:Real
    fit_single_peak_histogram(h, fitrange, fit_function = :Cauchy, weights = weights)
end

"""
    fit_single_peak_histogram_refined(h::StatsBase.Histogram, fitrange::Union{Tuple{T,T}, Missing} = missing; n_sig = 3.0, fit_function = :Gauss, weights = missing )::FitFunction{T, 1, NP} where {T <: AbstractFloat, NP}

This method is meant for quick and stable fitting. The estimation of initial parameters is done automatically, this allows and easy fit with just one line of code. The user is expected to pass a histogram containing only one peak or specify corresponding fitrange.
The method performs a least square fit with the model specified under  `fit_function = ...`. Currently supported are :Gauss (default), :Gauss_pol1, Cauchy. It then fits again using the results of the first fit as initial parameters.
The fitrange for the second fit can be customized using the `n_sig = ...` keyword.
"""
function fit_single_peak_histogram_refined(h::StatsBase.Histogram, fitrange::Union{Tuple{T,T}, Missing} = missing; n_sig = 4.0, fit_function::Symbol = :Gauss, weights = missing ) where T<:Real
    @assert fit_function in [:Gauss, :Gauss_pol1, :Cauchy] "Please specify 'fit_function = ' either :Gauss, :Gauss_pol1 or :Cauchy as keyword argument. Initial Parameter Estimation is attempted for you."
    fit_obj = fit_single_peak_histogram(h, fitrange, fit_function = fit_function, weights = weights)
    n_sig = minimum([abs(fit_obj.fitted_parameters[3] - fit_obj.bin_centers[1]) / abs(fit_obj.fitted_parameters[2]) , abs(fit_obj.bin_centers[end] - fit_obj.fitted_parameters[3]) / abs(fit_obj.fitted_parameters[2]), n_sig])
    set_fitranges!(fit_obj,([fit_obj.fitted_parameters[3] - n_sig*abs(fit_obj.fitted_parameters[2]), fit_obj.fitted_parameters[3] + n_sig*abs(fit_obj.fitted_parameters[2])],))
    set_initial_parameters!(fit_obj, fit_obj.fitted_parameters)
    lsqfit!(fit_obj, h, weights = weights)
    fit_obj
end

function fit_single_peak_histogram_w_gaussian_refined(h::StatsBase.Histogram, fitrange = missing; weights = missing) where T<:Real
    fit_single_peak_histogram_refined(h, fitrange, fit_function = :Gauss, weights = weights)
end

function fit_single_peak_histogram_w_cauchy_refined(h::StatsBase.Histogram, fitrange = missing; weights = missing) where T<:Real
    fit_single_peak_histogram_refined(h, fitrange, fit_function = :Cauchy, weights = weights)
end


text:
function generate_∇²(d::Domain, T = Float64)
    function L!(y::AbstractVector, x::AbstractVector)
        Δx² = d.Δx^2
        Δy² = d.Δy^2
        c = CartesianIndices(d)
        l = LinearIndices(d)
        @. y = 0
        # ∂x² (Dirichlet BC)
        for I = c[1,:,:]
            y[l[I]] += (-2x[l[I]] + x[l[I+dx]])/Δx²
        end
        for I = c[2:end-1,:,:]
            y[l[I]] += (x[l[I-dx]] - 2x[l[I]] + x[l[I+dx]])/Δx²
        end
        for I = c[end,:,:]
            y[l[I]] += (x[l[I-dx]] - 2x[l[I]])/Δx²
        end
        # + ∂y² (Dirichlet BC)
        for I = c[:,1,:]
            y[l[I]] += (-2x[l[I]] + x[l[I+dy]])/Δy²
        end
        for I = c[:,2:end-1,:]
            y[l[I]] += (x[l[I-dy]] - 2x[l[I]] + x[l[I+dy]])/Δy²
        end
        for I = c[:,end,:]
            y[l[I]] += (x[l[I-dy]] - 2x[l[I]])/Δy²
        end
    end
    return LinearMap{T}(L!, prod(size(d)); ismutating = true)
end

function generate_linearized_L(ψ0, ϕ0, d::Domain; T = Float64)

    function L!(y::AbstractVector, x::AbstractVector)

        # Set up indices and views
        Δx = d.Δx
        Δy = d.Δy
        Δz = d.Δz
        c = CartesianIndices(d)
        l = LinearIndices(d)
        oc = last(c)
        ol = last(l)
        ϕ′ = @view x[1:ol]
        ψ′ = @view x[ol+1:end]
        q′ = @view y[1:ol]
        nil = @view y[ol+1:end]

        # Compute q′
        @. q′ = 0
        # (1 + ∇²ψ0)*∂z²ϕ′ (Neumann BC)
        for I = c[:,:,1]
            q′[l[I]] += (1 + ∂x²(I,ψ0,d) + ∂y²(I,ψ0,d))*(-ϕ′[l[I]] + ϕ′[l[I+dz]])/Δz^2
        end
        for I = c[:,:,2:end-1]
            q′[l[I]] += (1 + ∂x²(I,ψ0,d) + ∂y²(I,ψ0,d))*(ϕ′[l[I-dz]] - 2ϕ′[l[I]] + ϕ′[l[I+dz]])/Δz^2
        end
        for I = c[:,:,end]
            q′[l[I]] += (1 + ∂x²(I,ψ0,d) + ∂y²(I,ψ0,d))*(ϕ′[l[I-dz]] - ϕ′[l[I]])/Δz^2
        end
        # ∂z²ϕ0*∂x²ψ′ (Dirichlet BC)
        for I = c[1,:,:]
            q′[l[I]] += ∂z²(I,ϕ0,d)*(-2ψ′[l[I]] + ψ′[l[I+dx]])/Δx^2
        end
        for I = c[2:end-1,:,:]
            q′[l[I]] += ∂z²(I,ϕ0,d)*(ψ′[l[I-dx]] - 2ψ′[l[I]] + ψ′[l[I+dx]])/Δx^2
        end
        for I = c[end,:,:]
            q′[l[I]] += ∂z²(I,ϕ0,d)*(ψ′[l[I-dx]] - 2ψ′[l[I]])/Δx^2
        end
        # ∂z²ϕ0*∂y²ψ′ (Dirichlet BC)
        for I = c[:,1,:]
            q′[l[I]] += ∂z²(I,ϕ0,d)*(-2ψ′[l[I]] + ψ′[l[I+dy]])/Δy^2
        end
        for I = c[:,2:end-1,:]
            q′[l[I]] += ∂z²(I,ϕ0,d)*(ψ′[l[I-dy]] - 2ψ′[l[I]] + ψ′[l[I+dy]])/Δy^2
        end
        for I = c[:,end,:]
            q′[l[I]] += ∂z²(I,ϕ0,d)*(ψ′[l[I-dy]] - 2ψ′[l[I]])/Δy^2
        end
        # -∂xzψ0*∂xzϕ′ (Neumann on z, Dirichlet on x)
        for I = c[1,:,1]         # x bottom, z bottom
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I+dx+dz]] - ϕ′[l[I+dx]])/(4*Δx*Δz)
        end
        for I = c[2:end-1,:,1]   # x interior, z bottom
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I+dx+dz]] + ϕ′[l[I-dx]] - ϕ′[l[I-dx+dz]] - ϕ′[l[I+dx]])/(4*Δx*Δz)
        end
        for I = c[end,:,1]       # x top, z bottom
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I-dx]] - ϕ′[l[I-dx+dz]])/(4*Δx*Δz)
        end
        for I = c[1,:,2:end-1]     # x bottom, z interior
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I+dx+dz]] - ϕ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[2:end-1,:,2:end-1] # x interior, z interior
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I+dx+dz]] + ϕ′[l[I-dx-dz]] - ϕ′[l[I-dx+dz]] - ϕ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[end,:,2:end-1]   # x top, z interior
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I-dx-dz]] - ϕ′[l[I-dx+dz]])/(4*Δx*Δz)
        end
        for I = c[1,:,end]       # x bottom, z top
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I+dx]] - ϕ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[2:end-1,:,end]   # x interior, z top
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I+dx]] + ϕ′[l[I-dx-dz]] - ϕ′[l[I-dx]] - ϕ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[end,:,end]     # x top, z top
            q′[l[I]] += -∂xz(I,ψ0,d)*(ϕ′[l[I-dx-dz]] - ϕ′[l[I-dx]])/(4*Δx*Δz)
        end
        # -∂xzϕ0*∂xzψ′ (Neumann on z, Dirichlet on x)
        for I = c[1,:,1]         # x bottom, z bottom
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I+dx+dz]] - ψ′[l[I+dx]])/(4*Δx*Δz)
        end
        for I = c[2:end-1,:,1]   # x interior, z bottom
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I+dx+dz]] + ψ′[l[I-dx]] - ψ′[l[I-dx+dz]] - ψ′[l[I+dx]])/(4*Δx*Δz)
        end
        for I = c[end,:,1]       # x top, z bottom
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I-dx]] - ψ′[l[I-dx+dz]])/(4*Δx*Δz)
        end
        for I = c[1,:,2:end-1]     # x bottom, z interior
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I+dx+dz]] - ψ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[2:end-1,:,2:end-1] # x interior, z interior
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I+dx+dz]] + ψ′[l[I-dx-dz]] - ψ′[l[I-dx+dz]] - ψ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[end,:,2:end-1]   # x top, z interior
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I-dx-dz]] - ψ′[l[I-dx+dz]])/(4*Δx*Δz)
        end
        for I = c[1,:,end]       # x bottom, z top
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I+dx]] - ψ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[2:end-1,:,end]   # x interior, z top
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I+dx]] + ψ′[l[I-dx-dz]] - ψ′[l[I-dx]] - ψ′[l[I+dx-dz]])/(4*Δx*Δz)
        end
        for I = c[end,:,end]     # x top, z top
            q′[l[I]] += -∂xz(I,ϕ0,d)*(ψ′[l[I-dx-dz]] - ψ′[l[I-dx]])/(4*Δx*Δz)
        end
        # -∂yzψ0*∂yzϕ′ (Neumann on z, Dirichlet on y)
        for I = c[:,1,1]         # y bottom, z bottom
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I+dy+dz]] - ϕ′[l[I+dy]])/(4*Δy*Δz)
        end
        for I = c[:,2:end-1,1]   # y interior, z bottom
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I+dy+dz]] + ϕ′[l[I-dy]] - ϕ′[l[I-dy+dz]] - ϕ′[l[I+dy]])/(4*Δy*Δz)
        end
        for I = c[:,end,1]       # y top, z bottom
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I-dy]] - ϕ′[l[I-dy+dz]])/(4*Δy*Δz)
        end
        for I = c[:,1,2:end-1]     # y bottom, z interior
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I+dy+dz]] - ϕ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,2:end-1,2:end-1] # y interior, z interior
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I+dy+dz]] + ϕ′[l[I-dy-dz]] - ϕ′[l[I-dy+dz]] - ϕ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,end,2:end-1]   # y top, z interior
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I-dy-dz]] - ϕ′[l[I-dy+dz]])/(4*Δy*Δz)
        end
        for I = c[:,1,end]       # y bottom, z top
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I+dy]] - ϕ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,2:end-1,end]   # y interior, z top
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I+dy]] + ϕ′[l[I-dy-dz]] - ϕ′[l[I-dy]] - ϕ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,end,end]     # y top, z top
            q′[l[I]] += -∂yz(I,ψ0,d)*(ϕ′[l[I-dy-dz]] - ϕ′[l[I-dy]])/(4*Δy*Δz)
        end
        # -∂yzϕ0*∂yzψ′ (Neumann on z, Dirichlet on y)
        for I = c[:,1,1]         # y bottom, z bottom
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I+dy+dz]] - ψ′[l[I+dy]])/(4*Δy*Δz)
        end
        for I = c[:,2:end-1,1]   # y interior, z bottom
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I+dy+dz]] + ψ′[l[I-dy]] - ψ′[l[I-dy+dz]] - ψ′[l[I+dy]])/(4*Δy*Δz)
        end
        for I = c[:,end,1]       # y top, z bottom
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I-dy]] - ψ′[l[I-dy+dz]])/(4*Δy*Δz)
        end
        for I = c[:,1,2:end-1]     # y bottom, z interior
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I+dy+dz]] - ψ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,2:end-1,2:end-1] # y interior, z interior
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I+dy+dz]] + ψ′[l[I-dy-dz]] - ψ′[l[I-dy+dz]] - ψ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,end,2:end-1]   # y top, z interior
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I-dy-dz]] - ψ′[l[I-dy+dz]])/(4*Δy*Δz)
        end
        for I = c[:,1,end]       # y bottom, z top
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I+dy]] - ψ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,2:end-1,end]   # y interior, z top
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I+dy]] + ψ′[l[I-dy-dz]] - ψ′[l[I-dy]] - ψ′[l[I+dy-dz]])/(4*Δy*Δz)
        end
        for I = c[:,end,end]     # y top, z top
            q′[l[I]] += -∂yz(I,ϕ0,d)*(ψ′[l[I-dy-dz]] - ψ′[l[I-dy]])/(4*Δy*Δz)
        end

        # Compute residual of nonlinear balance equation
        @. nil = 0
        # ∂x²ϕ′ - (1 + 2∂y²ψ0)∂x²ψ′ (Dirichlet BC)
        for I = c[1,:,:]
            nil[l[I]] += (
                (-2ϕ′[l[I]] + ϕ′[l[I+dx]])
              - (1 + 2*∂y²(I,ψ0,d))*(-2ψ′[l[I]] + ψ′[l[I+dx]])
            )/Δx^2
        end
        for I = c[2:end-1,:,:]
            nil[l[I]] += (
                (ϕ′[l[I-dx]] - 2ϕ′[l[I]] + ϕ′[l[I+dx]])
              - (1 + 2*∂y²(I,ψ0,d))*(ψ′[l[I-dx]] - 2ψ′[l[I]] + ψ′[l[I+dx]])
            )/Δx^2
        end
        for I = c[end,:,:]
            nil[l[I]] += (
                (ϕ′[l[I-dx]] - 2ϕ′[l[I]])
              - (1 + 2*∂y²(I,ψ0,d))*(ψ′[l[I-dx]] - 2ψ′[l[I]])
            )/Δx^2
        end
        # ∂y²ϕ′ - (1 + 2∂x²ψ0)∂y²ψ′ (Dirichlet BC)
        for I = c[:,1,:]
            nil[l[I]] += (
                (-2ϕ′[l[I]] + ϕ′[l[I+dy]])
              - (1 + 2*∂x²(I,ψ0,d))*(-2ψ′[l[I]] + ψ′[l[I+dy]])
            )/Δy^2
        end
        for I = c[:,2:end-1,:]
            nil[l[I]] += (
                (ϕ′[l[I-dy]] - 2ϕ′[l[I]] + ϕ′[l[I+dy]])
              - (1 + 2*∂x²(I,ψ0,d))*(ψ′[l[I-dy]] - 2ψ′[l[I]] + ψ′[l[I+dy]])
            )/Δy^2
        end
        for I = c[:,end,:]
            nil[l[I]] += (
                (ϕ′[l[I-dy]] - 2ϕ′[l[I]])
              - (1 + 2*∂x²(I,ψ0,d))*(ψ′[l[I-dy]] - 2ψ′[l[I]])
            )/Δy^2
        end
        # 4∂xyψ0*∂xyψ′ (Dirichlet BC)
        for I = c[1,1,:]             # x bottom, y bottom
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(ψ′[l[I+dx+dy]])/(4*Δx*Δy)
        end
        for I = c[2:end-1,1,:]       # x interior, y bottom
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(ψ′[l[I+dx+dy]] - ψ′[l[I-dx+dy]])/(4*Δx*Δy)
        end
        for I = c[end,1,:]           # x top, y bottom
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(-ψ′[l[I-dx+dy]])/(4*Δx*Δy)
        end
        for I = c[1,2:end-1,:]       # x bottom, y interior
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(ψ′[l[I+dx+dy]] - ψ′[l[I+dx-dy]])/(4*Δx*Δy)
        end
        for I = c[2:end-1,2:end-1,:] # x interior, y interior
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(ψ′[l[I+dx+dy]] + ψ′[l[I-dx-dy]] - ψ′[l[I-dx+dy]] - ψ′[l[I+dx-dy]])/(4*Δx*Δy)
        end
        for I = c[end,2:end-1,:]     # x top, y interior
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(ψ′[l[I-dx-dy]] - ψ′[l[I-dx+dy]])/(4*Δx*Δy)
        end
        for I = c[1,end,:]           # x bottom, y top
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(-ψ′[l[I+dx-dy]])/(4*Δx*Δy)
        end
        for I = c[2:end-1,end,:]     # x interior, y top
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(ψ′[l[I-dx-dy]] - ψ′[l[I+dx-dy]])/(4*Δx*Δy)
        end
        for I = c[end,end,:]         # x top, y top
            nil[l[I]] += 4*∂xy(I,ψ0,d)*(ψ′[l[I-dx-dy]])/(4*Δx*Δy)
        end

    end

    return LinearMap{T}(L!, 2*prod(size(d)); ismutating = true)
end

function generate_Lψ(ϕ, d::Domain; T = Float64)

    function L!(y::AbstractVector, x::AbstractVector)
        Δx² = d.Δx^2
        Δy² = d.Δy^2
        c = CartesianIndices(d)
        l = LinearIndices(d)
        # ∂x² (Dirichlet BC)
        for I = c[1,:,:]
            y[l[I]] = (-2x[l[I]] + x[l[I+dx]])/Δx²
        end
        for I = c[2:end-1,:,:]
            y[l[I]] = (x[l[I-dx]] - 2x[l[I]] + x[l[I+dx]])/Δx²
        end
        for I = c[end,:,:]
            y[l[I]] = (x[l[I-dx]] - 2x[l[I]])/Δx²
        end
        # + ∂y² (Dirichlet BC)
        for I = c[:,1,:]
            y[l[I]] += (-2x[l[I]] + x[l[I+dy]])/Δy²
        end
        for I = c[:,2:end-1,:]
            y[l[I]] += (x[l[I-dy]] - 2x[l[I]] + x[l[I+dy]])/Δy²
        end
        for I = c[:,end,:]
            y[l[I]] += (x[l[I-dy]] - 2x[l[I]])/Δy²
        end
        # * (1 + ∂z²ϕ)
        for I = c
            y[l[I]] *= (1 + ∂z²(I,ϕ,d))
        end
    end

    return LinearMap{T}(L!, prod(size(d)); ismutating = true)
end

function generate_Lϕ(ψ, d::Domain; T = Float64)

    function L!(y::AbstractVector, x::AbstractVector)
        Δx² = d.Δx^2
        Δy² = d.Δy^2
        Δz² = d.Δz^2
        c = CartesianIndices(d)
        l = LinearIndices(d)
        # ∂z² (Neumann BC)
        for I = c[:,:,1]
            y[l[I]] = (-x[l[I]] + x[l[I+dz]])/Δz²
        end
        for I = c[:,:,2:end-1]
            y[l[I]] = (x[l[I-dz]] - 2x[l[I]] + x[l[I+dz]])/Δz²
        end
        for I = c[:,:,end]
            y[l[I]] = (x[l[I-dz]] - x[l[I]])/Δz²
        end
        # * (1 + ∇²ψ)
        for I = c
            y[l[I]] *= (1 + ∂x²(I,ψ,d) + ∂y²(I,ψ,d))
        end
        # + ∂x² (Dirichlet BC)
        for I = c[1,:,:]
            y[l[I]] += (-2x[l[I]] + x[l[I+dx]])/Δx²
        end
        for I = c[2:end-1,:,:]
            y[l[I]] += (x[l[I-dx]] - 2x[l[I]] + x[l[I+dx]])/Δx²
        end
        for I = c[end,:,:]
            y[l[I]] += (x[l[I-dx]] - 2x[l[I]])/Δx²
        end
        # + ∂y² (Dirichlet BC)
        for I = c[:,1,:]
            y[l[I]] += (-2x[l[I]] + x[l[I+dy]])/Δy²
        end
        for I = c[:,2:end-1,:]
            y[l[I]] += (x[l[I-dy]] - 2x[l[I]] + x[l[I+dy]])/Δy²
        end
        for I = c[:,end,:]
            y[l[I]] += (x[l[I-dy]] - 2x[l[I]])/Δy²
        end
    end

    return LinearMap{T}(L!, prod(size(d)); ismutating = true)
end

∂x(I::CartesianIndex{3},f,d::Domain) = (f[I+dx] - f[I-dx])/(2d.Δx)
∂y(I::CartesianIndex{3},f,d::Domain) = (f[I+dy] - f[I-dy])/(2d.Δy)
∂z(I::CartesianIndex{3},f,d::Domain) = (f[I+dz] - f[I-dz])/(2d.Δz)
∂x²(I::CartesianIndex{3},f,d::Domain) = (f[I-dx] - 2f[I] + f[I+dx])/d.Δx^2
∂y²(I::CartesianIndex{3},f,d::Domain) = (f[I-dy] - 2f[I] + f[I+dy])/d.Δy^2
∂z²(I::CartesianIndex{3},f,d::Domain) = (f[I-dz] - 2f[I] + f[I+dz])/d.Δz^2
∂xy(I::CartesianIndex{3},f,d::Domain) = (
    f[I+dx+dy] + f[I-dx-dy] - f[I+dx-dy] - f[I-dx+dy]
)/(4*d.Δx*d.Δy)
∂xz(I::CartesianIndex{3},f,d::Domain) = (
    f[I+dx+dz] + f[I-dx-dz] - f[I+dx-dz] - f[I-dx+dz]
)/(4*d.Δx*d.Δz)
∂yz(I::CartesianIndex{3},f,d::Domain) = (
    f[I+dy+dz] + f[I-dy-dz] - f[I+dy-dz] - f[I-dy+dz]
)/(4*d.Δy*d.Δz)


text:
module TrueRandom

using Random
import Random: rand, rand!, _rand52

import Base.show

using Base.Threads

using HTTP
using JSON3

include("UInt16RNG.jl")

export ANUQRNG

include("ANUQRNG.jl")

end


text:
using Documenter
using ExtensibleScheduler

makedocs(
    format = :html,
    sitename = "ExtensibleScheduler.jl",
    pages = [
        "index.md",
        "getting_started.md",
    ]
)

deploydocs(
    repo = "github.com/scls19fr/ExtensibleScheduler.jl.git",
)


text:
type TraversalControls
    skip_children :: Bool
    call_after :: Bool
    called_after :: Bool
    args :: Any
    TraversalControls(args) = new(false, false, false, args)
end

function skip_children!(controls)
    controls.skip_children = true
end

function call_after!(controls)
    controls.call_after = true
end

function replace_args!(controls, args...)
    controls.args = args
end

function inspect_and_modify(func::Function, state, val, args...)

    controls = TraversalControls(args)

    state, new_val = func(controls, state, val, args...)

    children_args = controls.args

    if !controls.skip_children && new_val === val
        subnodes = to_subnodes(val)
        if !(subnodes === nothing)
            new_subnodes = Array{Any, 1}(length(subnodes))
            changed = false
            for (i, subnode) in enumerate(subnodes)
                state, new_subnode = inspect_and_modify(func, state, subnode, children_args...)
                changed |= !(new_subnode === subnode)
                new_subnodes[i] = new_subnode
            end

            new_val = changed ? from_subnodes(val, new_subnodes) : val
        end
    end

    if controls.call_after
        controls.called_after = true
        state, new_val = func(controls, state, new_val, args...)
    end

    state, new_val
end

function inspect(func::Function, state, val, args...)
    state, _ = inspect_and_modify(state, val, args...) do controls, state, val, args...
        func(controls, state, val, args...), val
    end
    state
end

function modify(func::Function, val, args...)
    _, new_val = inspect_and_modify(nothing, val, args...) do controls, state, val, args...
        state, func(controls, val, args...)
    end
    new_val
end


text:
function iterate!(job::BasicMCJob, ::Type{NUTS}, ::Type{Univariate})
  local a::Real
  local na::Integer

  if isa(job.tuner, DualAveragingMCTuner)
    job.sstate.count += 1
  end

  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  job.sstate.momentum = randn()

  job.sstate.oldhamiltonian = hamiltonian(job.pstate.logtarget, job.sstate.momentum)

  job.sstate.pstateplus.value = job.pstate.value
  job.sstate.pstateplus.gradlogtarget = job.pstate.gradlogtarget
  job.sstate.momentumplus = job.sstate.momentum
  job.sstate.pstateminus.value = job.pstate.value
  job.sstate.pstateminus.gradlogtarget = job.pstate.gradlogtarget
  job.sstate.momentumminus = job.sstate.momentum

  job.sstate.j = 0
  job.sstate.n = 1
  job.sstate.s = true
  if in(:accept, job.outopts[:diagnostics]) || job.tuner.verbose
    job.sstate.update = false
  end

  job.sstate.u = log(rand())+job.sstate.oldhamiltonian

  while (job.sstate.s && (job.sstate.j < job.sampler.maxndoublings))
    job.sstate.v = rand(Bool) ? 1 : -1

    if isa(job.tuner, VanillaMCTuner)
      if job.sstate.v == -1
        job.sstate.pstateminus,
        job.sstate.momentumminus,
        _,
        _,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime = build_tree!(
          job.sstate,
          job.sstate.pstateminus,
          job.sstate.momentumminus,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          VanillaMCTuner,
          job.outopts
        )
      else
        _,
        _,
        job.sstate.pstateplus,
        job.sstate.momentumplus,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime = build_tree!(
          job.sstate,
          job.sstate.pstateplus,
          job.sstate.momentumplus,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          VanillaMCTuner,
          job.outopts
        )
      end
    elseif isa(job.tuner, DualAveragingMCTuner)
      if job.sstate.v == -1
        job.sstate.pstateminus,
        job.sstate.momentumminus,
        _,
        _,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime,
        a,
        na = build_tree!(
          job.sstate,
          job.sstate.pstateminus,
          job.sstate.momentumminus,
          job.sstate.oldhamiltonian,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          DualAveragingMCTuner,
          job.outopts
        )
      else
        _,
        _,
        job.sstate.pstateplus,
        job.sstate.momentumplus,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime,
        a,
        na = build_tree!(
          job.sstate,
          job.sstate.pstateplus,
          job.sstate.momentumplus,
          job.sstate.oldhamiltonian,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          DualAveragingMCTuner,
          job.outopts
        )
      end
    end

    if (job.sstate.sprime && (rand() < job.sstate.nprime/job.sstate.n))
      job.pstate.value = job.sstate.pstateprime.value

      job.pstate.gradlogtarget = job.sstate.pstateprime.gradlogtarget
      if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
        job.pstate.gradloglikelihood = job.sstate.pstateprime.gradloglikelihood
      end
      if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
        job.pstate.gradlogprior = job.sstate.pstateprime.gradlogprior
      end

      job.pstate.logtarget = job.sstate.pstateprime.logtarget
      if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
        job.pstate.loglikelihood = job.sstate.pstateprime.loglikelihood
      end
      if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
        job.pstate.logprior = job.sstate.pstateprime.logprior
      end

      if in(:accept, job.outopts[:diagnostics]) || job.tuner.verbose
        job.sstate.update = true
      end
    end

    job.sstate.j += 1
    job.sstate.n += job.sstate.nprime
    job.sstate.s =
      job.sstate.sprime &&
      !uturn(job.sstate.pstateplus.value, job.sstate.pstateminus.value, job.sstate.momentumplus, job.sstate.momentumminus)
  end

  if !isempty(job.sstate.diagnosticindices)
    if haskey(job.sstate.diagnosticindices, :accept)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = job.sstate.update
    end

    if haskey(job.sstate.diagnosticindices, :ndoublings)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:ndoublings]] = job.sstate.j
    end

    if haskey(job.sstate.diagnosticindices, :a)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:a]] = a
    end

    if haskey(job.sstate.diagnosticindices, :na)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:na]] = na
    end
  end

  if job.tuner.verbose
    job.sstate.update && (job.sstate.tune.accepted += 1)
  end

  if isa(job.tuner, VanillaMCTuner) && job.tuner.verbose
    if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
      rate!(job.sstate.tune)

      if job.tuner.verbose
        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )
      end

      reset_burnin!(job.sstate.tune)
    end
  elseif isa(job.tuner, DualAveragingMCTuner)
    if (job.sstate.count <= job.tuner.nadapt)
      tune!(job.sstate.tune, job.tuner, job.sstate.count, a/na)
    else
      job.sstate.tune.step = job.sstate.tune.εbar
    end

    if job.tuner.verbose
      if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
        rate!(job.sstate.tune)

        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )

        reset_burnin!(job.sstate.tune)
      end
    end
  end
end

function iterate!(job::BasicMCJob, ::Type{NUTS}, ::Type{Multivariate})
  local a::Real
  local na::Integer

  if isa(job.tuner, DualAveragingMCTuner)
    job.sstate.count += 1
  end

  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  job.sstate.momentum[:] = randn(job.pstate.size)

  job.sstate.oldhamiltonian = hamiltonian(job.pstate.logtarget, job.sstate.momentum)

  job.sstate.pstateplus.value = copy(job.pstate.value)
  job.sstate.pstateplus.gradlogtarget = copy(job.pstate.gradlogtarget)
  job.sstate.momentumplus = copy(job.sstate.momentum)
  job.sstate.pstateminus.value = copy(job.pstate.value)
  job.sstate.pstateminus.gradlogtarget = copy(job.pstate.gradlogtarget)
  job.sstate.momentumminus = copy(job.sstate.momentum)

  job.sstate.j = 0
  job.sstate.n = 1
  job.sstate.s = true
  if in(:accept, job.outopts[:diagnostics]) || job.tuner.verbose
    job.sstate.update = false
  end

  job.sstate.u = log(rand())+job.sstate.oldhamiltonian

  while (job.sstate.s && (job.sstate.j < job.sampler.maxndoublings))
    job.sstate.v = rand(Bool) ? 1 : -1

    if isa(job.tuner, VanillaMCTuner)
      if job.sstate.v == -1
        job.sstate.pstateminus,
        job.sstate.momentumminus,
        _,
        _,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime = build_tree!(
          job.sstate,
          job.sstate.pstateminus,
          job.sstate.momentumminus,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          VanillaMCTuner,
          job.outopts
        )
      else
        _,
        _,
        job.sstate.pstateplus,
        job.sstate.momentumplus,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime = build_tree!(
          job.sstate,
          job.sstate.pstateplus,
          job.sstate.momentumplus,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          VanillaMCTuner,
          job.outopts
        )
      end
    elseif isa(job.tuner, DualAveragingMCTuner)
      if job.sstate.v == -1
        job.sstate.pstateminus,
        job.sstate.momentumminus,
        _,
        _,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime,
        a,
        na = build_tree!(
          job.sstate,
          job.sstate.pstateminus,
          job.sstate.momentumminus,
          job.sstate.oldhamiltonian,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          DualAveragingMCTuner,
          job.outopts
        )
      else
        _,
        _,
        job.sstate.pstateplus,
        job.sstate.momentumplus,
        job.sstate.pstateprime,
        job.sstate.nprime,
        job.sstate.sprime,
        a,
        na = build_tree!(
          job.sstate,
          job.sstate.pstateplus,
          job.sstate.momentumplus,
          job.sstate.oldhamiltonian,
          job.sstate.u,
          job.sstate.v,
          job.sstate.j,
          job.sstate.tune.step,
          job.parameter,
          job.sampler.maxδ,
          NUTS,
          DualAveragingMCTuner,
          job.outopts
        )
      end
    end

    if (job.sstate.sprime && (rand() < job.sstate.nprime/job.sstate.n))
      job.pstate.value = copy(job.sstate.pstateprime.value)

      job.pstate.gradlogtarget = copy(job.sstate.pstateprime.gradlogtarget)
      if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
        job.pstate.gradloglikelihood = copy(job.sstate.pstateprime.gradloglikelihood)
      end
      if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
        job.pstate.gradlogprior = copy(job.sstate.pstateprime.gradlogprior)
      end

      job.pstate.logtarget = job.sstate.pstateprime.logtarget
      if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
        job.pstate.loglikelihood = job.sstate.pstateprime.loglikelihood
      end
      if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
        job.pstate.logprior = job.sstate.pstateprime.logprior
      end

      if in(:accept, job.outopts[:diagnostics]) || job.tuner.verbose
        job.sstate.update = true
      end
    end

    job.sstate.j += 1
    job.sstate.n += job.sstate.nprime
    job.sstate.s =
      job.sstate.sprime &&
      !uturn(job.sstate.pstateplus.value, job.sstate.pstateminus.value, job.sstate.momentumplus, job.sstate.momentumminus)
  end

  if !isempty(job.sstate.diagnosticindices)
    if haskey(job.sstate.diagnosticindices, :accept)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = job.sstate.update
    end

    if haskey(job.sstate.diagnosticindices, :ndoublings)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:ndoublings]] = job.sstate.j
    end

    if haskey(job.sstate.diagnosticindices, :a)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:a]] = a
    end

    if haskey(job.sstate.diagnosticindices, :na)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:na]] = na
    end
  end

  if job.tuner.verbose
    job.sstate.update && (job.sstate.tune.accepted += 1)
  end

  if isa(job.tuner, VanillaMCTuner) && job.tuner.verbose
    if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
      rate!(job.sstate.tune)

      if job.tuner.verbose
        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )
      end

      reset_burnin!(job.sstate.tune)
    end
  elseif isa(job.tuner, DualAveragingMCTuner)
    if (job.sstate.count <= job.tuner.nadapt)
      tune!(job.sstate.tune, job.tuner, job.sstate.count, a/na)
    else
      job.sstate.tune.step = job.sstate.tune.εbar
    end

    if job.tuner.verbose
      if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
        rate!(job.sstate.tune)

        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )

        reset_burnin!(job.sstate.tune)
      end
    end
  end
end


text:
using Base.Test
using Permutations
p = Permutation([7:-1:1])
i = Permutation(7)
@test p*p == i
@test p == inv(p)
a =[2,3,1,6,7,8,5,4,9]
p = Permutation(a)
@test order(p) == 6
@test p^6 == Permutation(9)
@test p^-5 == p
@test array(p)==a

q = Permutation([1,5,3,9,4,8,6,7,2] )
@test p*q != q*p
@test q[2] == 5

P = matrix(p)
Q = matrix(q)
@test P*Q == matrix(q*p)
@test P' == matrix(inv(p))


text:
module Cryptoworld

# Write your package code here.

include("huobi.jl")
using .huobi

include("bithumb.jl")
using .bithumb
end


text:
module CytofDiffDensity

using Distributions
import Distributions.StatsBase: ecdf
using MCMC
using MCMC: Bijectors
using Random
using StatsFuns
using StatsPlots
using LaTeXStrings

export Gtilde, Pzero, infer
export MixSkewT, OrderedNormalMeanPrior, update

export SkewT, skewtlogpdf, skewtpdf, randskewt
export AltSkewT, skewfromaltskewt, scalefromaltskewt, fromaltskewt, toaltscale, toaltskew
export make_sampler, print_model_info, make_ordered_prior
export find_good_seed, make_callback
export loglike, simulate_gtilde, simulate_ftilde
export count_small_clusters
export HurdleModel, hellinger

include("hurdle.jl")
include("util.jl")
include("OrderedNormalMeanPrior.jl")
include("skewt.jl")
include("models/Gtilde.jl")
include("models/Gtilde-updates.jl")
include("models/Pzero.jl")
include("models/mixskewt-model.jl")
include("models/mixskewt-updates.jl")
include("simulate.jl")
include("metrics.jl")
include("postprocess.jl")

end # module


text:
# TODO: add kwarg handling for filtering by initial_time and len
function _read_realized_results(
    result_values::Dict{Symbol, DataFrames.DataFrame},
    names::Union{Nothing, Vector{Symbol}},
)
    existing_names = collect(keys(result_values))
    names = isnothing(names) ? existing_names : names
    PSI._validate_names(existing_names, names)
    return deepcopy(filter(p -> (p.first ∈ names), result_values))
end

function _read_results(
    result_values::Dict{Symbol, DataFrames.DataFrame},
    names::Union{Nothing, Vector{Symbol}},
    initial_time::Dates.DateTime,
)
    realized_results = _read_realized_results(result_values, names)
    results = PSI.FieldResultsByTime()
    for (name, df) in realized_results
        results[name] = PSI.ResultsByTime(initial_time => df)
    end
    return results
end

function PSI.read_realized_variables(
    res::PSI.ProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    kwargs...,
)
    variable_values = PSI.get_variables(res)
    return _read_realized_results(variable_values, names)
end

function PSI.read_realized_parameters(
    res::PSI.ProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    kwargs...,
)
    parameter_values = IS.get_parameters(res)
    return _read_realized_results(parameter_values, names)
end

function PSI.read_realized_duals(
    res::PSI.ProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    kwargs...,
)
    dual_values = get_duals(res)
    return _read_realized_results(dual_values, names)
end

function PSI.read_variables(
    res::PSI.ProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    kwargs...,
)
    result_values = PSI.get_variables(res)
    return _read_results(result_values, names, first(PSI.get_timestamps(res)))
end

function PSI.read_parameters(
    res::PSI.ProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    kwargs...,
)
    result_values = IS.get_parameters(res)
    return _read_results(result_values, names, first(PSI.get_timestamps(res)))
end

function read_duals(
    res::PSI.ProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    kwargs...,
)
    result_values = PSI.get_duals(res)
    return _read_results(result_values, names, first(PSI.get_timestamps(res)))
end

PSI.get_realized_timestamps(res::PSI.ProblemResults; kwargs...) = PSI.get_timestamps(res)


text:
module ImageEdgeDetection

using ColorVectorSpace
using ImageCore
using ImageFiltering
using MappedArrays
using Parameters: @with_kw # Same as Base.@kwdef but works on Julia 1.0

# TODO: port EdgeDetectionAPI to ImagesAPI
include("EdgeDetectionAPI/EdgeDetectionAPI.jl")
import .EdgeDetectionAPI: AbstractEdgeDetectionAlgorithm,
                          detect_edges, detect_edges!

 import .EdgeDetectionAPI: AbstractEdgeThinningAlgorithm,
                           thin_edges, thin_edges!
# TODO Relax this to all image color types
const GenericGrayImage = AbstractArray{<:Union{Number, AbstractGray}}

include("algorithms/canny.jl")
include("algorithms/nonmaxima_suppression.jl")

export
    # main types and functions
    Canny,
    NonmaximaSuppression,
    thin_edges,
    thin_edges!,
    detect_edges,
    detect_edges!

end # module


text:
"""
Julia expression wrapper
"""
struct Expression
    expr::Expr
    syms::Dict{Symbol,Int}
end
function Expression(ex::Expr)
    syms = Dict( s=>i for (i,s) in pairs(sort!(symbols(ex))) )
    Expression(ex, syms)
end
show(io::IO, e::Expression) = infix(io, e.expr)
(e::Expression)(vals::T...) where {T} = evaluate(e.expr, e.syms, vals...)

function symbols(ex::Expr)
    syms = Symbol[]
    for e in ex.args
        isa(e, Symbol) && push!(syms, e)
        isa(e, Expr) && append!(syms, symbols(e))
    end
    unique!(syms)
end

function compile(ex::Expr, params::Vector{Symbol})
    tprm = Expr(:tuple, params...)
    Expr(:->, tprm, ex)
end

height(ex) = isa(ex, Expr) ? maximum( height(e) for e in ex.args )+1 : 0
nodes(ex) = !isa(ex, Expr) ? 1 : length(ex.args) > 0 ? sum( nodes(e) for e in ex.args ) : 0
length(ex) = nodes(ex)

function depth(root, ex; d=0)
    return if root == ex
        d
    elseif isa(root, Expr)
        maximum( depth(e, ex, d=d+1) for e in root.args )
    else
        -1
    end
end

function copyto!(dst::Expr, src::Expr)
    dst.head = src.head
    dst.args = deepcopy(src.args)
    dst
end

function randsubexpr(ex)
    !isa(ex, Expr) && return ex
    csize = map(nodes, ex.args[2:end])
    cidx = rand(1:sum(csize)+1)-1
    if cidx == 0
        ex
    else
        nidx = findfirst(i->cidx<=i, cumsum(csize))
        randsubexpr(ex.args[nidx+1])
    end
end

function getindex(ex::Expr, idx::Int)
    if idx == 0
        ex
    else
        csize = cumsum(map(nodes, ex.args[2:end]))
        nidx = findfirst(i->i>=idx, csize)
        nex = ex.args[nidx+1]
        !isa(nex, Expr) ? nex : nex[csize[nidx]-idx]
    end
end

function setindex!(ex::Expr, subex, idx::Int)
    csize = cumsum(map(nodes, ex.args[2:end]))
    nidx = findfirst(i->i>=idx, csize)
    nex = ex.args[nidx+1]
    if !isa(nex, Expr) || csize[nidx]-idx == 0
        ex.args[nidx+1] = subex
    else
        nex[csize[nidx]-idx] = subex
    end
end

isnum(ex) = isa(ex, Number)
iszeronum(root) = isnum(root) && iszero(root)
isonenum(root) = isnum(root) && isone(root)
isdiv(ex) = ex in [/, div, pdiv, aq]
isexpr(ex) = isa(ex, Expr)
issym(ex) = isa(ex, Symbol)
isexprsym(ex) = isexpr(ex) || issym(ex)
isbinexpr(ex) = isexpr(ex) && length(ex.args) == 3

function evaluate(ex::Expr, psyms::Dict{Symbol,Int}, vals::T...)::T where {T}
    exprm = ex.args
    exvals = (isexpr(nex) || issym(nex) ? evaluate(nex, psyms, vals...) : nex for nex in exprm[2:end])
    exprm[1](exvals...)
end

function evaluate(ex::Symbol, psyms::Dict{Symbol,Int}, vals::T...)::T where {T}
    pidx = get(psyms, ex, 0)
    if pidx == 0
        @error "Undefined symbol: $ex"
    end
    return T(vals[pidx])
end

function simplifyunary!(root)
    fn, op = root.args
    if isexpr(op)
        fn2, op2 = op.args
        if fn == log && fn2 == exp
            return op2
        elseif fn == exp && fn2 == log
            return op2
        end
    end
    return root
end

function simplifybinary!(root)
    fn, op1, op2 = root.args

    # some elementary simplification rules
    if fn == (-) && op1 == op2 # additive inverse: x-x = 0
        root = 0
    elseif isdiv(fn) && (op1 == op2 || iszeronum(op2)) # look for multiplicative inverse x/x = 1
        root = 1
    elseif (fn == (*) || isdiv(fn)) &&  (iszeronum(op1) || iszeronum(op2)) # look for 0: x*0 = 0*x = 0
        root = 0
    elseif (fn == (+) || fn == (-)) && iszeronum(op2) # look for 0: x ± 0 = x
        root = op1
    elseif fn == (+) && iszeronum(op1) # look for 0: 0 + x = x
        root = op2
    elseif (fn == (*) || isdiv(fn)) && isonenum(op2) # x*1 = x || x/1 = x
        root = op1
    elseif fn == (*) && isonenum(op1) # 1*x = x
        root = op2
    elseif (fn == (+)) && (op1 == op2) # x+x = 2x
        root.args[1] = (*)
        root.args[2] = 2
    elseif (fn == (+) || fn == (-)) && (isbinexpr(op1) && isnum(op2))
        fn2, op11, op12 = op1.args
        if fn == (+) && fn2 == (-) && isnum(op12)
            # (x-m)+n = x+(n-m)
            root.args[2] = op11
            root.args[3] = op2-op12
        elseif fn == (-) && fn2 == (-) && isnum(op12)
            # (x-m)-n = x-(n+m)
            root.args[2] = op11
            root.args[3] = op12+op2
        elseif fn2 == (+) && (isnum(op11) || isnum(op12))
            # (m+x)±n = (x+m)±n = x+(n±m)
            var, n2 = isnum(op11) ? (op12, op11) : (op11, op12)
            n3 = fn(n2, op2)
            root.args[1] = fn2
            root.args[2] = var
            root.args[3] = n3
        elseif fn2 == (-) && isnum(op11)
            # (m-x)±n = (n±m)-x
            n3 = fn(op11, op2)
            root.args[1] = fn2
            root.args[2] = n3
            root.args[3] = op12
        end
    elseif fn == (+) && isnum(op1) && isbinexpr(op2)
        fn2, op21, op22 = op2.args
        if fn2 == (+) || fn2 == (-)
            if isnum(op21)
                # n+(m±x) = (n+m)±x
                root.args[1] = fn2
                root.args[2] = op1 + op21
                root.args[3] = op22
            elseif isnum(op22)
                # n+(x±m) = (n±m)+x
                root.args[1] = fn
                root.args[2] = fn2(op1, op22)
                root.args[3] = op21
            end
        end
    elseif fn == (-) && isnum(op1) && isbinexpr(op2)
        fn2, op21, op22 = op2.args
        if fn2 == (+) && (isnum(op21) || isnum(op22))
            # n-(x+m) = n-(m+x) = (n-m)-x
            var, n2 = isnum(op21) ? (op22, op21) : (op21, op22)
            root.args[2] = op1 - n2
            root.args[3] = var
        elseif fn2 == (-) && (isnum(op21) || isnum(op22))
            # n-(m-x) = (n-m)+x
            # n-(x-m) = (n+m)-x
            var, n2, f1p = isnum(op21) ? (op22, op21, true) : (op21, op22, false)
            root.args[1] = f1p ? (+) : (-)
            root.args[2] = f1p ? op1 - n2 : op1 + n2
            root.args[3] = var
        end
    elseif fn == (+) && (isexpr(op1) || isexpr(op2))
        if isbinexpr(op2)
            # x + (n - x) = n
            fn2, op21, op22 = op2.args
            if fn2 == (-) && op1 == op22
                root = op21
            end
        elseif isbinexpr(op1)
            # (n - x) + x = n
            fn2, op11, op12 = op1.args
            if fn2 == (-) && op12 == op2
                root = op11
            end
        end
    elseif fn == (-) && (isexpr(op1) || isexpr(op2))
        if isbinexpr(op2)
            fn2, op21, op22 = op2.args
            if op1 == op21
                # x - (x ± n) = -±n
                if fn2 == (-)
                    root = op22
                else
                    pop!(root.args)
                    root.args[end] = op22
                end
            elseif op1 == op22 && fn2 == (+)
                # x - (n + x) = -n
                pop!(root.args)
                root.args[end] = op21
            end
        elseif isbinexpr(op1)
            # (x ± n) - x = ±n
            fn2, op11, op12 = op1.args
            if op11 == op2
                if fn2 == (+)
                    root = op12
                else
                    root.args[1] = fn2
                    root.args[2] = op12
                    pop!(root.args)
                end
            end
        end
    end

    # evaluate numerical arithmetic
    isnum(op1) && isnum(op2) && return fn(op1, op2)

    return root
end

"""
    simplify!(expr)

Simplify an algebraic expression
"""
function simplify!(root)
    # return if not an expression
    !isa(root, Expr) && return root
    # return if not a function
    root.head != :call && return root

    n = length(root.args)
    # recursively simplify arguments
    for i in 2:n
        !isexpr(root.args[i]) && continue
        root.args[i] = simplify!(root.args[i])
    end

    n == 2 && return simplifyunary!(root)
    n == 3 && return simplifybinary!(root)
    return root
end

function infix(io::IO, root; digits=3)
    if isa(root, Number)
        print(io, round(root, digits=digits))
    elseif isa(root, Expr) && root.head == :call
        if root.args[1] ∈ [+, -, *, /, ^]
            print(io, "(")
            infix(io, root.args[2])
            infix(io, root.args[1])
            length(root.args)>2 && infix(io, root.args[3])
            print(io, ")")
        else
            infix(io, root.args[1])
            print(io, "(")
            for (i, ex) in enumerate(root.args[2:end])
                i > 1 && print(io, ", ")
                infix(io, ex)
            end
            print(io, ")")
        end
    else
        print(io, root)
    end
end
show(io::IO, ::MIME"text/html", e::Expression) = infix(io, e.expr)

latex(io::IO, root::Number; digits=3) = print(io, round(root, digits=digits))
latex(io::IO, root::Symbol; kwargs...) = print(io, root)
function latex(io::IO, root::Expr; digits=3)
    root.head != :call && print(io, expr)
    if root.args[1] == (/)
        print(io, "\\frac{")
        latex(io, root.args[2], digits=digits)
        print(io, "}{")
        latex(io, root.args[3], digits=digits)
        print(io, "}")
    elseif root.args[1] ∈ [+, -, *]
        print(io, "\\left(")
        latex(io, root.args[2], digits=digits)
        print(io, root.args[1])
        latex(io, root.args[3], digits=digits)
        print(io, "\\right)")
    elseif root.args[1] == (^)
        print(io, "{")
        latex(io, root.args[2], digits=digits)
        print(io, "}^{")
        latex(io, root.args[3], digits=digits)
        print(io, "}")
    else
        print(io, "\\")
        print(io, root.args[1])
        print(io, "\\left(")
        for ex in root.args[2:end]
            latex(io, ex, digits=digits)
        end
        print(io, "\\right)")
    end
end
show(io::IO, ::MIME"text/latex", e::Expression) = latex(io, e.expr)


text:
using JSON
using JLD
# Nproc=4
#
# rmprocs(workers())
# addprocs(Nproc)

@everywhere using Base.Test
@everywhere import Base.convert
@everywhere using MAT
@everywhere using ForwardDiff
@everywhere using Polynomials
@everywhere import ForwardDiff.DiffBase
@everywhere using Optim
@everywhere gc(true)

# local
@everywhere path="/home/emma/fit_behaviour_julia/"
# cluster
# @everywhere path="/home/hcli64/hcli64751/rats/fit_behaviour_julia/"
#
@everywhere include(path*"functions/functions_julia.jl")
@everywhere include(path*"functions/parameter_functions.jl")
@everywhere include(path*"functions/utils.jl")

rat_num = ARGS[1]
#rat_num = "24"
#dataset = ARGS[2]
dataset ="2"

path_load=path*"scripts_rats/data/"
if dataset == "1"
    filename="processed_data_rat"*rat_num*".json"
else
    filename = "processed_data_rat"*rat_num*"_dataset2_1.json"
end

# @everywhere include(path*"scripts/functions/functions_julia.jl")
# @everywhere include(path*"scripts/functions/parameter_functions.jl")
# @everywhere include(path*"scripts/functions/utils.jl")
@everywhere include(path*"/functions/simulations.jl")

# model_fit="PI_StConst"

# isubject=parse(Int,ARGS[1])
# #isubject=1
# SUBJECTS=[1,4,5,6,7,8,9,10,11,12,14,15,16,17,18,19]
# subject=SUBJECTS[isubject]

#
# path_data=path*"data/"
# data=load(path_data*"data_subject"*string(subject)*".jld")

data = Dict()
open(path_load*filename, "r") do f
    global data
    dicttxt = readstring(f)  # file information to string
    data=JSON.parse(dicttxt)  # parse and transform data
end

n_stim_10 = data["n_stim_10"]
n_stim_20 = data["n_stim_20"]

#Nstim=100
# stim=data["stim"][1:Nstim,:]
# choices=data["choices"][1:Nstim]
# i_sigma=data["i_sigma"][1:Nstim]
stim=data["stim_20"]
stim = permutedims(reshape(stim,(20,n_stim_20)))
choices = data["choices_20"]

# i_sigma=data["i_sigma"]
Nstim=size(stim)[1]

history_bias=zeros(Nstim) #so far no history bias

Nframes=20
Tframe=0.5

model_fit = "StLogPi"
model_f=Dict(:bias_d_f=>bias_d_const,
:hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
:ro_f=>read_out_perfect,:lambda=>0.0)

args=[ ("st",[1,2]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",0)]
consts=[("st",[3,4]),("x0",[1,2]),("bias_d",[1]),
("hbias_d",[1]),("sigma_a",[1])]

z=1.0*[0.0,0.0,0.0,0,0,0.0,1.0]

ub=[50.,50.]
lb=[0.0,-50.]

x_ini=zeros(length(ub))
# model_f=Dict(:bias_d_f=>bias_d_const,:hbias_d_f=>history_bias_d_const,
# :st_f=>st_const,:x0_f=>initial_condition_bias_hbias,
# :ro_f=>read_out_perfect,:lambda=>0.0)
#
#
# args=[ ("st",[1]),("x0",0),("bias_d",0),("sigma_a",0)]
# consts=[("st",0),("x0",[1,2]),("bias_d",[1]),("sigma_a",[1])]
#
# z=1.0*[0.0,0.0,0.0,1.0]
# #
# ub=[10.,10.]
# lb=[0.0,0.0]

for i in 1:300
        aux=2*(rand(8)-0.5)
        x_ini=[abs(aux[1])*10,aux[2]*10]
        #x_ini = [abs(aux[1])]
        param=make_dict(args,x_ini,consts,z)
        #PR=ProbRightPI(stim,history_bias,Tframe,param,model_f)
        #ll=ComputeLLPI(stim,choices,history_bias,Tframe,args,consts,x,z,model_f)
        control=true
        try
                global res=ModelFittingPI(stim,choices,history_bias, Tframe, args,consts,ub,lb, x_ini,z,model_f)
                println("Fit_done")

        catch
                println("something wrong fitting the model icond",i)
                control=false
        end

        if control
                j=1
                path_save_res=path*"results/model_fit"*model_fit*"/rat"*string(rat_num)*"_1sec/"
                mkpath(path_save_res)
                filename_save=path_save_res*"icond"*string(j)*".jld"
                param2=make_dict(args,res["x_bf"],consts,z)

                PR=ProbRightPI(stim,history_bias,Tframe,param2,model_f)

                while isfile(filename_save)
                        j=j+1
                        filename_save=path_save_res*"icond"*string(j)*".jld"
                end

                println(j," ",res["x_bf"])
                save(filename_save,"res",res,"PR_training",PR,"icond",j,"args",args,"consts",consts,"z",z,
                "param",param2)
                println("done")
        end
end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
#Install required packages
Pkg.add("Images")
Pkg.add("DataFrames")
using Images
using DataFrames

#typeData could be either "train" or "test.
#labelsInfo should contain the IDs of each image to be read
#The images in the trainResized and testResized data files
#are 20x20 pixels, so imageSize is set to 400.
#path should be set to the location of the data files.

function read_data(typeData, labelsInfo, imageSize, path)
 #Intialize x matrix
 x = zeros(size(labelsInfo, 1), imageSize)

 for (index, idImage) in enumerate(labelsInfo["ID"])
  #Read image file
  nameFile = "$(path)/$(typeData)Resized/$(idImage).Bmp"
  img = imread(nameFile)

  #Convert img to float values
  temp = float32sc(img)

  #Convert color images to gray images
  #by taking the average of the color scales.
  if ndims(temp) == 3
   temp = mean(temp.data, 1)
  end

  #Transform image matrix to a vector and store
  #it in data matrix
  x[index, :] = reshape(temp, 1, imageSize)
 end
 return x
end

imageSize = 400 # 20 x 20 pixels

#Set location of data files , folders
path = ...

#Read information about training data , IDs.
labelsInfoTrain = readtable("$(path)/trainLabels.csv")

#Read training matrix
xTrain = read_data("train", labelsInfoTrain, imageSize, path)

#Read information about test data ( IDs ).
labelsInfoTest = readtable("$(path)/sampleSubmission.csv")

#Read test matrix
xTest = read_data("test", labelsInfoTest, imageSize, path)

#Get only first character of string (convert from string to character).
#Apply the function to each element of the column "Class"
yTrain = map(x -> x[1], labelsInfoTrain["Class"])

#Convert from character to integer
yTrain = int(yTrain)

Pkg.add("DecisionTree")
using DecisionTree

#Train random forest with
#20 for number of features chosen at each random split,
#50 for number of trees,
#and 1.0 for ratio of subsampling.
model = build_forest(yTrain, xTrain, 20, 50, 1.0)

#Get predictions for test data
predTest = apply_forest(model, xTest)

#Convert integer predictions to character
labelsInfoTest["Class"] = char(predTest)

#Save predictions
writetable("$(path)/juliaSubmission.csv", labelsInfoTest, separator=',', header=true)

#Run 4 fold cross validation
accuracy = nfoldCV_forest(yTrain, xTrain, 20, 50, 4, 1.0);
println ("4 fold accuracy: $(mean(accuracy))")


text:
#import StatsBase: mean, var, std, sem, mean_and_std, mean_and_var
mean_and_sem(x) = (m = mean(x); s = sem(x); (m, s))

function check_l1loss(est)
    if isa(loss(est), L1Loss)
        mess = "variance is not well defined with L1Loss, use another loss or use `mad`."
        throw(ArgumentError(mess))
    end
end

function _model_from_univariate(
    est::AbstractEstimator,
    x::AbstractVector;
    method=:cg, # use :cg instead of :chol to avoid PosDefException
    dims=nothing, # explicit, so it is not passed to `rlm`
    kwargs...,
)

    check_l1loss(est)
    X = ones(eltype(x), (size(x, 1), 1))
    m = rlm(X, x, est; method=method, kwargs...)
end

function _mean(est::AbstractEstimator, x::AbstractVector; kwargs...)
    m = _model_from_univariate(est, x; kwargs...)
    coef(m)[1]
end

function _std(est::AbstractEstimator, x::AbstractVector; corrected::Bool=true, kwargs...)
    m = _model_from_univariate(est, x; kwargs...)
    n = dof_residual(m)
    if !corrected
        n += 1
    end
    dispersion(m.resp, n, false)
end

function _var(est::AbstractEstimator, x::AbstractVector; corrected::Bool=true, kwargs...)
    m = _model_from_univariate(est, x; kwargs...)
    n = dof_residual(m)
    if !corrected
        n += 1
    end
    dispersion(m.resp, n, true)
end

function _sem(est::AbstractEstimator, x::AbstractVector; kwargs...)
    m = _model_from_univariate(est, x; kwargs...)
    stderror(m)[1]
end

function _mean_and_std(
    est::AbstractEstimator,
    x::AbstractVector;
    corrected::Bool=true,
    kwargs...,
)
    m = _model_from_univariate(est, x; kwargs...)
    n = dof_residual(m)
    if !corrected
        n += 1
    end
    coef(m)[1], dispersion(m.resp, n, false)
end

function _mean_and_var(
    est::AbstractEstimator,
    x::AbstractVector;
    corrected::Bool=true,
    kwargs...,
)
    m = _model_from_univariate(est, x; kwargs...)
    n = dof_residual(m)
    if !corrected
        n += 1
    end
    coef(m)[1], dispersion(m.resp, n, true)
end

function _mean_and_sem(est::AbstractEstimator, x::AbstractVector; kwargs...)
    m = _model_from_univariate(est, x; kwargs...)
    coef(m)[1], stderror(m)[1]
end

## For arrays, along specific dimension with dims keyword argument

#Dims = Union{Colon,T,Tuple{Vararg{T}}} where T<:Int

function compatdims(array_ndims, dims)
    if isa(dims, Int)
        if 0 < dims <= array_ndims
            return dims
        end
    elseif isa(dims, Tuple)
        newdims = Tuple(d for d in dims if 0 < d <= array_ndims)
        if length(newdims) > 0
            return newdims
        end
    end
    return nothing
end

for fun in (:mean, :std, :var, :sem)
    _fun = Symbol("_$(fun)")
    @eval begin
        # `where Dims` to allow Colon
        function $(fun)(
            est::AbstractEstimator,
            x::AbstractArray;
            dims::Dims=:,
            kwargs...,
        ) where {Dims}
            check_l1loss(est)
            if dims === (:)
                return $(_fun)(est, vec(x); kwargs...)
            else
                dims = compatdims(ndims(x), dims)
                if isnothing(dims)
                    return x .* $(fun)(one(eltype(x)))
                else
                    return mapslices(r -> $(_fun)(est, vec(r); kwargs...), x; dims=dims)
                end
            end
        end
    end
end

for fun in (:mean_and_std, :mean_and_var, :mean_and_sem)
    _fun = Symbol("_$(fun)")
    @eval begin
        # `where Dims` to allow Colon
        function $(fun)(
            est::AbstractEstimator,
            x::AbstractArray;
            dims::Dims=:,
            kwargs...,
        ) where {Dims}
            check_l1loss(est)
            if dims === (:)
                ret = $(_fun)(est, vec(x); kwargs...)
                return ret
            else
                dims = compatdims(ndims(x), dims)
                if isnothing(dims)
                    m = mean(one(eltype(x)))
                    s = std(one(eltype(x)))  # NaN or NaN32
                    return (x * m, x * s)
                else
                    ret = mapslices(r -> $(_fun)(est, vec(r); kwargs...), x; dims=dims)
                    return (first.(ret), last.(ret))
                end
            end
        end
    end
end

## For iterators
for fun in (:mean, :std, :var, :sem, :mean_and_std, :mean_and_var, :mean_and_sem)
    _fun = Symbol("_$(fun)")
    @eval $(fun)(est::AbstractEstimator, itr; kwargs...) =
        $(_fun)(est, collect(itr); kwargs...)
end


text:
module MDToolbox

using Base.Threads: AtomicTypes
using FileIO
using NetCDF
#using NCDatasets
using NLsolve
using FFTW
#using Bio3DView
using CUDA
using ProgressMeter
using MetaGraphs, Graphs, GraphRecipes, Plots
using ChainRulesCore
import ChainRulesCore: rrule

# standard library
using Printf
using Statistics
using LinearAlgebra
using SparseArrays
using Dates
using Distributed
using Base.Threads

# package code goes here
export KB_kcalpermol, KB_kjpermol
export AbstractTrajectory, TrjArray, to3, select_atom
export mdload, mdsave, readdcd, writedcd, readnetcdf, writenetcdf, readpsf, writepsf, readpdb, writepdb, readnamdbin, writenamdbin, readcrd
export centerofmass, decenter, decenter!, orient!
export rotate, rotate!, rotate_with_matrix, superimpose, superimpose_serial, compute_rmsd, meanstructure, compute_rmsf
export compute_distance, compute_distancemap, compute_contactmap, compute_angle, compute_dihedral
export compute_qscore, compute_drms, compute_pairlist, compute_pairlist_bruteforce, compute_skrew
export clusterkcenters, clustercutoff, compute_cov, rsvd, pca, tica
export propagate_mcmc, propagate_md
export ksdensity, ksdensity_serial, compute_pmf
export wham, wham_iteration
export mbar, mbar_weight
export msmplot, msmgenerate, msmcountmatrix, msmtransitionmatrix, msmforward, msmforward_missing, msmbackward, msmbackward_missing, msmbaumwelch, msmbaumwelch_missing, msmviterbi, msmimpliedtime
export sp_delta_pmf, sp_design_matrix, sp_design_matrix_atom, sp_lsquares, sp_admm, sp_descent, sp_standardize!, sp_standardize, sp_cumulate_pmf, sp_cumulate_pmf_atom
export idilation, ierosion, itip_estimate!, surfing, afmize, AfmizeConfig, translateafm, getafmposterior, getposterior_parallel
export Asd, readasd
export get_residues, gpu, gpu32, cpu, cpu64, logsumexp

# constants
const KB_kcalpermol = 0.0019872041 #Boltzmann constant taken from wikipedia
const KB_kjpermol = 0.0083144621 #Boltzmann constant taken from wikipedia

# codes
include("trjarray.jl")
include("fileio.jl")
include("structure.jl")
include("reduction.jl")
include("dynamics.jl")
include("ksdensity.jl")
include("wham.jl")
include("mbar.jl")
include("msm.jl")
include("sparse.jl")
include("afm.jl")
include("asd.jl")
include("docking.jl")
include("utils.jl")

end # module


text:
module CuYao
using LuxurySparse, StaticArrays, LinearAlgebra, Base.Cartesian
using StatsBase
using BitBasis
using Reexport
import TupleTools
using Random

using YaoArrayRegister
using GPUArrays, CuArrays, CUDAnative
@reexport using Yao

const Ints = NTuple{<:Any, Int}

include("CUDApatch.jl")
include("GPUReg.jl")
include("gpuapplys.jl")
#include("gcompile.jl")

function __init__()
    CuArrays.allowscalar(false)
end

end


text:
#TODO make sure to renew V and S at each epoch

function layerUpdateParams!(
    model::Model,
    cLayer::CL,
    cnt::Integer = -1;
    tMiniBatch::Integer = 1,
    kwargs...,
) where {CL<:ConvLayer}

    optimizer = model.optimizer
    α = model.α
    β1, β2, ϵAdam = model.β1, model.β2, model.ϵAdam

    if cLayer.updateCount >= cnt
        return nothing
    end #if cLayer.updateCount >= cnt

    if !all(
        i -> (i.updateCount == cLayer.nextLayers[1].updateCount),
        cLayer.nextLayers,
    )
        return nothing
    end

    #initialize the needed variables to hold the corrected values
    #it is being init here cause these are not needed elsewhere

    if optimizer == :adam || optimizer == :momentum
        VCorrected = Dict(:dw => similar(cLayer.dW), :db => similar(cLayer.dB))

        if tMiniBatch == 1
            cLayer.V[:dw] .= 0
            cLayer.V[:db] .= 0
        end

        cLayer.V[:dw] .= β1 .* cLayer.V[:dw] .+ (1 - β1) .* cLayer.dW
        cLayer.V[:db] .= β1 .* cLayer.V[:db] .+ (1 - β1) .* cLayer.dB

        ##correcting
        VCorrected[:dw] .= cLayer.V[:dw] ./ (1 - β1^tMiniBatch)
        VCorrected[:db] .= cLayer.V[:db] ./ (1 - β1^tMiniBatch)

        if optimizer == :adam
            SCorrected =
                Dict(:dw => similar(cLayer.dW), :db => similar(cLayer.dB))

            if tMiniBatch == 1
                cLayer.S[:dw] .= 0
                cLayer.S[:db] .= 0
            end
            cLayer.S[:dw] .= β2 .* cLayer.S[:dw] .+ (1 - β2) .* (cLayer.dW .^ 2)
            cLayer.S[:db] .= β2 .* cLayer.S[:db] .+ (1 - β2) .* (cLayer.dB .^ 2)

            ##correcting
            SCorrected[:dw] .= cLayer.S[:dw] ./ (1 - β2^tMiniBatch)
            SCorrected[:db] .= cLayer.S[:db] ./ (1 - β2^tMiniBatch)

            ##update parameters with adam
            cLayer.W .-=
                (α .* (VCorrected[:dw] ./ (sqrt.(SCorrected[:dw]) .+ ϵAdam)))
            cLayer.B .-=
                (α .* (VCorrected[:db] ./ (sqrt.(SCorrected[:db]) .+ ϵAdam)))

            VCorrected = nothing
            SCorrected = nothing
        else#if optimizer==:momentum

            cLayer.W .-= (α .* VCorrected[:dw])
            cLayer.B .-= (α .* VCorrected[:db])
            VCorrected = nothing
        end #if optimizer==:adam
    else
        cLayer.W .-= (α .* cLayer.dW)
        cLayer.B .-= (α .* cLayer.dB)
    end #if optimizer==:adam || optimizer==:momentum

    cLayer.updateCount += 1

    # Base.GC.gc()

    return nothing
end #layerUpdateParams!

function layerUpdateParams!(
    model::Model,
    cLayer::PL,
    cnt::Integer = -1;
    tMiniBatch::Integer = 1,
    kwargs...,
) where {PL<:PoolLayer}

    if cLayer.updateCount >= cnt
        return nothing
    end #if cLayer.updateCount >= cnt

    if !all(
        i -> (i.updateCount == cLayer.nextLayers[1].updateCount),
        cLayer.nextLayers,
    )
        return nothing
    end

    cLayer.updateCount += 1

    return nothing
end #updateParams!


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/threadmessage.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_thread_message_queue_alloc,
    av_thread_message_queue_free,
    av_thread_message_queue_send,
    av_thread_message_queue_recv,
    av_thread_message_queue_set_err_send,
    av_thread_message_queue_set_err_recv

function av_thread_message_queue_alloc(mq,nelem::Integer,elsize::Integer)
    ccall((:av_thread_message_queue_alloc,libavutil),Cint,(Ptr{Ptr{AVThreadMessageQueue}},UInt32,UInt32),mq,nelem,elsize)
end

function av_thread_message_queue_free(mq)
    ccall((:av_thread_message_queue_free,libavutil),Void,(Ptr{Ptr{AVThreadMessageQueue}},),mq)
end

function av_thread_message_queue_send(mq,msg,flags::Integer)
    ccall((:av_thread_message_queue_send,libavutil),Cint,(Ptr{AVThreadMessageQueue},Ptr{Void},UInt32),mq,msg,flags)
end

function av_thread_message_queue_recv(mq,msg,flags::Integer)
    ccall((:av_thread_message_queue_recv,libavutil),Cint,(Ptr{AVThreadMessageQueue},Ptr{Void},UInt32),mq,msg,flags)
end

function av_thread_message_queue_set_err_send(mq,err::Integer)
    ccall((:av_thread_message_queue_set_err_send,libavutil),Void,(Ptr{AVThreadMessageQueue},Cint),mq,err)
end

function av_thread_message_queue_set_err_recv(mq,err::Integer)
    ccall((:av_thread_message_queue_set_err_recv,libavutil),Void,(Ptr{AVThreadMessageQueue},Cint),mq,err)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
import Random
using Plots

println("Pool Testing Optimization\nMonte Carlo Simulation")

min_total = 1
max_total = 100
min_infected = 0
iterations = 100

function run(total, infected, n_groups)
	# Trivial cases
	# No one is infected or testing everybody
	if infected==0 || n_groups==total
		return n_groups
	# Only 1 group (there is only 1 grouping step)
	elseif n_groups==1
		return total+1
	end

	people = zeros(Bool, total)
	people[1:infected] .= true
	Random.shuffle!(people)

	n_per_group = total ÷ n_groups
	n_plus_groups = total % n_groups

	groups = zeros(UInt, n_groups)

	for i in 1:infected
		j = 1+(Random.rand(UInt) % n_groups)
		groups[j] += 1
	end

	infected_groups = count(g -> g>0, groups)

	tests = n_groups

	for (index, value) in enumerate(groups)
		if value>0
			if index<=n_plus_groups
				tests += n_per_group+1
			else
				tests += n_per_group
			end
		end
	end

	#=
	println("Total   \t", total)
	println("Infected\t", infected)
	println("People\t", people)
	println("People per group\t", n_per_group)
	println("# groups with 1 more\t", n_plus_groups)
	println("Infected per group\t", groups)
	println("Number of infected groups\t", infected_groups)
	println("Number of tests\t", tests)
	=#

	return tests
end

n = binomial(max_total+2-min_infected, 2) - binomial(min_total+1-min_infected, 2)
results = Array{UInt, 2}(undef, n, 3)

i = 1
for total in min_total:max_total
	for infected in min_infected:total
		println("total ", total, "  infected ", infected)

		optimum_n_groups = total
		previous_sum_tests = iterations*total

		for n_groups in 1:total
			sum_tests = 0

			for iteration in 1:iterations
				tests = run(total, infected, n_groups)
				sum_tests += tests
				#println("total ", total, "  infected ", infected, "  n_groups ", n_groups, "  iteration ", iteration, "  tests ", tests)
			end

			if sum_tests<previous_sum_tests
				optimum_n_groups = n_groups
			end
		end

		global i
		results[i, :] = [total, infected, optimum_n_groups]
		i += 1
	end
end

println("Calculated")

matrix = zeros(Float32, (max_total+1-min_total, max_total+1-min_infected))
for k = 1:size(results, 1)
	local i = results[k, 1] + 1 - min_total
	local j = results[k, 2] + 1 - min_infected
	matrix[i, j] = results[k, 3]/results[k, 1]
end

#println("results ", results)
#println("matrix ", matrix)

heatmap(min_total:max_total, min_infected:max_total, matrix',
		xlabel="Total", ylabel="Infected",
		title="Optimal number of groups / Total",
		aspect_ratio=:equal,
		show=true)
println("Plotted")
readline()

println("End")


text:
"""
    `AbstractSpace{T,N_ARR,N_CDT}`
        with distance-type `T`,
        storable in array of dimension `N_ARR`,
        and with point-coordinates of dimension `N_CDT`.
    """
abstract type AbstractSpace{T,N_ARR,N_CDT} <: AbstractParameter{T} end

"""
    coordinates(space::AbstractSpace)

Return an object in the shape of the space where each element is the coordinate of that element.
"""
coordinates(space::AbstractSpace) = error("undefined.")

"""
    distances(space)

Return the distances between every pair of points in `space`
"""
function differences(space::AbstractSpace{T}) where T
    edges = Iterators.product(coordinates(space), coordinates(space))
    return difference.(Ref(space), edges)
end
function differences(space::AbstractSpace{T,N_ARR,N_CDT},
                     reference_location::NTuple{N_CDT,T}) where {T,N_ARR,N_CDT}
    edges = ((coord, reference_location) for coord in coordinates(space))
    return difference.(Ref(space), edges)
end

# Extend Base methods to AbstractSpace types
import Base: step, zero, length, size, ndims

zero(::Type{NTuple{N,T}}) where {N,T} = NTuple{N,T}(zero(T) for _ in 1:N)
zero(space::AbstractSpace{T}) where {T} = zeros(T,size(space)...)

ndims(space::AbstractSpace) = length(size(space))

function fft_center_idx(space::AbstractSpace)
    CartesianIndex(floor.(Ref(Int), size(space) ./ 2) .+ 1)
end
export fft_center_idx

# Include file definitions of various spaces
include("space/metrics.jl")
include("space/abstract_lattices.jl")
include("space/compact_lattices.jl")
include("space/periodic_lattices.jl")
include("space/lattice_augmentations.jl")