text:
# ---
# title: 417. Pacific Atlantic Water Flow
# id: problem417
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/pacific-atlantic-water-flow/description/>
# hidden: true
# ---
#
# Given an `m x n` matrix of non-negative integers representing the height of
# each unit cell in a continent, the "Pacific ocean" touches the left and top
# edges of the matrix and the "Atlantic ocean" touches the right and bottom
# edges.
#
# Water can only flow in four directions (up, down, left, or right) from a cell
# to another one with height equal or lower.
#
# Find the list of grid coordinates where water can flow to both the Pacific and
# Atlantic ocean.
#
# **Note:**
#
#   1. The order of returned grid coordinates does not matter.
#   2. Both _m_ and _n_ are less than 150.
#
#
#
# **Example:**
#
#
#
#     Given the following 5x5 matrix:
#
#       Pacific ~   ~   ~   ~   ~
#            ~  1   2   2   3  (5) *
#            ~  3   2   3  (4) (4) *
#            ~  2   4  (5)  3   1  *
#            ~ (6) (7)  1   4   5  *
#            ~ (5)  1   1   2   4  *
#               *   *   *   *   * Atlantic
#
#     Return:
#
#     [[0, 4], [1, 3], [1, 4], [2, 2], [3, 0], [3, 1], [4, 0]] (positions with parentheses in above matrix).
#
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function menu(savelocation::String = "."; compile::Bool = false)
    println('='^80)
    println("Policy Evaluation with Multiple Instrumental Variables")
    println("Mogstad, Torgovitsky, and Walters (2021)")
    println('='^80)
    println("What do you want to run?")
    println("\t 1. Illustration of mutual consistency (Figure 2)")
    println("\t 2. Numerical illustration for ATT (Figure 4)")
    println("\t 3. Numerical illustration for LATE(+20) (Figure 5)")
    println("\t 4. PRTE misspecification example (Figure 6)")
    println("\t 5. Everything")
    print("Enter choice: ")
    choice = readline()
    choice = parse(Int64, choice)

    if choice == 1
        savedir, _ = setup(savelocation, stub = "illustrate-mc")
        run_illustrate_mc(savedir, compile)
    elseif choice == 2
        savedir, _ = setup(savelocation, stub = "simulation-att")
        run_simulation_att(savedir, compile)
    elseif choice == 3
        savedir, _ = setup(savelocation, stub = "simulation-prte")
        run_simulation_prte(savedir, compile)
    elseif choice == 4
        savedir, _ = setup(savelocation, stub = "prte-misspecification")
        run_prte_misspecification(savedir, compile)
    elseif choice == 5
        savedir, _ = setup(savelocation, stub = "everything")
        results_mc = run_illustrate_mc(savedir, compile)
        results_att = run_simulation_att(savedir, compile)
        results_prte = run_simulation_prte(savedir, compile)
        results_prte_misspecification =
            run_prte_misspecification(savedir, compile)
        return Dict(:results_mc => results_mc,
                    :results_att => results_att,
                    :results_prte => results_prte)
    end
end
export menu

function run_illustrate_mc(savedir::String, compile::Bool = false)
    results = illustrate_mc()
    fn_results = joinpath(savedir, "illustrate-mc.csv")
    CSV.write(fn_results, results)
    if compile
        compile_latex(joinpath(savedir, "FigureMTRConsistent.tex"))
        compile_latex(joinpath(savedir, "FigureMTRInconsistent.tex"))
    end
    return results
end
export run_illustrate_mc

function run_simulation_att(savedir::String, compile::Bool = false)
    results = simulation_att()
    fn_results = joinpath(savedir, "simulation-att.csv")
    CSV.write(fn_results, results)
    if compile
        compile_latex(joinpath(savedir, "FigureATT.tex"))
    end
    return results
end
export run_simulation_att

function run_simulation_prte(savedir::String, compile::Bool = false)
    results = simulation_prte()
    fn_results = joinpath(savedir, "simulation-prte.csv")
    CSV.write(fn_results, results)
    if compile
        compile_latex(joinpath(savedir, "FigurePRTE.tex"))
    end
    return results
end
export run_simulation_prte

function run_prte_misspecification(savedir::String, compile::Bool = false)
    results = prte_misspecification()
    fn_results = joinpath(savedir, "prte-misspecification.csv")
    CSV.write(fn_results, results)
    if compile
        compile_latex(joinpath(savedir, "FigurePRTEMisspecification.tex"))
    end
    return results
end
export run_simulation_prte

function compile_latex(fn::String)
    oldwd = pwd()
    try
        cd(dirname(fn))
        cstr = `pdflatex -halt-on-error $(basename(fn)) "|" grep -a3 ^!`
        @suppress begin
            run(cstr)
            run(`latexmk -c`)
        end
        cd(oldwd)
    catch err
        cd(oldwd)
    end
end

ext: jl


text:
using SimIndex
using Base.Test
using Distances

println("Building an index on 1,000 1-dimensional points with Euclidean distance...")
si = SimIndex.Index(collect(1:1000), k=20, d=euclidean)
er = SimIndex.test_error_ratio(si)
println("Error ratio: $(er)")
println("")
@test er < 2.0

println("Building an index on 8,000 5-dimensional points with Euclidean distance...")
si = SimIndex.Index([string(i) => rand(5) for i=1:8000])
er = SimIndex.test_error_ratio(si)
println("Error ratio: $(er)")
println("")
@test er < 2.0

println("Building an index on 5,000 5-dimensional points with Cosine distance...")
coords = [string(i) => rand(5) for i=1:5000]
si = SimIndex.Index(coords, k=10, d=cosine_dist)
er = SimIndex.test_error_ratio(si)
println("Error ratio: $(er)")
println("")
@test er < 2.0

println("Testing incremental index construction...")
si = SimIndex.Index(collect(1:500), k=20, d=euclidean)
for i=501:1000
    push!(si, i)
end
compile!(si)
er = SimIndex.test_error_ratio(si)
println("Error ratio: $(er)")
println("")
@test er < 2.0

println("Testing index recompilation...")
si = SimIndex.Index(collect(1:1000), k=20, d=euclidean)
first_er = SimIndex.test_error_ratio(si)
println("Error ratio: $(er)")
compile!(si, 0.25)
second_er = SimIndex.test_error_ratio(si)
println("Error ratio: $(er)")
compile!(si, 0.05)
third_er = SimIndex.test_error_ratio(si)
println("Error ratio: $(er)")
@test first_er >= second_er
@test second_er >= third_er
@test third_er >= 1.0


text:
"""
`module  JAC.ImpactExcitationAutoion`
    ... a submodel of JAC that contains all methods for computing electron impact excitation cross sections and
        collision strengths.
"""
module ImpactExcitationAutoion

    using Printf, ..AngularMomentum, ..AutoIonization, ..Basics, ..ManyElectron, ..Radial, ..ImpactExcitation, ..TableStrings

    """
    `struct  ImpactExcitationAutoion.Settings`
        ... defines a type for the details and parameters of computing electron-impact excitation-autoionization pathways
            |i(N)>  --> |m(N)>  --> |f(N-1)>.

        + electronEnergies        ::Array{Float64,1}                   ... List of impact-energies of the incoming elecgtrons.
        + printBeforeComputation  ::Bool                               ... True, if all energies and lines are printed before their evaluation.
        + selectPathways          ::Bool                               ... True if particular pathways are selected for the computations.
        + selectedPathways        ::Array{Tuple{Int64,Int64,Int64},1}  ... List of list of pathways, given by tupels (inital, inmediate, final).
        + maxKappa                ::Int64                              ... Maximum kappa value of partial waves to be included.
    """
    struct Settings
        electronEnergies          ::Array{Float64,1}
        printBeforeComputation    ::Bool
        selectPathways            ::Bool
        selectedPathways          ::Array{Tuple{Int64,Int64,Int64},1}
        maxKappa                  ::Int64
    end

    """
    `ImpactExcitationAutoion.Settings()`  ... constructor for the default values of electron-impact excitation-autoionizaton settings.
    """
    function Settings()
       Settings( Float64[], false, false, Tuple{Int64,Int64,Int64}[], 0)
    end

    # `Base.show(io::IO, settings::ImpactExcitationAutoion.Settings)`
    #   ... prepares a proper printout of the variable settings::ImpactExcitationAutoion.Settings.
    function Base.show(io::IO, settings::ImpactExcitationAutoion.Settings)
        println(io, "electronEnergies:         $(settings.electronEnergies)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectPathways:           $(settings.selectPathways)  ")
        println(io, "selectedPathways:         $(settings.selectedPathways)  ")
        println(io, "maxKappa:                 $(settings.maxKappa)  ")
    end

    """
    `struct  ImpactExcitationAutoion.Channel`
        ... defines a type for a electron-impact excitaton & autoionization channel that specifies all quantum numbers, phases and
            amplitudes.

        + excitationChannel  ::ImpactExcitation.Channel      ... Channel that describes the electron-impact excitation process.
        + augerChannel       ::AutoIonization.Channel        ... Channel that describes the subsequent Auger/autoionization process.
    """
    struct  Channel
        excitationChannel    ::ImpactExcitation.Channel
        augerChannel         ::AutoIonization.Channel
    end

    """
    `struct  ImpactExcitationAutoion.Pathway`
        ... defines a type for a electron-impact excitation pathway that may include the definition of different excitation and
            autoionization channels and their corresponding amplitudes.

        + initialLevel        ::Level       ... initial-(state) level
        + intermediateLevel   ::Level       ... intermediate-(state) level
        + finalLevel          ::Level       ... final-(state) level
        + electronInEnergy    ::Float64     ... energy of the (incoming) electron
        + electronOutEnergy   ::Float64     ... energy of the (outgoing, scattered) electron
        + electronAugerEnergy ::Float64     ... energy of the (emitted Auger) electron
        + crossSection        ::Float64     ... total cross section of this pathway
        + hasChannels         ::Bool        ... Determines whether the individual excitation and autoionization
                                                channels are defined in terms of their free-electron kappa's, phases
                                                and the total angular momentum/parity as well as the amplitude, or not.
        + channels            ::Array{ImpactExcitationAutoion.Channel,1}  ... List of channels of this pathway.
    """
    struct  Pathway
        initialLevel          ::Level
        intermediateLevel     ::Level
        finalLevel            ::Level
        electronInEnergy      ::Float64
        electronOutEnergy     ::Float64
        electronAugerEnergy   ::Float64
        crossSection          ::Float64
        hasChannels           ::Bool
        channels              ::Array{ImpactExcitationAutoion.Channel,1}
    end

    """
    `ImpactExcitationAutoion.Pathway()`
        ... constructor for an electron-impact excitation-autoionization pathway between a specified initial, intermediate and
            final level.
    """
    function Line()
        Pathway(Level(), Level(), Level(), 0., 0., 0., 0., false, ImpactExcitationAutoion.Channel[] )
    end

    # `Base.show(io::IO, pathway::ImpactExcitationAutoion.Pathway)`
    #   ... prepares a proper printout of the variable pathway::ImpactExcitationAutoion.Pathway.
    function Base.show(io::IO, pathway::ImpactExcitationAutoion.Pathway)
        println(io, "initialLevel:               $(pathway.initialLevel)  ")
        println(io, "intermediateLevel:          $(pathway.intermediateLevel)  ")
        println(io, "finalLevel:                 $(pathway.finalLevel)  ")
        println(io, "electronInEnergy:           $(pathway.electronInEnergy)  ")
        println(io, "electronOutEnergy:          $(pathway.electronOutEnergy)  ")
        println(io, "electronAugerEnergy:        $(pathway.electronAugerEnergy)  ")
        println(io, "crossSection:               $(pathway.crossSection)  ")
        println(io, "hasChannels:                $(pathway.hasChannels)  ")
        println(io, "channels:                   $(pathway.channels)  ")
    end

    """
    `ImpactExcitationAutoion.computeAmplitudesProperties(pathway::ImpactExcitationAutoion.Pathway, grid::Radial.Grid,
                                                         settings::ImpactExcitationAutoion.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::ImpactExcitationAutoion.Line is returned for
            which the amplitudes and properties have now been evaluated.
    """
    function  computeAmplitudesProperties(pathway::ImpactExcitationAutoion.Pathway, grid::Radial.Grid, settings::ImpactExcitationAutoion.Settings)
        global JAC_counter
        newChannels = ImpactExcitationAutoion.Channel[]
        for channel in pathway.channels
            # Generate a continuum orbital
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("ImpactExcitationAutoion.computeAmplitudesProperties-aa: warning ... no cont. orb. is generated.")  end
            initialPhase  = 0.;    finalPhase = 0;    augerPhase = 0
            # Define a proper continuum basis from the finalLevel.basis and the continuum orbital
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("ImpactExcitationAutoion.computeAmplitudesProperties-ab: warning ... no cont. basis generated.") end
            # Compute the transition matrix for the continuum and the initial-state basis
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("ImpactExcitationAutoion.computeAmplitudesProperties-ac: warning ... no trans-matrix computed.") end
            # matrix    = PhotoExcitationAutoion.computeMatrix(channel.multipole, channel.gauge, line.omega, line.finalLevel.basis,
            #                                                      line.initialLevel.basis, grid, settings)
            # amplitude = line.finalLevel.mc * matrix * line.initialLevel.mc
            amplitude = 1.0
            exChannel = channel.excitationChannel
            eChannel  = ImpactExcitation.Channel( exChannel.initialKappa, exChannel.finalKappa, exChannel.symmetry,
                                                  initialPhase, finalPhase, Complex(0.))
            aChannel  = AutoIonization.Channel( channel.augerChannel.kappa, channel.augerChannel.symmetry, augerPhase, Complex(0.))
            push!( newChannels, ImpactExcitationAutoion.Channel(eChannel, aChannel) )
        end
        # Calculate the totalRate
        JAC_counter = JAC_counter + 1
        if   JAC_counter < 20   println("ImpactExcitationAutoion.computeAmplitudesProperties-ba: warning ... crossSection set to -1.") end
        crossSection = 1.
        pathway = ImpactExcitationAutoion.Pathway( pathway.initialLevel, pathway.intermediateLevel, pathway.finalLevel,
                                                   pathway.electronInEnergy, pathway.electronOutEnergy, pathway.electronAugerEnergy,
                                                   crossSection, true, newChannels)
        return( pathway )
    end

    """
    `ImpactExcitationAutoion.computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                             grid::Radial.Grid, settings::ImpactExcitationAutoion.Settings; output=true)`
        ... to compute the electron-impact-excitation-autoionization amplitudes and all properties as requested by the given settings.
            A list of lines::Array{ImpactExcitationAutoion.Lines} is returned.
    """
    function  computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                              settings::ImpactExcitationAutoion.Settings; output=true)
        println("")
        printstyled("ImpactExcitationAutoion.computePathways(): The computation of e-impact excitation-autoionization cross sections starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------------------------------------------ \n", color=:light_green)
        println("")
        #
        pathways = ImpactExcitationAutoion.determinePathways(finalMultiplet, intermediateMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBeforeComputation    ImpactExcitationAutoion.displayPathways(pathways)    end
        # Calculate all amplitudes and requested properties
        newPathways = ImpactExcitationAutoion.Pathway[]
        for  pathway in pathways
            newPathway = ImpactExcitationAutoion.computeAmplitudesProperties(pathway, grid, settings)
            push!( newPathways, newPathway)
        end
        # Print all results to screen
        ImpactExcitationAutoion.displayResults(pathways)
        #
        if    output    return( pathways )
        else            return( nothing )
        end
    end

    """
    `ImpactExcitationAutoion.determinePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                               settings::ImpactExcitationAutoion.Settings)`
        ... to determine a list of electron-impact excitation-autoionization pathways between the levels from the given initial-,
            intermediate- and final-state multiplets and by taking into account the particular selections and settings for this
            computation; an Array{ImpactExcitationAutoion.Line,1} is returned. Apart from the level specification, all physical
            properties are set to zero during the initialization process.
    """
    function  determinePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                settings::ImpactExcitationAutoion.Settings)
        if    settings.selectPathways    selectPathways = true;   selectedPathways = Basics.determine("selected pathways", settings.selectedPathways)
        else                             selectPathways = false
        end

        pathways = ImpactExcitationAutoion.Pathway[]
        for  i = 1:length(initialMultiplet.levels)
            for  n = 1:length(intermediateMultiplet.levels)
                for  f = 1:length(finalMultiplet.levels)
                    if  selectPathways  &&  !(haskey(selectedPathways, (i,n,f)) )    continue   end
                    for  en in settings.electronEnergies
                        eInEnergy  = en
                        eOutEnergy = en - (intermediateMultiplet.levels[n].energy - initialMultiplet.levels[i].energy)
                        eAEnergy   = intermediateMultiplet.levels[n].energy - finalMultiplet.levels[f].energy
                        if  eOutEnergy < 0.   ||   eAEnergy < 0    continue    end

                        channels = ImpactExcitationAutoion.determineChannels(finalMultiplet.levels[f], intermediateMultiplet.levels[n],
                                                                                initialMultiplet.levels[i], settings)
                        push!( pathways, ImpactExcitationAutoion.Pathway(initialMultiplet.levels[i], intermediateMultiplet.levels[n],
                                                                         finalMultiplet.levels[f], eInEnergy, eOutEnergy, eAEnergy,
                                                                         0., true, channels) )
                    end
                end
            end
        end
        return( pathways )
    end

    """
    `ImpactExcitationAutoion.determineChannels(finalLevel::Level, intermediateLevel::Level, initialLevel::Level,
                                                   settings::ImpactExcitationAutoion.Settings)`
        ... to determine a list of electron-impact excitation-autoionization Channels for a pathway from the initial to an
            intermediate and to a final level, and by taking into account the particular settings of for this computation;
            an Array{ImpactExcitationAutoion.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, intermediateLevel::Level, initialLevel::Level, settings::ImpactExcitationAutoion.Settings)
        symi      = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        symn      = LevelSymmetry(intermediateLevel.J, intermediateLevel.parity)
        # Determine first the electron-impact excitation channels
        eChannels = ImpactExcitation.Channel[];
        for  inKappa = -(settings.maxKappa+1):settings.maxKappa
            if  inKappa == 0    continue    end
            symtList = AngularMomentum.allowedTotalSymmetries(symi, inKappa)
            for  symt in symtList
                outKappaList = AngularMomentum.allowedKappaSymmetries(symt, symn)
                for  outKappa in outKappaList
                    push!(eChannels, ImpactExcitation.Channel( inKappa, outKappa, symt, 0., 0.,Complex(0.)) )
                end
            end
        end

        # Determine next the AutoIonization channels
        aChannels = AutoIonization.Channel[];
        kappaList = AngularMomentum.allowedKappaSymmetries(symn, symf)
        for  kappa in kappaList
            push!(aChannels, AutoIonization.Channel(kappa, symn, 0., Complex(0.)) )
        end

        # Now combine all these channels
        channels  = ImpactExcitationAutoion.Channel[];
        for    ec in eChannels
            for    a in aChannels
                push!(channels,  ImpactExcitationAutoion.Channel(ec, a) )
            end
        end

        return( channels )
    end

    """
    `ImpactExcitationAutoion.displayPathways(pathways::Array{ImpactExcitationAutoion.Line,1})`
        ... to display a list of pathways and channels that have been selected due to the prior settings. A neat table of all
            selected transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayPathways(pathways::Array{ImpactExcitationAutoion.Pathway,1})
        println(" ")
        println("  Selected electron-impact excitation-autoionization pathways:")
        println(" ")
        println("  ", TableStrings.hLine(179))
        sa = "     ";   sb = "     "
        sa = sa * TableStrings.center(23, "Levels"; na=2);            sb = sb * TableStrings.center(23, "i  --  m  --  f"; na=2);
        sa = sa * TableStrings.center(23, "J^P symmetries"; na=3);    sb = sb * TableStrings.center(23, "i  --  m  --  f"; na=3);
        sa = sa * TableStrings.center(58, "Energies  " * TableStrings.inUnits("energy");           na=4);
        sb = sb * TableStrings.center(58, "  m--i        m--f     e^-(in)     e^-(out)    e^-(Auger)"; na=4)
        sa = sa * TableStrings.flushleft(57, "List of partial waves and total symmetries"; na=4)
        sb = sb * TableStrings.flushleft(57, "(partial_in -> partial_out) J_total^P [partial_Auger] "; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(179))
        #
        for  pathway in pathways
            sa  = "  ";    isym = LevelSymmetry( pathway.initialLevel.J,      pathway.initialLevel.parity)
                           msym = LevelSymmetry( pathway.intermediateLevel.J, pathway.intermediateLevel.parity)
                           fsym = LevelSymmetry( pathway.finalLevel.J,        pathway.finalLevel.parity)
            sa = sa * TableStrings.center(23, TableStrings.levels_imf(pathway.initialLevel.index, pathway.intermediateLevel.index,
                                                                              pathway.finalLevel.index); na=2)
            sa = sa * TableStrings.center(23, TableStrings.symmetries_imf(isym, msym, fsym);  na=4)
            en_mi = pathway.intermediateLevel.energy - pathway.initialLevel.energy
            en_mf = pathway.intermediateLevel.energy - pathway.finalLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", en_mi))   * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", en_mf))   * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronInEnergy))    * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronOutEnergy))   * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronAugerEnergy)) * "    "
            kappaSymmetryList = Tuple{Int64,Int64,Int64,LevelSymmetry}[]
            for  i in 1:length(pathway.channels)
                eChannel = pathway.channels[i].excitationChannel;    aChannel = pathway.channels[i].augerChannel;
                push!( kappaSymmetryList, (eChannel.initialKappa, eChannel.finalKappa, aChannel.kappa, eChannel.symmetry) )
            end
            wa = TableStrings.kappaKappaKappaSymmetryTupels(57, kappaSymmetryList)
            if  length(wa) > 0    sb = sa * wa[1];    println( sb )    end
            for  i = 2:length(wa)
                sb = TableStrings.hBlank( length(sa) )
                if       i == 7   sb = sb * "... " * string( length(kappaSymmetryList) ) * " in total"
                elseif   i >  7   continue
                else     sb = sb * wa[i]
                end
                println( sb )
            end
        end
        println("  ", TableStrings.hLine(179))
        #
        return( nothing )
    end

    """
    `ImpactExcitationAutoion.displayResults(pathways::Array{ImpactExcitationAutoion.Line,1})`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but nothing is returned
            otherwise.
    """
    function  displayResults(pathways::Array{ImpactExcitationAutoion.Pathway,1})
        println(" ")
        println("  Electron-impact excitation-autoionization cross sections:")
        println(" ")
        println("  ", TableStrings.hLine(133))
        sa = "     ";   sb = "     "
        sa = sa * TableStrings.center(23, "Levels"; na=2);            sb = sb * TableStrings.center(23, "i  --  m  --  f"; na=2);
        sa = sa * TableStrings.center(23, "J^P symmetries"; na=3);    sb = sb * TableStrings.center(23, "i  --  m  --  f"; na=3);
        sa = sa * TableStrings.center(58, "Energies  " * TableStrings.inUnits("energy");           na=3);
        sb = sb * TableStrings.center(58, "  m--i        m--f     e^-(in)     e^-(out)    e^-(Auger)"; na=3)
        sa = sa * TableStrings.center(16, "Cross sections"; na=2);
        sb = sb * TableStrings.center(16, TableStrings.inUnits("cross section"); na=2)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(133))
        #
        for  pathway in pathways
            sa  = "  ";    isym = LevelSymmetry( pathway.initialLevel.J,      pathway.initialLevel.parity)
                           msym = LevelSymmetry( pathway.intermediateLevel.J, pathway.intermediateLevel.parity)
                           fsym = LevelSymmetry( pathway.finalLevel.J,        pathway.finalLevel.parity)
            sa = sa * TableStrings.center(23, TableStrings.levels_imf(pathway.initialLevel.index, pathway.intermediateLevel.index,
                                                                              pathway.finalLevel.index); na=2)
            sa = sa * TableStrings.center(23, TableStrings.symmetries_imf(isym, msym, fsym);  na=4)
            en_mi = pathway.intermediateLevel.energy - pathway.initialLevel.energy
            en_mf = pathway.intermediateLevel.energy - pathway.finalLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", en_mi))   * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", en_mf))   * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronInEnergy))      * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronOutEnergy))     * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronAugerEnergy))   * "     "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", pathway.crossSection))   * "  "
            println(sa)
        end
        println("  ", TableStrings.hLine(133))
        #
        return( nothing )
    end

end # module


text:
function gausspdf(x, mu, sigma)
	k = 1/(sigma * sqrt(2*pi))
	p = k * exp(-((x-mu)^2/(2*sigma^2)))
end

function initProbs(w, mu, sigma)
	f(x) = gausspdf(x, mu, sigma)
	init = map(f, [1:w])
	init / sum(init)
end

function possiblePrevPositions(prev_len, delta, i)
	if i >= 2 && i <= prev_len + delta
		[ max(1, i-delta) : 1 : min(i-1, prev_len) ]
	else
		nothing
	end
end

function forwardProbs(prevProbs, w, mu, sigma)
	delta = w
	n_prev = length( prevProbs )
	nextProbs = zeros(n_prev + w - 1)
	n_next = length(nextProbs)
	for i = 2 : (n_next+1)
		prevPositions = possiblePrevPositions( n_prev, delta, i)
		if prevPositions != nothing
			prevProbs1 = prevProbs[prevPositions]
			deltas = map(x -> i - x, prevPositions)
			prob_deltas = map(x -> gausspdf(x, mu, sigma), deltas)
			nextProbs[i-1] = sum(prob_deltas .* prevProbs1)
		end
	end
	nextProbs = nextProbs / sum(nextProbs)
end

function pccf(num_events, w, mu, sigma)
	tic()
	allM = zeros(num_events * w, num_events)
	allM[1:w, 1] = initProbs(w, mu, sigma)
	for i = 2:num_events
		prev_vect = allM[(i-1):((i-1)*w), i-1]
		next_vec_of_probs = forwardProbs(prev_vect, w, mu, sigma)
		allM[(i:(i*w)), i] = next_vec_of_probs
	end
	nn, mm = size(allM)
	cc = zeros(nn)
	for i = 1:nn
		s = 0.
		for j = 1:mm
			s += allM[i,j]
		end
		cc[i] = s
	end
	toc()
	cc = cc / sum(cc)
end

w = 150
m = 25.0
s = 3.2
r = pccf(30, w, m, s)

# Write results to output file
outfile = open("C:/.../pccfjl.txt", "w")
write(outfile, join(r, "\n"))
close(outfile)

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
using DSGE, HDF5, JLD
using BenchmarkTools
using Base.Test

path = dirname(@__FILE__)

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:saveroot, tempdir())
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m <= Setting(:forecast_uncertainty_override, Nullable(false))
m <= Setting(:use_population_forecast, true)
m <= Setting(:compute_shockdec_bands, true)

estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "metropolis_hastings.h5")

output_vars = add_requisite_output_vars([:histpseudo, :histobs,
                                         :hist4qpseudo, :hist4qobs,
                                         :forecastpseudo, :forecastobs,
                                         :forecast4qpseudo, :forecast4qobs,
                                         :shockdecpseudo, :shockdecobs,
                                         :irfpseudo, :irfobs])

##########
# Modal
##########
fc_trial = @benchmark forecast_one(m, :mode, :none, output_vars, verbose = :none) gcsample = true
mb_trial = @benchmark compute_meansbands(m, :mode, :none, output_vars;
                                         compute_shockdec_bands = true, verbose = :none) gcsample = true
mb2m_trial = @benchmark meansbands_to_matrix(m, :mode, :none, output_vars; verbose =
                                             :none) gcsample = true

trials = [fc_trial, mb_trial, mb2m_trial]
trial_names = [:forecast_one, :compute_meansbands, :meansbands_to_matrix]
group = construct_trial_group(trials, trial_names)
group_name = "modal_meansbands"

# write_ref_trial_group(group, trial_names, group_name)
print_all_benchmarks(group, "../reference/$group_name.jld", group_name, trial_names)

#######################
# Full-distribution
#######################
@everywhere using DSGE
m <= Setting(:forecast_block_size, 5)
fc_trial = @benchmark forecast_one(m, :full, :none, output_vars, verbose = :none) gcsample = true
mb_trial = @benchmark compute_meansbands(m, :full, :none, output_vars;
                                         compute_shockdec_bands = true,
                                         verbose = :none) gcsample = true
mb2m_trial = @benchmark meansbands_to_matrix(m, :full, :none, output_vars;
                                             verbose = :none) gcsample = true

trials = [fc_trial, mb_trial, mb2m_trial]
trial_names = [:forecast_one, :compute_meansbands, :meansbands_to_matrix]
group = construct_trial_group(trials, trial_names)
group_name = "full_meansbands"

# write_ref_trial_group(group, trial_names, group_name)
print_all_benchmarks(group, "../reference/$group_name.jld", group_name, trial_names)


text:
using AxisKeys
using DataFrames, CSV
using StanSample
using Statistics, Test

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"), DataFrame);

# Define the Stan language model

stan10_4 = "
data{
    int N;
    int N_actors;
    int pulled_left[N];
    int prosoc_left[N];
    int condition[N];
    int actor[N];
}
parameters{
    vector[N_actors] a;
    real bp;
    real bpC;
}
model{
    vector[N] p;
    bpC ~ normal( 0 , 10 );
    bp ~ normal( 0 , 10 );
    a ~ normal( 0 , 10 );
    for ( i in 1:504 ) {
        p[i] = a[actor[i]] + (bp + bpC * condition[i]) * prosoc_left[i];
        p[i] = inv_logit(p[i]);
    }
    pulled_left ~ binomial( 1 , p );
}
";

data = (N = size(df, 1), N_actors = length(unique(df.actor)),
    actor = df.actor, pulled_left = df.pulled_left,
    prosoc_left = df.prosoc_left, condition = df.condition);

# Sample using cmdstan

m10_4s = SampleModel("m10.4s", stan10_4)
rc10_4s = stan_sample(m10_4s; data);

# Result rethinking

rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
bp    0.84 0.26  0.43  1.26  2271    1
bpC  -0.13 0.29 -0.59  0.34  2949    1

a[1] -0.74 0.27 -1.16 -0.31  3310    1
a[2] 10.88 5.20  4.57 20.73  1634    1
a[3] -1.05 0.28 -1.52 -0.59  4206    1
a[4] -1.05 0.28 -1.50 -0.60  4133    1
a[5] -0.75 0.27 -1.18 -0.32  4049    1
a[6]  0.22 0.27 -0.22  0.65  3877    1
a[7]  1.81 0.39  1.22  2.48  3807    1
";

# Update sections

@testset "KeyedArray manipulations" begin
    if success(rc10_4s)
        ka = read_samples(m10_4s, :keyedarray)
        kb = ka(Symbol("a.1"))

        # Other manipulations

        @test Tables.istable(ka) == true
        # All draws of :a_1
        @test size(vcat(kb...)) == (4000,)
        # All of parameters @testset "Basic HelpModel" begin
        kar = reshape(ka.data, 4000, 9);
        @test size(kar) == (4000, 9)
        # Axes ranges
        @test axes(ka) == (Base.OneTo(1000), Base.OneTo(4), Base.OneTo(9))
        # Axes keys
        @test axiskeys(ka) == (
            1:1000, 1:4,
            [Symbol("a.1"), Symbol("a.2"), Symbol("a.3"), Symbol("a.4"),
                Symbol("a.5"), Symbol("a.6"), Symbol("a.7"), :bp, :bpC]
        )
        # A single axis
        @test axiskeys(ka, :param) == vcat([Symbol("a.$i") for i in 1:7], :bp, :bpC)

        # Test combining vector param 'a'
        ma = matrix(ka, "a");
        rma = reshape(ma.data, 4000, size(ma, 3))
        @test mean(rma, dims=1) ≈ [-0.7 10.9 -1 -1 -0.7 0.2 1.8] atol=0.7
    end
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
module WfnUtils

using LinearAlgebra
using ..MapViews

export applygate!

"""
Generates a mapping that moves qBit k in N qBits to the top
To be used internally by multibubble
O(2^N) complexity
"""
function bubble(k::Int, N::Int)
    index_list = collect(1:2^N)
    part_size = 2^(N - abs(k))
    partitions = collect(Iterators.partition(index_list, part_size))

    odds  = 1:2:length(partitions)
    evens = 2:2:length(partitions)

    if k > 0
        part_map = collect(Iterators.flatten((odds, evens)))
        # if k < 0 generate a map correspond to the qBit in question being NOT'ed
    else
        part_map = collect(Iterators.flatten((evens, odds)))
    end

    mapped_parts = MapView(part_map, partitions)
    collect(Iterators.flatten(mapped_parts))
end

"""
Generates a mapping such that the qbits in ks are at the top of N qBits
With the top to bottom ordering corresponds to begin to end of ks
Negative k indicate control-on-zero
Example:
    Let there be 5 qBit lines in order:
        ψ₁ ψ₂ ψ₃ ψ₄ ψ₅
  multibubble([2, 4], 5) would generate a mapping that corresponds to:
        ψ₂ ψ₄ ψ₁ ψ₃ ψ₅

Complexity: O(n * 2^N + n^2)
where:
        N is the number of qBits
        n is the number of qBits to be bubbled
"""
function multibubble(ks::Vector{Int}, N::Int)
    final_map = collect(1:2^N)

    push_down_list = []
    for k in reverse(ks)
        # separate sign from value
        (sign, k) = (div(k, abs(k)), abs(k))

        push_down = 0
        for i in push_down_list
            if k < i
                push_down += 1
            elseif k == i
                throw("bubbling list has dublicate entry")
            end
        end
        push!(push_down_list, k)

        map = bubble(sign * (k + push_down), N)
        final_map = collect(MapView(map, final_map))
    end

    final_map
end

"""
TODO: document
"""
function applygate!(
        ψ::V,
        Gate::M,
        target::Vector{Int};
        control::Vector{Int} = Vector{Int}()
        ) where {V <: AbstractVector{<: Number}, M <: AbstractMatrix{<: Number}}

    N = Int(log2(length(ψ)))
    n = length(target)

    @assert size(Gate)[1] == size(Gate)[2] "Invalid gate"
    @assert size(Gate)[1] == 2^n "Mismatch gate and # of input"

    top_qbit_list = [control; target]
    partition_size = div(length(ψ), 2^length(top_qbit_list))

    @assert max(top_qbit_list...) <= N "qBit number larger than number of qBits"

    # Generate a mapping that bubbles all relevent qBits to the top
    # Since in controlled gates' operations,
    #   only the bottom right of the matrix (the Gate application) is non-trivial,
    #   only takes the end side of ψ that corresponds to the Gate
    bubble_map = multibubble(top_qbit_list, N)
    gate_index = length(ψ) - size(Gate)[1] * partition_size + 1
    map_ψ = MapView(bubble_map[gate_index:end], ψ)

    # Apply Gate on ψ
    map_ψ_part = collect(Iterators.partition(map_ψ, partition_size))
    map_ψ .= collect(Iterators.flatten(Gate * map_ψ_part)) # ψ mutated here
end

end


text:
module CMSSPs

# For re-exporting
using Reexport

# Stdlib Requirements
using Logging
using Printf
using Random
using Statistics

# Common non-stdlib requirements
using DataStructures
using LinearAlgebra
using StaticArrays
using Distributions
using PDMats

# Parsing and support
using TOML

# Include submodule files
include("models/CMSSPModel.jl")
include("hhpc/HHPC.jl")
include("domains/CMSSPDomains.jl")

@reexport using CMSSPs.CMSSPModel
@reexport using CMSSPs.HHPC
@reexport using CMSSPs.CMSSPDomains

end # module


text:
module NetCDF

using NetCDF_jll
using Formatting
using Base.Cartesian
import DiskArrays: readblock!, writeblock!, AbstractDiskArray, eachchunk, GridChunks,
       estimate_chunksize, haschunks, Chunked, Unchunked

include("netcdf_constants.jl")
include("netcdf_c.jl")

import Base.show

export NcDim,
       NcVar,
       NcFile,
       ncread,
       ncread!,
       ncwrite,
       nccreate,
       ncinfo,
       ncclose,
       ncputatt,
       NC_BYTE,
       NC_SHORT,
       NC_INT,
       NC_FLOAT,
       NC_DOUBLE,
       NC_STRING,
       ncgetatt,
       NC_NOWRITE,
       NC_WRITE,
       NC_CHAR,
       NC_CLOBBER,
       NC_NOCLOBBER,
       NC_CLASSIC_MODEL,
       NC_64BIT_OFFSET,
       NC_NETCDF4,
       NC_UNLIMITED,
       nc_char2string,
       nc_string2char

@deprecate ncclose() nothing
@deprecate ncsync() nothing
@deprecate ncclose(::AbstractString) nothing

NC_VERBOSE = false
#Some constants
#Define a short ASCIIChar type to map to NC_CHAR,
#We can not use Julias Char because it has 4 bytes size
#create a new AbstractChar type to test the fallbacks,
#this is re-used from the julialang tests
#https://github.com/JuliaLang/julia/blob/master/test/char.jl#L242
primitive type ASCIIChar <: AbstractChar 8 end
ASCIIChar(c::UInt8) = reinterpret(ASCIIChar, c)
ASCIIChar(c::UInt32) = ASCIIChar(UInt8(c))
Base.codepoint(c::ASCIIChar) = reinterpret(UInt8, c)
Base.zero(::Type{ASCIIChar}) = ASCIIChar(0)

global const nctype2jltype = Dict(
    NC_BYTE => Int8,
    NC_UBYTE => UInt8,
    NC_SHORT => Int16,
    NC_USHORT => UInt16,
    NC_INT => Int32,
    NC_UINT => UInt32,
    NC_INT64 => Int64,
    NC_UINT64 => UInt64,
    NC_FLOAT => Float32,
    NC_DOUBLE => Float64,
    NC_CHAR => ASCIIChar,
    NC_STRING => String,
)

global const nctype2string = Dict(
    NC_BYTE => "BYTE",
    NC_UBYTE => "UBYTE",
    NC_SHORT => "SHORT",
    NC_USHORT => "USHORT",
    NC_INT => "INT",
    NC_UINT => "UINT",
    NC_INT64 => "INT64",
    NC_UINT64 => "UINT64",
    NC_FLOAT => "FLOAT",
    NC_DOUBLE => "DOUBLE",
    NC_STRING => "STRING",
    NC_CHAR => "CHAR",
)

function jl2nc(t::DataType)
    popfirst!(collect(keys(nctype2jltype))[findall(
        e -> (t <: e),
        collect(values(nctype2jltype)),
    )])
end
jl2nc(t::Type{UInt8}) = NC_UBYTE

getJLType(t::DataType) = t
getJLType(t::Int) = nctype2jltype[t]

getNCType(t::DataType) = jl2nc(t)
getNCType(t::Int) = Int(t)

"""
    NcDim

Represents a NetCDF dimension of name `name` optionally holding the dimension values.
"""
mutable struct NcDim
    ncid::Int32
    dimid::Int32
    varid::Int32
    name::String
    dimlen::UInt
    vals::AbstractArray
    atts::Dict
    unlim::Bool
end

"""
    NcDim(name::String,dimlength::Integer)

Creates a NetCDF dimension object named `name` with length `dimlength` in memory.

### Keyword arguments

* `values=[]` optional for adding dimension values which can be written to a dimension variable
* `atts` a Dict of attributes associated to the dimension variable
* `unlimited` is this the record dimension, defaults to `false`
"""
function NcDim(
    name::AbstractString,
    dimlength::Integer;
    values::Union{AbstractArray,Number} = [],
    atts::Dict = Dict{Any,Any}(),
    unlimited = false,
)
    length(values) > 0 &&
    dimlength != 0 &&
    length(values) != dimlength &&
    error("Dimension value vector must have the same length as dimlength!")
    NcDim(-1, -1, -1, string(name), dimlength, collect(values), atts, unlimited)
end

"""
NcDim(name::String,values::AbstractVector)

Creates a NetCDF dimension object named `name` and associated `values` in memory.

### Keyword arguments

* `atts` a Dict of attributes associated to the dimension variable
* `unlimited` is this the record dimension, defaults to `false`
"""
NcDim(
    name::AbstractString,
    values::AbstractArray;
    atts::Dict = Dict{Any,Any}(),
    unlimited = false,
) = NcDim(name, length(values), values = values, atts = atts, unlimited = unlimited)
NcDim(name::AbstractString, values::AbstractArray, atts::Dict; unlimited = false) =
    NcDim(name, length(values), values = values, atts = atts, unlimited = unlimited)

"""
    NcVar

`NcVar{T,N,M}` represents a NetCDF variable. It is a subtype of `AbstractArray{T,N}`, so normal indexing using `[]`
will work for reading and writing data to and from a NetCDF file. `NcVar` objects are returned by `NetCDF.open`, by
indexing an `NcFile` object (e.g. `myfile["temperature"]`) or, when creating a new file, by its constructor. The type parameter `M`
denotes the NetCDF data type of the variable, which may or may not correspond to the Julia Data Type.
"""
mutable struct NcVar{T,N,M} <: AbstractDiskArray{T,N}
    ncid::Int32
    varid::Int32
    ndim::Int32
    natts::Int32
    nctype::Int32
    name::String
    dimids::Vector{Int32}
    dim::Vector{NcDim}
    atts::Dict
    compress::Int32
    chunksize::NTuple{N,Int32}
    ncfile #Kept to prevent closing through finalizer
end

Base.convert(::Type{NcVar{T,N,M}}, v::NcVar{S,N,M}) where {S,T,N,M} = NcVar{T,N,M}(
    v.ncid,
    v.varid,
    v.ndim,
    v.natts,
    v.nctype,
    v.name,
    v.dimids,
    v.dim,
    v.atts,
    v.compress,
    v.chunksize,
    v.ncfile,
)
Base.convert(::Type{NcVar{T}}, v::NcVar{S,N,M}) where {S,T,N,M} = NcVar{T,N,M}(
    v.ncid,
    v.varid,
    v.ndim,
    v.natts,
    v.nctype,
    v.name,
    v.dimids,
    v.dim,
    v.atts,
    v.compress,
    v.chunksize,
    v.ncfile
)

# Implement the DiskArrays interface
function readblock!(v::NcVar, aout, r::AbstractUnitRange...)
  readvar!(v,aout,start = [first(i) for i in r], count = [length(i) for i in r])
end
function writeblock!(v::NcVar, a, r::AbstractUnitRange...)
  putvar(v,a,start = [first(i) for i in r], count = [length(i) for i in r])
  sync(v.ncfile)
end
getchunksize(v::NcVar) = getchunksize(haschunks(v),v)
getchunksize(::Chunked, v::NcVar) = reverse(map(Int64,v.chunksize))
getchunksize(::Unchunked, v::NcVar) = estimate_chunksize(v)
eachchunk(v::NcVar) = GridChunks(v, getchunksize(v))
haschunks(v::NcVar) = all(iszero,v.chunksize) ? Unchunked() : Chunked()

"""
    NcVar(name::AbstractString,dimin::Union{NcDim,Array{NcDim,1}}

Creates a new NetCDF variable `name` in memory. The variable is associated to a
list of NetCDF dimensions specified by `dimin`.

### Keyword arguments

* `atts` Dictionary representing the variables attributes
* `t` either a Julia type, (one of `Int16`, `Int32`, `Float32`, `Float64`, `String`) or a NetCDF data type (`NC_SHORT`, `NC_INT`, `NC_FLOAT`, `NC_DOUBLE`, `NC_CHAR`, `NC_STRING`) determines the data type of the variable. Defaults to -1
* `compress` Integer which sets the compression level of the variable for NetCDF4 files. Defaults to -1 (no compression). Compression levels of 1..9 are valid
"""
function NcVar(
    name::AbstractString,
    dimin::Union{NcDim,Array{NcDim,1}};
    atts::Dict = Dict{Any,Any}(),
    t::Union{DataType,Integer} = Float64,
    compress::Integer = -1,
    chunksize::Tuple = ntuple(i -> zero(Int32), isa(dimin, NcDim) ? 1 : length(dimin)),
)
    dim = isa(dimin, NcDim) ? NcDim[dimin] : dimin
    chunksize = map(Int32, chunksize)
    return NcVar{getJLType(t),length(dim),getNCType(t)}(
        Int32(-1),
        Int32(-1),
        Int32(length(dim)),
        Int32(length(atts)),
        getNCType(t),
        name,
        fill(Int32(-1), length(dim)),
        dim,
        atts,
        Int32(compress),
        chunksize,
        nothing,
    )
end

#Array methods
@generated function Base.size(a::NcVar{T,N}) where {T,N}
    :(@ntuple($N, i -> Int(a.dim[i].dimlen)))
end
const IndR = Union{Integer,UnitRange,Colon}
const ArNum = Union{AbstractArray,Number}

"""
    NcFile

Represents a link to a NetCDF file. Is returned by `NetCDF.open` or by `NetCDF.create`.
"""
mutable struct NcFile
    ncid::Int32
    nvar::Int32
    ndim::Int32
    ngatts::Int32
    vars::Dict{String,NcVar}
    dim::Dict{String,NcDim}
    gatts::Dict
    unlimdimids::Vector{Int32}
    name::String
    omode::UInt16
    in_def_mode::Bool
end
#Define getindex method to retrieve a variable
Base.getindex(nc::NcFile, i::AbstractString) = haskey(nc.vars, i) ? nc.vars[i] :
error("NetCDF file $(nc.name) does not have a variable named $(i)")
Base.keys(nc::NcFile) = keys(nc.vars)

include("netcdf_helpers.jl")

readvar!(
    nc::NcFile,
    varname::AbstractString,
    retvalsa::AbstractArray;
    start::Vector = defaultstart(nc[varname]),
    count::Vector = defaultcount(nc[varname]),
) = readvar!(nc[varname], retvalsa, start = start, count = count)

"""
    NetCDF.readvar!(v::NcVar, d::Array;start::Vector=ones(UInt,ndims(d)),count::Vector=size(d))

Reads the values from the file associated to the `NcVar` object `v` into the array `d`. By default the whole variable is read

### Keyword arguments

* `start` Vector of length `ndim(v)` setting the starting index for each dimension
* `count` Vector of length `ndim(v)` setting the count of values to be read along each dimension. The value -1 is treated as a special case to read all values from this dimension

### Example

Assume `v` is a NetCDF variable with dimensions (3,3,10).

    x = zeros(3,1,10)
    NetCDF.readvar(v, x, start=[1,2,1], count=[3,1,-1])

This reads all values from the first and last dimension and only the second value from the second dimension.
"""
function readvar!(
    v::NcVar,
    retvalsa::AbstractArray;
    start::Vector = defaultstart(v),
    count::Vector = defaultcount(v),
)

    isa(retvalsa, Array) ||
    Base.iscontiguous(retvalsa) || error("Can only read into contiguous pieces of memory")

    length(start) == v.ndim || error("Length of start ($(length(start))) must equal the number of variable dimensions ($(v.ndim))")
    length(count) == v.ndim || error("Length of start ($(length(count))) must equal the number of variable dimensions ($(v.ndim))")

    p = preparestartcount(start, count, v)

    length(retvalsa) != p && error("Size of output array ($(length(retvalsa))) does not equal number of elements to be read ($p)!")

    nc_get_vara_x!(v, gstart, gcount, retvalsa)

    retvalsa
end

readvar(
    nc::NcFile,
    varname::AbstractString;
    start::Vector = defaultstart(nc[varname]),
    count::Vector = defaultcount(nc[varname]),
) = readvar(nc[varname], start = start, count = count)

"""
    NetCDF.readvar(v::NcVar;start::Vector=ones(UInt,ndims(d)),count::Vector=size(d))

Reads the values from the file associated to the `NcVar` object `v` and returns them. By default the whole variable is read

### Keyword arguments

* `start` Vector of length `ndim(v)` setting the starting index for each dimension
* `count` Vector of length `ndim(v)` setting the count of values to be read along each dimension. The value -1 is treated as a special case to read all values from this dimension

### Example

Assume `v` is a NetCDF variable with dimensions (3,3,10).

    x = NetCDF.readvar(v, start=[1,2,1], count=[3,1,-1])

This reads all values from the first and last dimension and only the second value from the second dimension.
"""
function readvar(
    v::NcVar{T,N};
    start::Vector = defaultstart(v),
    count::Vector = defaultcount(v),
) where {T,N}
    s = [count[i] == -1 ? size(v, i) - start[i] + 1 : count[i] for i = 1:length(count)]
    retvalsa = Array{T}(undef, s...)
    readvar!(v, retvalsa, start = start, count = count)
    return retvalsa
end

"""
    NetCDF.readvar{T,N}(v::NcVar{T,N},I::Union{Integer, UnitRange, Colon}...)

Reads data from a NetCDF file with array-style indexing. `Integer`s and `UnitRange`s and `Colon`s are valid indices for each dimension.
"""
function readvar(v::NcVar{T,N}, I::IndR...) where {T,N}
    count = ntuple(i -> counti(I[i], v.dim[i].dimlen), length(I))
    retvalsa = Array{T}(undef, count...)
    readvar!(v, retvalsa, I...)
    return retvalsa
end

# Here are some functions for array-style indexing readvar
#For single indices
@generated function readvar(v::NcVar{T,N}, I::Integer...) where {T,N}
    N == length(I) || error("Dimension mismatch")
    if N == 0
        quote
            checkbounds(v, I...)
            gstart[1] = 0
            nc_get_var1_x(v, gstart, T)::T
        end
    else
        quote
            checkbounds(v, I...)
            @nexprs $N i -> gstart[v.ndim+1-i] = I[i] - 1
            nc_get_var1_x(v, gstart, T)::T
        end
    end
end

firsti(i::Integer, l::Integer) = i - 1
counti(i::Integer, l::Integer) = 1
firsti(r::UnitRange, l::Integer) = first(r) - 1
counti(r::UnitRange, l::Integer) = length(r)
firsti(r::Colon, l::Integer) = 0
counti(r::Colon, l::Integer) = Int(l)

# For ranges
"""
    NetCDF.readvar!{T,N}(v::NcVar{T,N},d::AbstractArray, I::Union{Integer, UnitRange, Colon}...)

Reads data from a NetCDF file with array-style indexing and writes them to d. `Integer`s and `UnitRange`s and `Colon`s are valid indices for each dimension.
"""
@generated function readvar!(v::NcVar{T,N}, retvalsa::AbstractArray, I::IndR...) where {T,N}

    N == length(I) || error("Dimension mismatch")

    quote
        checkbounds(v, I...)

        isa(retvalsa, Array) ||
        Base.iscontiguous(retvalsa) ||
        error("Can only read into contiguous pieces of memory")

        @nexprs $N i -> gstart[v.ndim+1-i] = firsti(I[i], v.dim[i].dimlen)
        @nexprs $N i -> gcount[v.ndim+1-i] = counti(I[i], v.dim[i].dimlen)
        p = 1
        @nexprs $N i -> p = p * gcount[v.ndim+1-i]
        length(retvalsa) != p && error(string(
            "Size of output array ($(length(retvalsa))) does not equal number of elements to be read (",
            p,
            ")!",
        ))
        nc_get_vara_x!(v, gstart, gcount, retvalsa)
        retvalsa
    end
end

for (t, ending, arname) in funext
    fname = Symbol("nc_get_vara_$ending")
    fname1 = Symbol("nc_get_var1_$ending")
    arsym = Symbol(arname)
    @eval nc_get_vara_x!(
        v::NcVar{$t},
        start::Vector{UInt},
        count::Vector{UInt},
        retvalsa::AbstractArray{$t},
    ) = $fname(v.ncid, v.varid, start, count, retvalsa)
    @eval nc_get_var1_x(v::NcVar{$t}, start::Vector{UInt}, ::Type{$t}) = begin
        $fname1(v.ncid, v.varid, start, $(arsym))
        $(arsym)[1]
    end
end

nc_get_vara_x!(
    v::NcVar{ASCIIChar,N,NC_CHAR},
    start::Vector{UInt},
    count::Vector{UInt},
    retvalsa::AbstractArray{<:ASCIIChar},
) where {N} = nc_get_vara_text(v.ncid, v.varid, start, count, retvalsa)
nc_get_var1_x!(
    v::NcVar{ASCIIChar,N,NC_CHAR},
    start::Vector{UInt},
    retvalsa::AbstractArray{<:ASCIIChar},
) where {N} = nc_get_var1_text(v.ncid, v.varid, start, retvalsa)

function nc_get_vara_x!(
    v::NcVar{String,N,NC_STRING},
    start::Vector{UInt},
    count::Vector{UInt},
    retvalsa::AbstractArray{<:String},
) where {N}
    @assert length(retvalsa) == prod(view(count, 1:N))
    retvalsa_c = fill(Ptr{UInt8}(0), length(retvalsa))
    nc_get_vara_string(v.ncid, v.varid, start, count, retvalsa_c)
    for i = 1:length(retvalsa)
        retvalsa[i] = unsafe_string(retvalsa_c[i])
    end
    nc_free_string(length(retvalsa_c), retvalsa_c)
    retvalsa
end

function nc_get_var1_x(
    v::NcVar{String,N,NC_STRING},
    start::Vector{UInt},
    ::Type{String},
) where {N}
    retvalsa_c = fill(Ptr{UInt8}(0), 1)
    nc_get_var1_string(v.ncid, v.varid, start, retvalsa_c)
    retval = unsafe_string(retvalsa_c[1])
    nc_free_string(1, retvalsa_c)
    retval
end

function putatt(ncid::Integer, varid::Integer, atts::Dict)
    for a in atts
        name = a[1]
        val = a[2]
        nc_put_att(ncid, varid, name, val)
    end
end

function putatt(nc::NcFile, atts::Dict)
    putatt(nc.ncid, NC_GLOBAL, atts)
    merge!(nc.gatts, getatts(nc.ncid, NC_GLOBAL, collect(keys(atts))))
end

function putatt(var::NcVar, atts::Dict)
    putatt(var.ncid, var.varid, atts)
    merge!(var.atts, getatts(var.ncid, var.varid, collect(keys(atts))))
end

"""
    NetCDF.putatt(nc::NcFile, varname::String, atts::Dict)

Writes the attributes defined in `atts` to the variable `varname` in the NetCDF file
`nc`. Existing attributes are overwritten. If varname is not a valid variable name,
a global attribute will be written.
"""
function putatt(nc::NcFile, varname::AbstractString, atts::Dict)
    chdef = false
    if !nc.in_def_mode
        chdef = true
        nc_redef(nc.ncid)
    end
    if haskey(nc.vars, varname)
        putatt(nc.vars[varname], atts)
    else
        putatt(nc, atts)
    end
    chdef && nc_enddef(nc.ncid)
end

"""
    ncputatt(nc::String,varname::String,atts::Dict)

Writes the attributes defined in `atts` to the variable `varname` for the given NetCDF file name
`nc`. Existing attributes are overwritten. If varname is not a valid variable name,
a global attribute will be written.
"""
function ncputatt(ncs::AbstractString, varname::AbstractString, atts::Dict)
    open(ncs, mode = NC_WRITE) do nc
        if (nc.omode == NC_NOWRITE)
            fil = nc.name
            close(nc)
            println("reopening file in WRITE mode")
            open(fil, mode = NC_WRITE)
        end
        putatt(nc, varname, atts)
    end
end

putvar(
    nc::NcFile,
    varname::AbstractString,
    vals::AbstractArray;
    start = ones(Int, length(size(vals))),
    count = [size(vals)...],
) = putvar(nc[varname], vals, start = start, count = count)

"""
    NetCDF.putvar(v::NcVar,vals::Array;start::Vector=ones(Int,length(size(vals))),count::Vector=[size(vals)...])

Writes the values from the array `vals` to a NetCDF file. `v` is the `NcVar` handle of the respective variable and `vals` an array
with the same dimension as the variable in the NetCDF file.

### Keyword arguments

* `start` Vector of length `ndim(v)` setting the starting index for each dimension
* `count` Vector of length `ndim(v)` setting the count of values to be read along each dimension. The value -1 is treated as a special case to read all values from this dimension

"""
function putvar(
    v::NcVar,
    vals::AbstractArray;
    start::Vector = ones(Int, length(size(vals))),
    count::Vector = [size(vals)...],
)
    isa(vals, Array) ||
    Base.iscontiguous(vals) || error("Can only write from contiguous pieces of memory")
    p = preparestartcount(start, count, v)
    nc_put_vara_x(v, gstart, gcount, vals)
end

"""
    NetCDF.putvar(v::NcVar, val, i...)

Writes the value(s) `val` to the variable `v` while the indices are given in in an array-style indexing manner.
"""
@generated function putvar(v::NcVar{T,N}, val::Any, I::IndR...) where {T,N}

    N == length(I) || error("Dimension mismatch")
    quote
        isa(val, Number) ||
        isa(val, Array) ||
        Base.iscontiguous(val) || error("Can only write from contiguous pieces of memory")
        @nexprs $N i -> gstart[v.ndim+1-i] = firsti(I[i], v.dim[i].dimlen)
        @nexprs $N i -> gcount[v.ndim+1-i] = counti(I[i], v.dim[i].dimlen)
        checkboundsNC(v)
        p = 1
        @nexprs $N i -> p = p * gcount[v.ndim+1-i]
        length(val) != p && error(string(
            "Size of output array ($(length(retvalsa))) does not equal number of elements to be read (",
            p,
            ")!",
        ))
        nc_put_vara_x(v, gstart, gcount, val)
    end
end

@generated function putvar(v::NcVar{T,N}, val::Any, I::Integer...) where {T,N}

    N == length(I) || error("Dimension mismatch")
    if N == 0
        quote
            gstart[1] = 0
            nc_put_var1_x(v, gstart, val[1])
        end
    else
        quote
            @nexprs $N i -> gstart[v.ndim+1-i] = I[i] - 1
            @nall(
                $N,
                d -> ((I[d] <= v.dim[d].dimlen && I[d] > 0) || v.dim[d].unlim),
            ) || throw(BoundsError(v, I))
            nc_put_var1_x(v, gstart, val[1])
        end
    end

end

for (t, ending, arname) in funext
    fname = Symbol("nc_put_vara_$ending")
    fname1 = Symbol("nc_put_var1_$ending")
    arsym = Symbol(arname)
    @eval nc_put_vara_x(
        v::NcVar,
        start::Vector{UInt},
        count::Vector{UInt},
        vals::AbstractArray{$t},
    ) = $fname(v.ncid, v.varid, start, count, vals)
    @eval nc_put_var1_x(v::NcVar, start::Vector{UInt}, val::$t) = begin
        $(arsym)[1] = val
        $fname1(v.ncid, v.varid, start, $(arsym))
    end
end

nc_put_vara_x(
    v::NcVar{ASCIIChar,N,NC_CHAR},
    start::Vector{UInt},
    count::Vector{UInt},
    vals::AbstractArray{ASCIIChar},
) where {N} = nc_put_vara_text(v.ncid, v.varid, start, count, vals)
nc_put_var1_x(
    v::NcVar{ASCIIChar,N,NC_CHAR},
    start::Vector{UInt},
    val::ASCIIChar,
) where {N} = begin
    vals = UInt8[val]
    nc_put_var1_text(v.ncid, v.varid, start, count, vals)
end

function nc_put_vara_x(
    v::NcVar{String,N,NC_STRING},
    start,
    count,
    vals::AbstractArray{String},
) where {N}
    vals_p = map(x -> pointer(x), vals)
    nc_put_vara_string(v.ncid, v.varid, start, count, vals_p)
end

function nc_put_var1_x(
    v::NcVar{String,N,NC_STRING},
    start::Vector{UInt},
    val::String,
) where {N}
    val_p = [pointer(val)]
    nc_put_var1_string(v.ncid, v.varid, start, val_p)
end

function Base.push!(v::NcVar, a::AbstractArray)
    sold = size(v)
    N = ndims(v)
    iunlim = findall(map(x -> x.unlim, v.dim))
    length(iunlim) == 1 || error("You can only push to a NetCDF variable with one unlimited dimension")
    st = fill(1, N)
    st[iunlim[1]] = sold[iunlim[1]] + 1
    co = fill(-1, N)
    if ndims(v) == ndims(a)
        co[iunlim[1]] = size(a, iunlim[1])
    elseif ndims(v) == ndims(a) + 1
        co[iunlim[1]] = 1
    else
        error("You can only push variables that have equal or one fewer dimension than the NetCDF Variable")
    end
    NetCDF.putvar(v, a, start = st, count = co)
end

function Base.push!(v::NcVar{T,1}, a::Number) where {T}
    push!(v, collect(a))
end

"""
    NetCDF.sync(nc::NcFile)

Synchronizes the changes made to the file and writes changes to the disk.
"""
function sync(nc::NcFile)
    nc_sync(nc.ncid)
end

function setcompression(v::NcVar, mode)
    if v.compress > -1
        if (NC_NETCDF4 & mode) == 0
            warn("Compression only possible for NetCDF4 file format. Compression will be ingored.")
            v.compress = -1
        else
            v.compress = min(v.compress, 9)
            nc_def_var_deflate(v.ncid, v.varid, Int32(1), Int32(1), v.compress)
        end
    end
end

"""
    NetCDF.create(name::String,varlist::Array{NcVar};gatts::Dict=Dict{Any,Any}(),mode::UInt16=NC_NETCDF4)

Creates a new NetCDF file. Here, `name`
 is the name of the file to be created and `varlist` an array of `NcVar` holding the variables that should appear in the file.

### Keyword arguments

* `gatts` a Dict containing global attributes of the NetCDF file
* `mode` NetCDF file type (`NC_NETCDF4`, `NC_CLASSIC_MODEL` or `NC_64BIT_OFFSET`), defaults to `NC_NETCDF4`
"""
function create(
    name::AbstractString,
    varlist::Array{<:NcVar};
    gatts::Dict = Dict{Any,Any}(),
    mode::UInt16 = NC_NETCDF4,
    add_finalizer = true,
)

    #Create the file
    id = nc_create(name, mode)
    # Collect Dimensions and set NetCDF ID
    vars = Dict{String,NcVar}()
    varlist = collect(NcVar, varlist)
    dims = Set{NcDim}()
    for v in varlist
        v.ncid = id
        for d in v.dim
            push!(dims, d)
        end
    end
    nunlim = 0
    ndim = Int32(length(dims))

    #Create the NcFile Object
    nc = NcFile(
        id,
        Int32(length(vars)),
        ndim,
        zero(Int32),
        vars,
        Dict{String,NcDim}(),
        Dict{Any,Any}(),
        Int32[],
        name,
        NC_WRITE,
        true,
    )

    for d in dims
        create_dim(nc, d)
        if d.unlim
            push!(nc.unlimdimids, d.dimid)
        end
        if (length(d.vals) > 0) & (!haskey(nc.vars, d.name))
            elt = eltype(d.vals)
            push!(
                varlist,
                NcVar{elt,1,jl2nc(elt)}(
                    id,
                    varida[1],
                    1,
                    length(d.atts),
                    jl2nc(elt),
                    d.name,
                    [d.dimid],
                    [d],
                    d.atts,
                    -1,
                    (zero(Int32),),
                    nc
                ),
            )
        end
    end

    # Create variables in the file
    for v in varlist
        create_var(nc, v, mode)
        v.ncfile = nc
    end

    # Put global attributes
    if !isempty(gatts)
        putatt(nc, gatts)
    end

    # Leave define mode
    nc_enddef(nc)

    for d in nc.dim
        #Write dimension variable
        if length(d[2].vals) > 0
            putvar(nc, d[2].name, d[2].vals)
        end
    end

    add_finalizer && finalizer(_close,nc)

    return nc
end

create(
    name::AbstractString,
    varlist::NcVar...;
    gatts::Dict = Dict{Any,Any}(),
    mode::UInt16 = NC_NETCDF4,
    add_finalizer = true,
) = create(name, NcVar[varlist[i] for i = 1:length(varlist)]; gatts = gatts, mode = mode, add_finalizer = add_finalizer)

"""
    NetCDF._close(nc::NcFile)

closes a NetCDF file handle
"""
function _close(nco::Union{NcFile, NcVar})
    #Close file
    nc_close(nco.ncid)
    nothing
end

function close(x::Union{NcFile, NcVar})
  println(stacktrace())
  @warn "NetCDF.close is deprecated, because closing files is done with finalizers from now on. "
end

"""
    NetCDF.open(fil::AbstractString,v::AbstractString)

opens a NetCDF variable `v` in the NetCDF file `fil` and returns an `NcVar` handle that implements the `AbstractArray` interface for reading and writing.
Note that it is in the user's responsibility to close the file after usage using `NetCDF.close`.

### Keyword arguments

* `mode` mode in which the file is opened, defaults to `NC_NOWRITE`, choose `NC_WRITE` for write access
* `readdimvar` determines if dimension variables will be read into the file structure, default is `false`
"""
function open(
    fil::AbstractString,
    v::AbstractString;
    mode::Integer = NC_NOWRITE,
    readdimvar::Bool = false,
    add_finalizer = true,
)
    nc = open(fil, mode = mode, readdimvar = readdimvar, add_finalizer=add_finalizer)
    nc.vars[v]
end

"""
    NetCDF.open(fil::AbstractString)

opens the NetCDF file `fil` and returns a `NcFile` handle. Note that it is in the user's responsibility to close the file after usage using `NetCDF.close`.

### Keyword arguments

* `mode` mode in which the file is opened, defaults to `NC_NOWRITE`, choose `NC_WRITE` for write access
* `readdimvar` determines if dimension variables will be read into the file structure, default is `false`
"""
function open(fil::AbstractString; mode::Integer = NC_NOWRITE, readdimvar::Bool = false, add_finalizer=true)

    # Open netcdf file
    ncid = nc_open(fil, mode)

    #Get initial information
    ndim, nvar, ngatt, nunlimdimid = nc_inq(ncid)

    #Get unlimited dimensions
    nunlimdims = Ref(Int32(0))
    nc_inq_unlimdims(ncid, nunlimdims, C_NULL)
    unlimdims = fill(Int32(0), nunlimdims[])
    nc_inq_unlimdims(ncid, nunlimdims, unlimdims)

    #Create ncdf object
    ncf = NcFile(
        ncid,
        Int32(nvar - ndim),
        ndim,
        ngatt,
        Dict{String,NcVar}(),
        Dict{String,NcDim}(),
        Dict{Any,Any}(),
        unlimdims,
        abspath(fil),
        mode,
        false,
    )

    #Read global attributes
    ncf.gatts = getatts_all(ncid, NC_GLOBAL, ngatt)

    #Read dimensions
    for dimid = 0:ndim-1
        (name, dimlen) = nc_inq_dim(ncid, dimid)
        ncf.dim[name] = NcDim(
            ncid,
            dimid,
            -1,
            name,
            dimlen,
            [],
            Dict{Any,Any}(),
            in(dimid, unlimdims) ? true : false,
        )
    end

    #Read variable information
    for varid = 0:(nvar-1)
        (name, nctype, dimids, natts, vndim, isdimvar, chunksize) = nc_inq_var(ncf, varid)
        if (isdimvar)
            ncf.dim[name].varid = varid
        end
        atts = getatts_all(ncid, varid, natts)
        vdim = Array{NcDim}(undef, length(dimids))
        for (i, did) in enumerate(dimids)
            vdim[i] = ncf.dim[getdimnamebyid(ncf, did)]
        end
        ncf.vars[name] = NcVar{nctype2jltype[nctype],Int(vndim),Int(nctype)}(
            ncid,
            Int32(varid),
            vndim,
            natts,
            nctype,
            name,
            dimids[vndim:-1:1],
            vdim[vndim:-1:1],
            atts,
            0,
            chunksize,
            ncf,
        )
    end
    readdimvar == true && _readdimvars(ncf)
    add_finalizer && finalizer(_close,ncf)
    return ncf
end

"""
    NetCDF.open(f::Function, args...;kwargs...)

Opens a NetCDF file, applies the function f on the resulting file or variable handle and properly closes the file after usage.
This is convenient to use together with the `do` block syntax, for example:

data = open("myfile.nc","myvar") do v
  NetCDF.readvar(v,start=[1,1,1], count=[-1,-1,1])
end
"""
function open(f::Function, args...; kwargs...)
    io = open(args...;add_finalizer=false, kwargs...)
    try
        f(io)
    finally
        _close(io)
    end
end

"""
    NetCDF.create(f::Function, args...;kwargs...)

Creates a NetCDF file, applies the function f on the resulting file or variable handle and properly closes the file after usage.
This is convenient to use together with the `do` block syntax, for example:

d = NcDim("time",1:10)
v = NcVar("obs",d);
NetCDF.create("newfile.nc",v) do nc
  nc["obs"][:] = rand(10)
end
"""
function create(f::Function, args...; kwargs...)
    io = create(args...; add_finalizer=false, kwargs...)
    try
        f(io)
    finally
        _close(io)
    end
end

# Define some high-level functions
"""
    ncread(filename, varname)

reads the values of the variable varname from file filename and returns the values in an array.

### Keyword arguments

* `start` Vector of length `ndim(v)` setting the starting index for each dimension
* `count` Vector of length `ndim(v)` setting the count of values to be read along each dimension. The value -1 is treated as a special case to read all values from this dimension

### Example

To read the second slice of a 3D NetCDF variable one can write:

    ncread("filename","varname", start=[1,1,2], count = [-1,-1,1])

"""
function ncread(
    fil::AbstractString,
    vname::AbstractString;
    start::Array{T} = Array{Int}(undef, 0),
    count::Array{T} = Array{Int}(undef, 0),
) where {T<:Integer}
    open(fil) do nc
        length(start) == 0 && (start = defaultstart(nc[vname]))
        length(count) == 0 && (count = defaultcount(nc[vname]))
        readvar(nc[vname], start = start, count = count)
    end
end
ncread(
    fil::AbstractString,
    vname::AbstractString,
    start::Array{T,1},
    count::Array{T,1},
) where {T<:Integer} = ncread(fil, vname, start = start, count = count)

"""
    ncread!(filename, varname, d)

reads the values of the variable varname from file filename and writes the results to the pre-allocated array `d`

### Keyword arguments

* `start` Vector of length `ndim(v)` setting the starting index for each dimension
* `count` Vector of length `ndim(v)` setting the count of values to be read along each dimension. The value -1 is treated as a special case to read all values from this dimension

### Example

To read the second slice of a 3D NetCDF variable one can write:

    d = zeros(10,10,1)
    ncread!("filename","varname", d, start=[1,1,2], count = [-1,-1,1])

"""
function ncread!(
    fil::AbstractString,
    vname::AbstractString,
    vals::AbstractArray;
    start::Vector{Int} = ones(Int, ndims(vals)),
    count::Vector{Int} = [size(vals, i) for i = 1:ndims(vals)],
)
    open(fil) do nc
        readvar!(nc, vname, vals, start = start, count = count)
    end
end

"""
    ncinfo()

prints information on the variables, dimension and attributes conatained in the file
"""
function ncinfo(fil::AbstractString)
    open(show, fil)
end

iswritable(nc::NcFile) = (nc.omode & NC_WRITE) != zero(UInt16)

#High-level functions for writing data to a file
"""
    ncwrite(x::Array,fil::AbstractString,vname::AbstractString)

Writes the array `x` to the file `fil` and variable `vname`.

### Keyword arguments

* `start` Vector of length `ndim(v)` setting the starting index for writing for each dimension
* `count` Vector of length `ndim(v)` setting the count of values to be written along each dimension. The value -1 is treated as a special case to write all values from this dimension. This is usually inferred by the given array size.
"""
function ncwrite(
    x::Array,
    fil::AbstractString,
    vname::AbstractString;
    start = ones(Int, length(size(x))),
    count = [size(x)...],
)
    open(fil, mode = NC_WRITE) do nc
        if (nc.omode == NC_NOWRITE)
            close(nc)
            println("reopening file in WRITE mode")
            open(fil, mode = NC_WRITE)
        end
        putvar(nc, vname, x, start = start, count = count)
    end
end
ncwrite(x::Array, fil::AbstractString, vname::AbstractString, start::Array) =
    ncwrite(x, fil, vname, start = start)

"""
    ncgetatt(filename, varname, attname)

This reads a NetCDF attribute `attname` from the specified file and variable. To read global attributes, set varname to `Global`.
"""
function ncgetatt(fil::AbstractString, vname::AbstractString, att::AbstractString)
    open(fil) do nc
        haskey(nc.vars, vname) ? get(nc.vars[vname].atts, att, nothing) :
        get(nc.gatts, att, nothing)
    end
end

#High-level function for creating files and variables
#
# if the file does not exist, it will be created
# if the file already exists, the variable will be added to the file

function create_dim(nc, dim)
    nc_redef(nc)
    nc_def_dim(nc.ncid, dim.name, dim.dimlen, dima)
    dim.dimid = dima[1]
    nc.dim[dim.name] = dim
end

function create_var(nc, v, mode)
    nc_redef(nc)

    v.dimids = Int32[v.dim[i].dimid for i = 1:length(v.dim)]
    for i = 1:v.ndim
        dumids[i] = v.dimids[v.ndim+1-i]
    end
    nc_def_var(nc.ncid, v.name, v.nctype, v.ndim, dumids, vara)
    v.varid = vara[1]
    if any(i -> i > 0, v.chunksize)
        chunk_sizea[1:v.ndim] .= reverse(v.chunksize)
        nc_def_var_chunking(nc.ncid, v.varid, NC_CHUNKED, chunk_sizea)
    end
    nc.vars[v.name] = v
    putatt(nc.ncid, v.varid, v.atts)
    setcompression(v, mode)
end

"""
    nccreate (filename, varname, dimensions ...)

Create a variable in an existing NetCDF file or generates a new file. `filename` and `varname` are strings.
After that follows a list of dimensions. Each dimension entry starts with a dimension name (a String), and
may be followed by a dimension length, an array with dimension values or a Dict containing dimension attributes.
Then the next dimension is entered and so on. Have a look at examples/high.jl for an example use.

### Keyword arguments

- **atts** Dict of attribute names and values to be assigned to the variable created
- **gatts** Dict of attribute names and values to be written as global attributes
- **compress** Integer [0..9] setting the compression level of the file, only valid if `mode=NC_NETCDF4`
- **t** variable type, currently supported types are: const `NC_BYTE`, `NC_CHAR`, `NC_SHORT`, `NC_INT`, `NC_FLOAT`, `NC_LONG`, `NC_DOUBLE`
- **mode** file creation mode, only valid when new file is created, choose one of: `NC_NETCDF4`, `NC_CLASSIC_MODEL`, `NC_64BIT_OFFSET`
"""
function nccreate(
    fil::AbstractString,
    varname::AbstractString,
    dims...;
    atts::Dict = Dict{Any,Any}(),
    gatts::Dict = Dict{Any,Any}(),
    compress::Integer = -1,
    t::Union{DataType,Integer} = NC_DOUBLE,
    mode::UInt16 = NC_NETCDF4,
    chunksize = (0,),
)
  # Checking dims argument for correctness
    dim = parsedimargs(dims)
  # Check chunksize
    chunksize = chunksize[1] == 0 ? ntuple(i -> 0, length(dim)) : chunksize
  # create the NcVar object
    v = NcVar(varname, dim, atts = atts, compress = compress, t = t, chunksize = chunksize)
  # Test if the file already exists
    if isfile(fil)
        nc = open(fil, mode = NC_WRITE) do nc
            v.ncid = nc.ncid
            v.ncfile = nc
            haskey(
                nc.vars,
                varname,
            ) && error("Variable $varname already exists in file $fil")
      # Check if dimensions exist, if not, create
      # Remember if dimension was created

            dcreate = falses(length(dim))
            for i = 1:length(dim)
                if !haskey(nc.dim, dim[i].name)
                    create_dim(nc, dim[i])
                    v.dimids[i] = dim[i].dimid
                    if length(dim[i].vals) > 0
                        elt = eltype(dim[i].vals)
                        create_var(
                            nc,
                            NcVar{elt,1,jl2nc(elt)}(
                                nc.ncid,
                                0,
                                1,
                                length(dim[i].atts),
                                jl2nc(elt),
                                dim[i].name,
                                [dim[i].dimid],
                                [dim[i]],
                                dim[i].atts,
                                -1,
                                (0,),
                                nc,
                            ),
                            mode,
                        )
                    end
                    dcreate[i] = true
                else
                    v.dimids[i] = nc.dim[dim[i].name].dimid
                    v.dim[i] = nc.dim[dim[i].name]
                    dcreate[i] = false
                end
            end

      # Create variable
            create_var(nc, v, mode)
            nc_enddef(nc)
            for i = 1:length(dim)
                if dcreate[i] & !isempty(dim[i].vals)
                    ncwrite(dim[i].vals, fil, dim[i].name)
                end
            end
        end
    else
        create(fil, v, gatts = gatts, mode = mode | NC_NOCLOBBER) do nc
          for d in dim
            !isempty(d.vals) && ncwrite(d.vals, fil, d.name)
          end
        end
    end
    return nothing
end

#show{T<:Any,N}(io::IO,a::NcVar{T,N})=println(io,a.name)
#showcompact{T<:Any,N}(io::IO,a::NcVar{T,N})=println(io,a.name)
function show(io::IO, nc::NcFile)
    nrow, ncol = Base.displaysize(io)
    hline = repeat("-", ncol)
    l1 = div(ncol, 3)
    l2 = 2 * l1
    println(io, "")
    println(io, "##### NetCDF File #####")
    println(io, "")
    println(io, nc.name)
    println(io, "")
    println(io, "##### Dimensions #####")
    println(io, "")
    println(io, tolen("Name", l2), tolen("Length", l1))
    println(io, hline)
    for d in nc.dim
        println(
            io,
            tolen(d[2].name, l2),
            tolen(
                d[2].unlim ? string("UNLIMITED (", d[2].dimlen, " currently)") :
                d[2].dimlen,
                l1,
            ),
        )
    end
    l1 = div(ncol, 5)
    l2 = 2 * l1
    println(io, "")
    println(io, "##### Variables #####")
    println(io, "")
    println(io, tolen("Name", l2), tolen("Type", l1), tolen("Dimensions", l2))
    println(io, hline)
    for v in nc.vars
        s1 = string(tolen(v[2].name, l2))
        s2 = string(tolen(nctype2string[Int(v[2].nctype)], l1))
        s3 = ""
        for d in v[2].dim
            s3 = string(s3, d.name, " ")
        end
        println(io, s1, s2, tolen(s3, l2))
    end
    l1 = div(ncol, 4)
    l2 = 2 * l1
    println(io, "")
    println(io, "##### Attributes #####")
    println(io, "")
    println(io, tolen("Variable", l1), tolen("Name", l1), tolen("Value", l2))
    println(io, hline)
    for a in nc.gatts
        println(io, tolen("global", l1), tolen(a[1], l1), tolen(a[2], l2))
    end
    for v in nc.vars
        for a in v[2].atts
            println(io, tolen(v[2].name, l1), tolen(a[1], l1), tolen(a[2], l2))
        end
    end
end

tolen(s::Any, l::Number) = tolen(string(s), round(Int, l))
function tolen(s::AbstractString, l::Integer)
    ls = length(s)
    if ls < l
        return rpad(s, l)
    elseif ls > l
        return string(s[1:prevind(s, l - 1)], "..")
    else
        return s
    end
end
function tolen(s::Array{String,1}, l::Number)
    cs = "["
    for se in s
        cs *= string(se, ", ")
    end
    cs = isempty(s) ? string(cs, ']') : string(cs[1:end-2], ']')
    return tolen(cs, l)
end

end # Module


text:
# Hicks-Moorsteen-Bjurek (HMB) total factor productivity index as defined in Bjurek (SJE, 1996)
immutable HMB{S<:AbstractDataEnvelopment,T<:RS}
  In0::Union{FDH,DEA}
  In1::Union{FDH,DEA}
  Out0::Union{FDH,DEA}
  Out1::Union{FDH,DEA}
  K::Int64

  function HMB{S,T}(X0::Array,Y0::Array,X1::Array,Y1::Array) where {S<:AbstractDataEnvelopment,T<:RS}
    if(size(X0) != size(X1) || size(Y0) != size(Y1))
      error("The number of DMUs must be the same in both periods!")
    end
    if(eltype(HMB{S,T})[1] == Convex)
      new(DEA{T}(X0,Y0,true),DEA{T}(X1,Y1,true),DEA{T}(X0,Y0,false),DEA{T}(X1,Y1,false),size(X0,1))
    elseif(eltype(HMB{S,T})[1] == FreeDisposal)
      new(FDH{T}(X0,Y0,true),FDH{T}(X1,Y1,true),FDH{T}(X0,Y0,false),FDH{T}(X1,Y1,false),size(X0,1))
    else
      # Not defined yet
    end
  end

end

# Return parametric types as a tuple
eltype{S,T}(::Type{HMB{S,T}}) = (S, T)

# Compute the HMB index
function (H::HMB)()
  MO0 = Array{Float64}(H.K)
  MO1 = Array{Float64}(H.K)
  MI0 = Array{Float64}(H.K)
  MI1 = Array{Float64}(H.K)
  for k=1:H.K
    X0,Y0 = getdata(H.Out0)[k]
    X1,Y1 = getdata(H.Out1)[k]

    MO0[k] = geteff(H.Out0(X0,Y1))/geteff(H.Out0(X0,Y0))
    MO1[k] = geteff(H.Out1(X1,Y1))/geteff(H.Out1(X1,Y0))

    X0,Y0 = getdata(H.In0)[k]
    X1,Y1 = getdata(H.In1)[k]

    MI0[k] = geteff(H.In0(X0,Y0))/geteff(H.In0(X1,Y0))
    MI1[k] = geteff(H.In1(X0,Y1))/geteff(H.In1(X1,Y1))
  end
  return sqrt((MO0./MI0).*(MO1./MI1))
end

# Technical efficiency change
function TEI(H::HMB)
  return geteff(H.Out1())./geteff(H.Out0())
end

# Technical change
function TC(H::HMB)
  TC0 = Array{Float64}(H.K)
  TC1 = Array{Float64}(H.K)
  for k=1:H.K
    X0,Y0 = getdata(H.Out0)[k]
    TC0[k] = geteff(H.Out1(X0,Y0))
    X1,Y1 = getdata(H.Out1)[k]
    TC1[k] = geteff(H.Out0(X1,Y1))
  end
  TC0 = geteff(H.Out0())./TC0
  TC1 = TC1./geteff(H.Out1())
  return real(sqrt(complex(TC0.*TC1)))
end

# Scale efficiency change
function SEC(H::HMB)
  return H()./(TC(H).*TEI(H))
end

max_stars_repo_name: kepiej/NEPA


text:
cd(@__DIR__); include("setups/grid7x3.jl");

## ROT1 + pmf2
D = natural_eigdist(𝚽, 𝛌, Q; α = 0.1, input_format = :pmf2, distance = :ROT)
E = transform(fit(MDS, D, maxoutdim=2, distances=true))
E[1, :] .*= sign(E[1, 2]); E[2, :] .*= -sign(E[1, 2])
plt = grid7x3_mds_heatmaps(E, 𝚽; backend = :pyplot)
savefig(plt, "../figs/Grid7x3_MDS_ROT1_pmf2_alpha01.png")
display(plt)


text:
## Load necessary code

using Pkg
Pkg.activate(@__DIR__)
pkg"instantiate"
pkg"precompile"

using FileIO

include(joinpath(@__DIR__, "CommonUtils.jl"))
# include(joinpath(@__DIR__, "Plotting.jl"))

## load specific case

pargs = if !isdefined(Main, :parsed_args)
  parse_commandline()
else
  parsed_args
end

pargs["reportDir"] = if !haskey(pargs, "reportDir") && isdefined(Main, :fg)
  getLogPath(fg)
else
  # "/tmp/caesar/2020-02-03T03:07:45.938/fg_after_x1181.tar.gz"
  # "/tmp/caesar/2020-02-22T02:21:20.777/fg_after_x781.tar.gz"
  # "/tmp/caesar/2020-02-23T01:43:32.222/fg_after_x1391.tar.gz"
  # "/tmp/caesar/2020-02-23T01:43:32.222/fg_final.tar.gz"
  # "/tmp/caesar/2020-02-23T01:43:32.222/fg_after_x1391.tar.gz"
  pargs["reportDir"]
end

# @show pargs["reportDir"]
# @show splitpath(pargs["reportDir"])

# load the factor graph needed
fg = if true || !isdefined(Main, :fg)
  println("going to load fg")
  fg = LightDFG{SolverParams}(params=SolverParams())
  @show pathElem = splitpath(pargs["reportDir"])
  @show getSolverParams(fg).logpath = joinpath(pathElem[1:end-1]...)
  loadDFG(pargs["reportDir"], Main, fg)
  fg
else
  fg
end

# random temp data
drt_data, TTm, XXm, YYm, XXf, YYf = loadResultsDRT(fg)

# load here due to name conflict with Point2
using MakieLayout
using Makie

include(joinpath(@__DIR__, "MakiePlotsFG.jl"))

##

DRT = true
PPE = true
REF = false

include(joinpath(@__DIR__, "MakiePlotImg.jl"))

scene
##

@show scenepath = joinLogPath(fg, "withBackround.png")
Base.rm(scenepath, force=true)
Makie.save(scenepath, scene)

##

using Glob

# pargs["plotSeriesBeliefs"] = 0

# plot a series of frames
if 0 <= pargs["plotSeriesBeliefsNew"]

global frame = 0 # still finiky, leave at 0
pattern = "fg_x"
ext = ".tar.gz"
files = glob("$(pattern)*$ext", getLogPath(fg))
indiv = splitpath.(files) .|> x->x[end]
presort = (x->x[length(pattern)+1:end-length(ext)]).(indiv)
sorted = parse.(Int, presort) |> sortperm
# reduce the number of frames, if requested
indiv = pargs["plotSeriesBeliefsNew"] == 0 ? indiv[sorted] : indiv[sorted[(frame+1):pargs["plotSeriesBeliefsNew"]]]
# get range from last file in glob
fg = LightDFG{SolverParams}(params=SolverParams(logpath=getLogPath(fg)))
loadDFG(joinLogPath(fg, indiv[end]), Main, fg)
# drt_data, TTm, XXm, YYm, XXf, YYf = loadResultsDRT(fg)
minmax = getRangeCartesian(fg,r"x\d",digits=-1)
# xmin=minmax[1,1]-10;xmax=minmax[1,2]+10;ymin=minmax[2,1]-10;ymax=minmax[2,2]+10;
Base.mkpath(joinLogPath(fg, "background"))
for ind in indiv
  global fg
  global frame += 1
  println("frame $frame of $(length(indiv))")
  fname = split(ind, '.')[1]
  if pargs["skip"] && isfile(joinLogPath(fg,"background/$fname.png"))
    @info "skip $fname since its already there"
    continue
  end
  fg = LightDFG{SolverParams}(params=SolverParams(logpath=getLogPath(fg)))
  loadDFG(joinLogPath(fg, ind), Main, fg)
  nvars = length(ls(fg, r"x\d"))
  # if nvars < 15
  #   continue
  # end
  try

    include(joinpath(@__DIR__, "MakiePlotImg.jl"))

    Base.rm(joinLogPath(fg,"background/$fname.png"), force=true)
    Makie.save(joinLogPath(fg,"background/$fname.png"), scene)
  catch ex
    @error ex
  end
  GC.gc()
end

end


text:
module Utils
# ジョルダン曲線: ねじれのない閉曲線のこと.
"""
ジョルダン閉曲線であるかどうか
"""
function is_jordan_curve(_gamma::Array)
    progression_of_differences = [sqrt((_gamma[i, 1] - _gamma[i+1, 1])^2 + (_gamma[i, 2] - _gamma[i+1, 2])^2) for i = 1:(length(_gamma[:, 1])-1)]
    ave_distance = sum(progression_of_differences) / length(progression_of_differences)
    if ave_distance * 2 < abs(_gamma[1, 1] - _gamma[length(_gamma[:, 1]), 1])
        return true
    else
        return false
    end
end

# FIXME: 倍々ではなく, linearに補完数を指定できるように。

"n行目とn+1行目の間のデータを補完,境界条件を付与 (x, 2)->(x*2-1, 2)"
function complement_p(array::Array, multiple::Bool, point_space::Float64)
    (x, y) = size(array)
    return_value = Array{Float64}(undef, 2 * x, y)
    if multiple
        for i = 1:x-1
            # 点の間隔が平均以下なので、補間する -> 同一の曲線とみなす
            if sqrt(sum((array[i, :] .- array[i+1, :]) .^ 2)) < point_space * 5 # NOTE: ここの加減が難しい. Nに依存?
                return_value[i*2-1, :] = array[i, :]
                return_value[i*2, :] = (array[i, :] .+ array[i+1, :]) ./ 2
                # 点の間隔が離れているので違う曲線とみる
            else
                return_value[i*2-1, :] = array[i, :]
                return_value[i*2, :] = array[i, :]
            end
        end
        return_value[x*2-1, :] = array[x, :]
        return_value[x*2, :] = array[x, :]
    else
        for i = 1:x-1
            # 点の間隔が平均以下なので、補間する -> 同一の曲線とみなす
            return_value[i*2-1, :] = array[i, :]
            return_value[i*2, :] = (array[i, :] .+ array[i+1, :]) ./ 2
        end
        return_value[x*2-1, :] = array[x, :]
        return_value[x*2, :] = array[1, :] # Note: _gamma += _gamma's head line coz boundary condition. size: (N+1,2)
    end
    return return_value
end
# precompile()

function remove_same_point(array::Any)
    return_value = Array{Any}(undef, 0, 2)
    array_length = length(array[:, 1])
    for i = 1:array_length-1
        if array[i, :] != array[i+1, :]
            return_value = vcat(return_value, array[i, :]')
        end
    end
    return_value = vcat(return_value, array[array_length, :]')
    # println(return_value)
    return return_value
end

"""
    times回, 倍の数だけデータを補完する
"""
function interpolation(array::Array, times::Int, multiple = false)
    if times > 0
        tmp = []
        (x, y) = size(array)
        point_space = sum([sqrt(sum((array[i, :] .- array[i+1, :]) .^ 2)) for i = 1:x-1]) / x
        for _ = 1:times
            tmp = complement_p(array, multiple, point_space)
            array = tmp
        end
        return remove_same_point(tmp)
    end
end
precompile(interpolation, (Array, Int, Bool))

export is_jordan_curve, interpolation, remove_same_point
end

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
using AvailablePotentialEnergyFramework
using Test
using NCDatasets: Dataset, variable
using Statistics: mean
using Unitful: @u_str, unit, ustrip, Quantity

# @test integrate_vertically(1:10,dz=2) == 110
# @test integrate_vertically(1:10,dz=2,weight=2) == 220
# testmat_2d= repeat(1:10,1,10)
# testweights_1d = reverse(1:10)
# testweights_2d = (1:10)'.*repeat(reverse(1:10),1,10)
# @test integrate_vertically(testmat_2d) == reshape(10*(1:10),10,1)
# @test integrate_vertically(testmat_2d,weight=testweights_1d) == reshape(10*(1:10)).*reverse((1:10),10,1)
# @test integrate_vertically(testmat_2d,weight=testweights_2d) == reshape(55*(1:10)).*reverse(1:10),10,1)
@testset "AvailablePotentialEnergyFramework" begin

    include("compositehelperfunctions.jl")
    include("physicsfunctions.jl")
    include("useful_diagnostics.jl")
    include("apebudgets.jl")
#    include("ape_computation_from_julia_output.jl")
    #include("arrayfiltering.jl")

end


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

@doc raw"""FieldsV1 stores a set of fields in a data structure like a Trie, in JSON format.  Each key is either a '.' representing the field itself, and will always map to an empty set, or a string representing a sub-field or item. The string will follow one of these four formats: 'f:<name>', where <name> is the name of a field in a struct, or key in a map 'v:<value>', where <value> is the exact json formatted value of a list item 'i:<index>', where <index> is position of a item in a list 'k:<keys>', where <keys> is a map of  a list item's key fields to their unique values If a key maps to an empty Fields value, the field that key represents is part of the set.  The exact format is defined in sigs.k8s.io/structured-merge-diff

    IoK8sApimachineryPkgApisMetaV1FieldsV1(;
    )

"""
mutable struct IoK8sApimachineryPkgApisMetaV1FieldsV1 <: SwaggerModel

    function IoK8sApimachineryPkgApisMetaV1FieldsV1(;)
        o = new()
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1FieldsV1

const _property_map_IoK8sApimachineryPkgApisMetaV1FieldsV1 = Dict{Symbol,Symbol}()
const _property_types_IoK8sApimachineryPkgApisMetaV1FieldsV1 = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1FieldsV1 }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1FieldsV1))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1FieldsV1 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1FieldsV1[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1FieldsV1 }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1FieldsV1[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1FieldsV1)
    true
end

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


text:
using LinearAlgebra
using GPFlux
using Zygote
using Random
using GaussianProcesses
using Test

@testset "negative logliklihood" begin
	Random.seed!(4)
	d = 4
	n = 10
	X = rand(d, n)
	y = rand(n)

	# GaussianProcesses GP
	gp_zmean = MeanZero()
	gp_se_kernel = SE(4.0, 3.0)
	gp_lognoise = -1.0
	gp_gp = GP(X, y, gp_zmean, gp_se_kernel, gp_lognoise)
	PT = GaussianProcesses.FullCovariancePrecompute(n)
	CST = GaussianProcesses.FullCovariance()
	GaussianProcesses.update_mll_and_dmll!(gp_gp, PT)
	gp_nll = -gp_gp.mll
	gp_grads = -gp_gp.dmll

	# GPFlux GP
	zmean = ConstantMean()
	se_kernel = IsoGaussKernel([4.0], [3.0])
	noise2 = [-1.0]
	gp = GaussProcess(zmean, se_kernel, noise2)
	ps = params(gp)
	my_nll = negloglik(gp, X, y)
	gs = gradient(()->negloglik(gp, X, y), ps)
	my_grads = [gs[p][1] for p in ps if p!=[0.0]]

	@show gp_nll
	@show my_nll
	@show sort(gp_grads)
	@show sort(my_grads)

	@test isapprox(my_nll, gp_nll, rtol=1e-3, atol=1e-8)
	@test isapprox(sort(gp_grads), sort(my_grads), rtol=1e-4, atol=1e-8)
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
include("GLNS.jl")

using  DataFrames, CSV
DataFrames_ = DataFrames;
function mkdire_(dire_)
	if !isdir(dire_)
		mkdir(dir_)
	end
end

ranges = 20:5:60
kk = 5
mm = 3
# print("k ", kk, " m ", mm,"\n")
# k_m = "k_" * string(kk) * "_m_" * string(mm)

itera = 5
ins_folder = "GLNSLIB_RANDOM_TPP/"

function runSimulation()
    scale = 1000
    for num_clusters = ranges
        for iter = 1:2
            file_name = "F:/Acad/research/Carlsson/GLNS/TPP_results/result_$num_clusters" * "_" * "$iter.txt"
            instance_name = "GLNSLIB/input_" * string(num_clusters) * "_" * string(iter) * ".gtsp"
            ret = GLNS.solver("$instance_name", mode = "slow")
            sol_obj = ret[1]/scale
            sol_time = ret[2]
            f = open(file_name, "w")
            write(f, "Optimal Objective: $sol_obj \r\n")
            write(f, "Solution Time: $sol_time \r\n")
            close(f)
        end
    end
end
# runSimulation()

function runProofTest()
    scale = 1000
    file_name = "result" * k_m * ".txt"
    f = open(file_name, "w")
    for num_clusters = ranges
        println("num_clusters = $num_clusters")
        sol_obj = 0
        divided = mm * sqrt(2 * num_clusters / kk)
        for iter = 1:itera
            instance_name = "PROOFLIB_" * k_m *"/input_" * string(num_clusters) * "_" * string(iter) * ".gtsp"
            ret = GLNS.solver("$instance_name", mode = "slow")
            sol_obj += ret[1]/scale
        end
        sol_obj = sol_obj / itera
        println("obj = $sol_obj")
        retVal = sol_obj / divided
        println("normed_obj = $retVal")
        write(f, "$retVal,")
    end
    close(f)
end
# runProofTest()

function write_res(ins, n, k, m, seed_, obj)

	df = DataFrames_.DataFrame(ins=ins, n=n,k=k,m=m, seed_=seed_, obj=obj)
    CSV.write("res.csv", df, append=true)
end

function runGTSPApx(ins_folder, ins_name)
    scale = 1000
	ins_na = "\n\n\n"*ins_name
	println(ins_na)

	ins_info = split(ins_name,"_")
	nn = parse(Int, ins_info[end-3])
	kk = parse(Int, ins_info[end-2])
	mm = parse(Int, ins_info[end-1])
	seed_ = parse(Int, split(ins_info[end],".")[1])
	print("n $nn, k $kk, m $mm , seed $seed_")
	sol_obj = 0
	divided = mm * sqrt(2 * nn / kk)
	for iter = 1:itera
		ret = GLNS.solver(ins_folder* ins_name, mode = "fast")
		sol_obj += ret[1]/scale
	end
	sol_obj = sol_obj / itera
	println("obj = $sol_obj \n")
	retVal = sol_obj / divided
	println("normed_obj = $retVal \n")
	write_res(ins_name, nn,kk,mm,seed_, sol_obj)

end

if size(ARGS)[1]>0
	ins_folder = ARGS[1]
	ins_name = ARGS[2]
	runGTSPApx(ins_folder, ins_name)
	exit()
end

function runGTSPApx_folder()
    for ins in readdir(ins_folder)
		runGTSPApx(ins_folder, ins)
    end
end
runGTSPApx_folder()

max_stars_repo_name: Azizimj/GTSPbnd

max_issues_repo_name: Azizimj/GTSPbnd

max_forks_repo_name: Azizimj/GTSPbnd