text:
# Helper types

struct PVTKArgs
    part::Int
    nparts::Int
    ismain::Bool
    ghost_level::Int
end

struct PVTKFile <: VTKFile
    pvtkargs::PVTKArgs
    xdoc::XMLDocument
    vtk::DatasetFile
    path::String
end

# Returns true if the arguments do *not* contain any cell vectors.
_pvtk_is_structured(x::AbstractVector{<:AbstractMeshCell}, args...) = Val(false)
_pvtk_is_structured(x, args...) = _pvtk_is_structured(args...)
_pvtk_is_structured() = Val(true)

_pvtk_nparts(structured::Val{false}; nparts::Integer, etc...) = nparts
_pvtk_nparts(structured::Val{true}; extents::AbstractArray, etc...) = length(extents)

_pvtk_extents(structured::Val{false}; etc...) = nothing
_pvtk_extents(structured::Val{true}; extents::AbstractArray, etc...) = extents

# Filter keyword arguments to be passed to vtk_grid.
_remove_parallel_kwargs(; nparts = nothing, extents = nothing, kws...) =
    NamedTuple(kws)

# Determine whole extent from the local extents associated to each process.
function compute_whole_extent(extents::AbstractArray{<:Extent})
    # Compute the minimum and maximum across each dimension
    begins = mapreduce(ext -> map(first, ext), min, extents)
    ends = mapreduce(ext -> map(last, ext), max, extents)
    map((b, e) -> b:e, begins, ends)
end

compute_whole_extent(::Nothing) = nothing

"""
    pvtk_grid(
        filename, args...;
        part, nparts, extents, ismain = (part == 1), ghost_level = 0,
        kwargs...,
    )

Returns a handler representing a parallel VTK file, which can be
eventually written to file with `vtk_save`.

Positional and keyword arguments in `args` and `kwargs` are passed to `vtk_grid`
verbatim.
Note that serial filenames are automatically generated from `filename` and from
the process id `part`.

The following keyword arguments only apply to parallel VTK file formats.

Mandatory ones are:

- `part`: current (1-based) part id,
- `nparts`: total number of parts (only required for **unstructured** grids),
- `extents`: array specifying the partitioning of a **structured** grid across
  different processes (see below for details).

Optional ones are:

- `ismain`: `true` if the current part id `part` is the main (the only one that will write the `.pvtk` file),
- `ghost_level`: ghost level.

## Specifying extents for a structured grid

For structured grids, the partitioning of the dataset across different processes
must be specified via the `extents` argument.
This is an array where each element represents the data extent associated to a
given process.

For example, for a dataset of global dimensions ``15×12×4`` distributed across 4
processes, this array may look like the following:

```julia
extents = [
    ( 1:10,  1:5, 1:4),  # process 1
    (10:15,  1:5, 1:4),  # process 2
    ( 1:10, 5:12, 1:4),  # process 3
    (10:15, 5:12, 1:4),  # process 4
]
```

Some important notes:

- the `extents` argument must be the same for all processes;
- the extents **must overlap**, or VTK / ParaView will complain when trying to
  open the files;
- the length of the `extents` array gives the number of processes.
  Therefore, the `nparts` argument is redundant and does not need to be passed.

"""
function pvtk_grid(
        filename::AbstractString, args...;
        part, ismain = (part == 1), ghost_level = 0, kwargs...,
    )
    is_structured = _pvtk_is_structured(args...)
    nparts = _pvtk_nparts(is_structured; kwargs...)
    extents = _pvtk_extents(is_structured; kwargs...)

    mkpath(filename)
    bname = basename(filename)
    prefix = joinpath(filename, bname)
    fn = _serial_filename(part, nparts, prefix, "")

    vtk = let kws_vtk = _remove_parallel_kwargs(; kwargs...)
        kws = if extents === nothing
            kws_vtk
        else
            (; kws_vtk..., extent = extents[part])
        end
        vtk_grid(fn, args...; kws...)
    end

    pvtkargs = PVTKArgs(part, nparts, ismain, ghost_level)
    xdoc = XMLDocument()
    _, ext = splitext(vtk.path)
    path = filename * ".p" * ext[2:end]
    pvtk = PVTKFile(pvtkargs, xdoc, vtk, path)
    _init_pvtk!(pvtk, extents)

    pvtk
end

# Add point and cell data as usual
function Base.setindex!(
        pvtk::PVTKFile, data, name::AbstractString, args...,
    )
    pvtk.vtk[name, args...] = data
end

# Save as usual
function vtk_save(pvtk::PVTKFile)
    outfiles = String[]
    if isopen(pvtk)
        if pvtk.pvtkargs.ismain
            _update_pvtk!(pvtk)
            save_file(pvtk.xdoc, pvtk.path)
            push!(outfiles, pvtk.path)
        end
        append!(outfiles, vtk_save(pvtk.vtk))
        close(pvtk)
    end
    outfiles
end

# Helper functions
function _serial_filename(part, nparts, prefix, extension)
    p = lpad(part, ndigits(nparts), '0')
    prefix * "_$p" * extension
end

function _init_pvtk!(pvtk::PVTKFile, extents)
    # Recover some data
    vtk = pvtk.vtk
    pvtkargs = pvtk.pvtkargs
    pgrid_type = "P" * vtk.grid_type
    npieces = pvtkargs.nparts
    pref, _ = splitext(pvtk.path)
    _, ext = splitext(vtk.path)
    prefix = joinpath(pref, basename(pref))

    # VTKFile (root) node
    pvtk_root = create_root(pvtk.xdoc, "VTKFile")
    set_attribute(pvtk_root, "type", pgrid_type)
    set_attribute(pvtk_root, "version", "1.0")
    if IS_LITTLE_ENDIAN
        set_attribute(pvtk_root, "byte_order", "LittleEndian")
    else
        set_attribute(pvtk_root, "byte_order", "BigEndian")
    end

    # Grid node
    pvtk_grid = new_child(pvtk_root, pgrid_type)
    set_attribute(pvtk_grid, "GhostLevel", string(pvtkargs.ghost_level))

    # Pieces (i.e. Pointers to serial files)
    for piece ∈ 1:npieces
        pvtk_piece = new_child(pvtk_grid, "Piece")
        fn = _serial_filename(piece, npieces, prefix, ext)
        set_attribute(pvtk_piece, "Source", fn)

        # Add local extent if necessary
        if extents !== nothing
            set_attribute(pvtk_piece, "Extent", extent_attribute(extents[piece]))
        end
    end

    # Add whole extent if necessary
    whole_extent = compute_whole_extent(extents)
    if whole_extent !== nothing
        set_attribute(pvtk_grid, "WholeExtent", extent_attribute(whole_extent))
    end

    # Getting original grid informations
    # Recover point type and number of components
    vtk_root = root(vtk.xdoc)
    vtk_grid = find_element(vtk_root, vtk.grid_type)

    # adding origin and spacing if necessary
    origin = attribute(vtk_grid, "Origin")
    if origin !== nothing
        set_attribute(pvtk_grid, "Origin", origin)
    end

    spacing = attribute(vtk_grid, "Spacing")
    if spacing !== nothing
        set_attribute(pvtk_grid, "Spacing", spacing)
    end

    # Getting original piece informations
    vtk_piece = find_element(vtk_grid, "Piece")

    # If serial VTK has points
    vtk_points = find_element(vtk_piece, "Points")
    if vtk_points !== nothing
        vtk_data_array = find_element(vtk_points, "DataArray")
        point_type = attribute(vtk_data_array, "type")
        Nc = attribute(vtk_data_array, "NumberOfComponents")
        ## PPoints node
        pvtk_ppoints = new_child(pvtk_grid, "PPoints")
        pvtk_pdata_array = new_child(pvtk_ppoints, "PDataArray")
        set_attribute(pvtk_pdata_array, "type", point_type)
        set_attribute(pvtk_pdata_array, "Name", "Points")
        set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)
    end

    # If serial VTK has coordinates
    vtk_coordinates = find_element(vtk_piece, "Coordinates")
    if vtk_coordinates !== nothing
        pvtk_pcoordinates = new_child(pvtk_grid, "PCoordinates")
        for c in child_elements(vtk_coordinates)
            pvtk_pdata_array = new_child(pvtk_pcoordinates, "PDataArray")
            set_attribute(pvtk_pdata_array, "type", attribute(c, "type"))
            set_attribute(pvtk_pdata_array, "Name", attribute(c, "Name"))
            set_attribute(pvtk_pdata_array, "NumberOfComponents", attribute(c, "NumberOfComponents"))
        end
    end

    pvtk
end

function _update_pvtk!(pvtk::PVTKFile)
    vtk = pvtk.vtk
    vtk_root = root(vtk.xdoc)
    vtk_grid = find_element(vtk_root, vtk.grid_type)
    vtk_piece = find_element(vtk_grid, "Piece")

    pgrid_type = "P" * vtk.grid_type
    pvtk_root = root(pvtk.xdoc)
    pvtk_grid = find_element(pvtk_root, pgrid_type)

    # Generate PPointData
    vtk_point_data = find_element(vtk_piece, "PointData")
    if vtk_point_data !== nothing
        pvtk_ppoint_data = new_child(pvtk_grid, "PPointData")
        for vtk_data_array in child_elements(vtk_point_data)
            t = attribute(vtk_data_array, "type")
            name = attribute(vtk_data_array, "Name")
            Nc = attribute(vtk_data_array, "NumberOfComponents")
            pvtk_pdata_array = new_child(pvtk_ppoint_data, "PDataArray")
            set_attribute(pvtk_pdata_array, "type", t)
            set_attribute(pvtk_pdata_array, "Name", name)
            set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)
        end
    end

    # Generate PCellData
    vtk_cell_data = find_element(vtk_piece, "CellData")
    if vtk_cell_data !== nothing
        pvtk_pcell_data = new_child(pvtk_grid, "PCellData")
        for vtk_data_array in child_elements(vtk_cell_data)
            t = attribute(vtk_data_array, "type")
            name = attribute(vtk_data_array, "Name")
            Nc = attribute(vtk_data_array, "NumberOfComponents")
            pvtk_pdata_array = new_child(pvtk_pcell_data, "PDataArray")
            set_attribute(pvtk_pdata_array, "type", t)
            set_attribute(pvtk_pdata_array, "Name", name)
            set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)
        end
    end

    pvtk
end


text:
# __precompile__()

module ClimateTools

# External modules
using NetCDF
using Reexport
@reexport using NCDatasets
using Shapefile
using AxisArrays
using NaNMath
const axes = Base.axes
using ArgCheck
using PyCall
using PyPlot
using Interpolations
using Images
using ProgressMeter
using Polynomials
using IterTools
using Statistics
using Dates
using GeoStats
using InverseDistanceWeighting
import Base.vcat
import Base.getindex
import Base.show
import Base.size
import Base.setindex!
import Base.similar
import Base.write
import Statistics.minimum
import Statistics.maximum
import Statistics.std
import Statistics.var
import Statistics.mean
import Base: +
import Base: -
import Base: *
import Base: /
import Base.findmax

const basemap = PyNULL()
const mpl = PyNULL()
const cmocean = PyNULL()
# const scipy = PyNULL()

function __init__()
  copy!(mpl, pyimport_conda("matplotlib", "matplotlib"))
  copy!(basemap, pyimport_conda("mpl_toolkits.basemap", "basemap"))
  copy!(cmocean, pyimport_conda("cmocean", "cmocean", "conda-forge"))
  # copy!(scipy, pyimport_conda("scipy.interpolate", "scipy"))
end

# TYPES

"""
    ClimGrid{A <: AxisArray}

In-memory representation of Climate Forecast netCDF files.

struct ClimGrid\n
  data::AxisArray # Data\n
  longrid::AbstractArray{N,2} where N # the longitude grid\n
  latgrid::AbstractArray{N,2} where N # the latitude grid\n
  msk::Array{N, 2} where N # Data mask (NaNs and 1.0)\n
  grid_mapping::Dict#{String, Any} # bindings for native grid\n
  dimension_dict::Dict\n
  model::String\n
  frequency::String # Day, month, years\n
  experiment::String # Historical, RCP4.5, RCP8.5, etc.\n
  run::String\n
  project::String # CORDEX, CMIP5, etc.\n
  institute::String # UQAM, DMI, etc.\n
  filename::String # Path of the original file\n
  dataunits::String # Celsius, kelvin, etc.\n
  latunits::String # latitude coordinate unit\n
  lonunits::String # longitude coordinate unit\n
  variable::String # Type of variable (i.e. can be the same as "typeofvar", but it is changed when calculating indices)\n
  typeofvar::String # Variable type (e.g. tasmax, tasmin, pr)\n
  typeofcal::String # Calendar type\n
  timeattrib::Dict # Time attributes (e.g. days since ... )\n
  varattribs::Dict # Variable attributes dictionary\n
  globalattribs::Dict # Global attributes dictionary\n
end\n
"""
struct ClimGrid{A <: AxisArray}
  data::A
  longrid::Array{N,T} where T where N
  latgrid::Array{N,T} where T where N
  msk::Array{N,T} where T where N
  grid_mapping::Dict # information of native grid
  dimension_dict::Dict
  timeattrib::Dict
  model::String
  frequency::String
  experiment::String
  run::String
  project::String
  institute::String
  filename::String
  dataunits::String
  latunits::String # of the coordinate variable
  lonunits::String # of the coordinate variable
  variable::String # Type of variable
  typeofvar::String # Variable type (e.g. tasmax, tasmin, pr)
  typeofcal::String # Calendar type
  varattribs::Dict # Variable attributes
  globalattribs::Dict # Global attributes

end

"""
    ClimGrid(data; longrid=[], latgrid=[], msk=[], grid_mapping=Dict(), dimension_dict=Dict(), model="NA", frequency="NA", experiment="NA", run="NA", project="NA", institute="NA", filename="NA", dataunits="NA", latunits="NA", lonunits="NA", variable="NA", typeofvar="NA", typeofcal="NA", varattribs=Dict(), globalattribs=Dict())

Constructor of the ClimGrid function. Data is an AxisArray. Everything else is optional, but usually needed for further processing (mapping, interpolation, etc...).

struct ClimGrid\n
  data::AxisArray # Data \n
  longrid::AbstractArray{N,2} where N # the longitude grid \n
  latgrid::AbstractArray{N,2} where N # the latitude grid \n
  msk::Array{N, 2} where N # Data mask (NaNs and 1.0) \n
  grid_mapping::Dict#{String, Any} # bindings for native grid \n
  dimension_dict::Dict\n
  model::String\n
  frequency::String # Day, month, years\n
  experiment::String # Historical, RCP4.5, RCP8.5, etc.\n
  run::String\n
  project::String # CORDEX, CMIP5, etc.\n
  institute::String # UQAM, DMI, etc.\n
  filename::String # Path of the original file\n
  dataunits::String # Celsius, kelvin, etc.\n
  latunits::String # latitude coordinate unit\n
  lonunits::String # longitude coordinate unit\n
  variable::String # Type of variable (i.e. can be the same as "typeofvar", but it is changed when calculating indices)\n
  typeofvar::String # Variable type (e.g. tasmax, tasmin, pr)\n
  typeofcal::String # Calendar type\n
  timeattrib::Dict # Time attributes (e.g. days since ... )\n
  varattribs::Dict # Variable attributes dictionary\n
  globalattribs::Dict # Global attributes dictionary\n
end\n
"""
function ClimGrid(data; longrid=[], latgrid=[], msk=[], grid_mapping=Dict(), dimension_dict=Dict(), timeattrib=Dict(), model="NA", frequency="NA", experiment="NA", run="NA", project="NA", institute="NA", filename="NA", dataunits="NA", latunits="NA", lonunits="NA", variable="NA", typeofvar="NA", typeofcal="NA", varattribs=Dict(), globalattribs=Dict())

    if isempty(dimension_dict)
        dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
    end

    if isempty(msk)
        msk = Array{Float64}(ones((size(data, 1), size(data, 2))))
    end

    ClimGrid(data, longrid, latgrid, msk, grid_mapping, dimension_dict, timeattrib, model, frequency, experiment, run, project, institute, filename, dataunits, latunits, lonunits, variable, typeofvar, typeofcal, varattribs, globalattribs)
end

"""
    TransferFunction(itp::Array, method::String, detrend::Bool)

Transfer function used during quantile-quantile mapping bias correction.
"""
struct TransferFunction
    itp::Array
    method::String
    detrend::Bool
end

# Included files
include("functions.jl")
include("indices.jl")
include("indicators.jl")
include("extract.jl")
include("interface.jl")
include("mapping.jl")
include("biascorrect.jl")
include("export.jl")
include("time.jl")
include("spatial.jl")
include("analysis.jl")

# Exported functions
export ClimGrid
# export uconvert
export inpoly, inpolygrid, meshgrid, inpolyvec, ndgrid
export findmax, findmin
export frostdays, summerdays, icingdays, tropicalnights
export daysabove10 #, daysbelow0, degdaysabove, degdaysbelow
export customthresover, customthresunder, annualmax, annualmin
export annualmean, annualsum, prcp1
export drought_dc
export ensemble_mean, ensemble_std, ensemble_max, ensemble_min
export load, load2D
export regrid, applymask, TransferFunction
export shapefile_coords, shapefile_coords_poly
export periodmean, resample, temporalsubset, spatialsubset
export qqmap, qqmaptf
export permute_west_east
export getdim_lat, getdim_lon, isdefined, extractpoly
export polyfit, polyval
export @isdefined
export plot, mapclimgrid, hist
export merge, vaporpressure, approx_surfacepressure
export wbgt, diurnaltemperature, meantemperature
export minimum, maximum, std, var, mean
export get_timevec, daymean, daysum
export monthmean, monthsum, temporalmean
export yearmonthdayhour
export write

end #module

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Conduit"
version = v"0.8.0"
sources = [
    ArchiveSource("https://github.com/LLNL/conduit/releases/download/v$(version)/conduit-v$(version)-src-with-blt.tar.gz",
		  "0607dcf9ced44f95e0b9549f5bbf7a332afd84597c52e293d7ca8d83117b5119")
]

# Note: we need to build the (optional) webserver by default until v0.8.1 is released
script = raw"""
cd ${WORKSPACE}/srcdir/conduit*
rm -rf build && mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      -DBUILD_SHARED_LIBS=ON \
      -DENABLE_TESTS=OFF \
      -DENABLE_EXAMPLES=OFF \
      -DENABLE_UTILS=OFF \
      -DENABLE_DOCS=OFF \
      -DENABLE_RELAY_WEBSERVER=ON \
      -DENABLE_COVERAGE=OFF \
      -DENABLE_PYTHON=OFF \
      -DENABLE_FORTRAN=OFF \
      -DENABLE_MPI=OFF \
      ../src
make -j${nproc}
make install
"""

platforms = expand_cxxstring_abis(supported_platforms())

products = [
    LibraryProduct("libconduit", :libconduit),
    LibraryProduct("libconduit_blueprint", :libconduit_blueprint),
    LibraryProduct("libconduit_relay", :libconduit_relay),
]

dependencies = Dependency[
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
	       julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
abstract type AbstractLineProperties end
struct LineProperties{T<:AF} <: AbstractLineProperties
    species::String
    mass::T             # amu
    depth::T
    λrest::T            # angstroms
    geff::T
    height::T           # km
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
struct Reference
    m::Module
    identifier::Union{Nothing, Symbol}
    fullname::String
    kind::Symbol
    ispublic::Bool
end

function Reference(m::Module, identifier::Symbol)
    t = referencetype(m, identifier)
    # Find module where reference was defined
    if t in (:function, :struct, :type, :pstruct)
        m = parentmodule(getfield(m, identifier))
    end

    return Reference(
        m,
        identifier,
        join((string(m), string(identifier)), '.'),
        t,
        Base.isexported(m, identifier),
    )
end

function referencetype(m::Module, identifier::Symbol)
    T = getfield(m, identifier)
    return referencetype(T)
end

function referencetype(x::DataType)
    isconcretetype(x) ? :struct : :type
end

function referencetype(x::UnionAll)
    isconcretetype(x) ? :pstruct : :ptype
end

function referencetype(::T) where T<:Function
    return :function
end

referencetype(x) = :const

referencetype(m::Module, ::Nothing = nothing) = :module

"""
    resolveidentifier(name, modules = ()) -> Union{Nothing, Reference}

Resolve an identifier `name` to a `Reference`. Return `nothing` if the identifier
can't be resolved in `Main` or `modules`.

```julia
resolveidentifier("sum", ("Base",)) == Reference(Base, :sum)
resolveidentifier("Base.sum") == Reference(Base, :sum)
```
"""
function resolveidentifier(name, modules = ())
    modulename, bindingname = splitidentifier(name)
    bindingsymbol = Symbol(bindingname)
    if modulename == ""
        for m in modules
            if isdefined(m, bindingsymbol)
                return Reference(m, bindingsymbol)

            end
        end
        return nothing
    else
        for m in modules
            # TODO: make generic to any level of submodule
            if isdefined(m, Symbol(modulename))
                submodule = getfield(m, Symbol(modulename))
                if isdefined(submodule, bindingsymbol)
                    return Reference(submodule, bindingsymbol)
                end
            end
        end
        m = getmodule(modulename)
        return isnothing(m) ? nothing : Reference(m, bindingsymbol)
    end
end

function populatereferences!(references, doc::XNode, linkfn = nothing, modules = ())
    sel = SelectTag(:a) & Pollen.SelectAttrEq(:href, "#")
    return cata(doc, sel) do x
        refname = gettext(x)
        refname == "" && return x
        ref = resolveidentifier(refname, modules)
        if isnothing(ref)
            @info "Could not resolve reference $refname in modules $modules."
            return x
        else
            references[ref.fullname] = ref
            return XNode(
                :a,
                Dict(:href => "/" * linkfn(ref.fullname)),
                [XNode(:code, [XLeaf(reflinkname(ref, modules))])],
            )
        end
    end
end

function reflinkname(ref::Reference, modules = ())
    if ref.ispublic && ref.m in modules
        return String(ref.identifier)
    else
        return ref.fullname
    end
end


text:
module CBinding
	import Libdl
	using Todo: @todo_str

	export Clongdouble, Caggregate, Cstruct, Cunion, Carray, Cenum, Clibrary, Cglobal, Cglobalconst, Cfunction, Cconvention, Calignment, Cconst, Caccessor
	export STDCALL, CDECL, FASTCALL, THISCALL
	export @ctypedef, @cstruct, @cunion, @carray, @calign, @cenum, @cextern, @cbindings
	export propertytypes

	# provide a temporary placeholder for 128-bit floating point primitive
	primitive type Clongdouble <: AbstractFloat sizeof(Cdouble)*2*8 end

	abstract type Cstruct end
	abstract type Cunion end
	const Caggregate = Union{Cstruct, Cunion}

	abstract type Cenum <: Integer end

	const Copaques = Union{Caggregate, Cenum}

	# alignment strategies
	struct Calignment{SymT}
	end

	const ALIGN_NATIVE = Calignment{:native}
	const ALIGN_PACKED = Calignment{:packed}

	# calling conventions
	struct Cconvention{SymT}
	end

	const STDCALL  = Cconvention{:stdcall}
	const CDECL    = Cconvention{:cdecl}
	const FASTCALL = Cconvention{:fastcall}
	const THISCALL = Cconvention{:thiscall}

	include("ctypespec.jl")
	include("clibrary.jl")
	include("cbindings.jl")
	include("cenum.jl")
	include("carray.jl")
	include("cconst.jl")
	include("caggregate.jl")
	include("cglobal.jl")
	include("cfunction.jl")
	include("cetc.jl")
	include("ctypelayout.jl")
	include("caccessor.jl")
end


text:
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("NeidSolarScripts")   end
 using Pkg
 Pkg.activate(".")

verbose = true
 make_plots = true
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using EchelleInstruments, EchelleInstruments.NEID
 using EchelleCCFs
 using RvSpectML
 if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, StatsBase, Statistics, Dates
 using JLD2, FileIO
 using NaNMath

target_subdir = "good_days"   # USER: Replace with directory of your choice
  fits_target_str = "Sun"
  output_dir = "output/"
  outputs = Dict{String,Any}()
  paths_to_search_for_param = [pwd(),"examples",joinpath(pkgdir(RvSpectMLBase),"..","RvSpectML","examples"), "/gpfs/scratch/jpn23/"]
  # NOTE: make_manifest does not update its paths_to_search when default_paths_to_search is defined here, so if you change the line above, you must also include "paths_to_search=default_paths_to_search" in the make_manifest() function call below
  pipeline_plan = PipelinePlan()
  dont_make_plot!(pipeline_plan, :movie)

reset_all_needs!(pipeline_plan)
#if need_to(pipeline_plan,:read_spectra)
if verbose println("# Finding what data files are avaliable.")  end
if isfile("manifest.csv")
    df_files  = CSV.read("manifest.csv", DataFrame)
    @assert size(df_files,1) >= 1
    @assert hasproperty(df_files,:Filename)
    @assert hasproperty(df_files,:target)
    @assert hasproperty(df_files,:bjd)
    @assert hasproperty(df_files,:ssbz)
    @assert hasproperty(df_files,:exptime)
else
    eval(read_data_paths(paths_to_search=paths_to_search_for_param))
    @assert isdefined(Main,:neid_data_path)
    df_files = make_manifest(neid_data_path, target_subdir, NEID )
    CSV.write("manifest.csv", df_files)
end

if verbose println("# Reading in customized parameters from param.jl.")  end
   if !@isdefined(idx_day_to_use)
       idx_day_to_use = 1
   end
   eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param))
   if match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1] ==  match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1]
      date_str = match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]
    else
      date_str = string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]) * "-" * string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1])
   end
   outputs["df_files_use"] = df_files_use

   outputs_filename = joinpath(output_dir,"solar_" * date_str * ".jld2")
   if isfile(outputs_filename) && false
     times_already_processed = load(outputs_filename, "times")
     files_in_day_to_process = size(df_files_solar_by_day.data[idx_day_to_use],1)
      if files_in_day_to_process == length(times_already_processed)
         println("# Already processed all ", length(times_already_processed), " files for ", date_str)
         exit()
      end
   end

alt_sun = map(i->calc_solar_alt(df_files_use[i,:bjd]),1:size(df_files_use,1))
obs_idx_max_alt_sun = argmax(alt_sun)
time_max_alt_sun = df_files_use[obs_idx_max_alt_sun,:bjd]
df_files_use.Δt = df_files_use[!,:bjd].-time_max_alt_sun
max_Δt = 0.5/24
df_files_use = df_files_use |> @filter( abs(_.Δt) <= max_Δt ) |> DataFrame

#for idx_day_to_use in 1:size(df_files_solar_by_day,1)
   #df_files_use = df_files_solar_by_day[idx_day_to_use,:data] |> @orderby(_.bjd) |> @take(max_spectra_to_use) |> DataFrame

if verbose println("# Reading in ", size(df_files_use,1), " FITS files for ", date_str, ".")  end
    @time all_spectra = map(row->NEID.read_data(row), eachrow(df_files_use))
    GC.gc()
    dont_need_to!(pipeline_plan,:read_spectra)

good_orders = 56:111 # Everything with λ and flux
 # good_orders = orders_to_use_default(NEID2D())
 #order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=max_orders, recalc=true )
 #order_list_timeseries = extract_orders(all_spectra,pipeline_plan,  orders_to_use=good_orders,  remove_bad_chunks=false, recalc=true )

df_orders_pixels = NEID.make_good_orders_pixels_df(all_spectra ; orders=good_orders)
chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_from_orders_pixels_df(all_spectra,NEID2D(), df_orders_pixels)

need_to!(pipeline_plan, :template)
 if need_to(pipeline_plan, :template)  # Compute order CCF's & measure RVs
    if verbose println("# Making template spectra.")  end
    #@assert !need_to(pipeline_plan,:extract_orders)
    map(i->chunk_list_timeseries.metadata[i][:rv_est] = 0.0, 1:length(chunk_list_timeseries) )
    GC.gc()   # run garbage collector for deallocated memory
    #map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_ccf[i]-mean(rvs_ccf), 1:length(order_list_timeseries) )
    # Smothing is broken with GP interpolation.  Need to fix.  In mean time, here's a Sinc interpolation workaround
    @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(chunk_list_timeseries, smooth_factor=2.0)
    #@time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries, alg=:Sinc)
    if save_data(pipeline_plan, :template)
       using JLD2, FileIO
       save(joinpath(output_dir,target_subdir * "_matrix.jld2"), Dict("λ"=>spectral_orders_matrix.λ,"spectra"=>spectral_orders_matrix.flux,"var_spectra"=>spectral_orders_matrix.var,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
    end
    dont_need_to!(pipeline_plan, :template);
 end

#=
chunkid = 21
 idx = spectral_orders_matrix.chunk_map[chunkid]
 plt = plot(spectral_orders_matrix.λ[idx],f_mean[idx],markersize=1.0,label="Template")
 plot!(plt,chunk_list_timeseries[1][chunkid].λ,chunk_list_timeseries[1][chunkid].flux, markersize=1.2,label="Obs  1")
 plot!(plt,chunk_list_timeseries[10][chunkid].λ,chunk_list_timeseries[10][chunkid].flux, markersize=1.2,label="Obs 10")
 plot!(plt,chunk_list_timeseries[20][chunkid].λ,chunk_list_timeseries[20][chunkid].flux, markersize=1.2,label="Obs 20")
 #plot!(plt,spectral_orders_matrix.λ[idx],deriv[idx]./mean(abs.(deriv[idx]))/2,markersize=1.1,label="Deriv")
 #plot!(plt,spectral_orders_matrix.λ[idx],deriv2[idx]./mean(abs.(deriv2[idx]))/2,markersize=1.1,label="Deriv 2")
 xlabel!("λ (Å)")
 ylabel!("f(λ)")
 title!("Template spectrum for chunk " * string(chunkid) )
 xlims!(spectral_orders_matrix.λ[idx[floor(Int,1+0.2*length(idx))]],spectral_orders_matrix.λ[idx[floor(Int,1+0.35*length(idx))]])

chunkid = 21
 idx = spectral_orders_matrix.chunk_map[chunkid]
 plt = plot(spectral_orders_matrix.λ[idx],(f_mean[idx].-1.0)./mean(abs.((f_mean[idx].-1.0))),markersize=1.0,label="Template")
 plot!(plt,spectral_orders_matrix.λ[idx],deriv[idx]./mean(abs.(deriv[idx])),markersize=1.1,label="Deriv")
 plot!(plt,spectral_orders_matrix.λ[idx],deriv2[idx]./mean(abs.(deriv2[idx])),markersize=1.1,label="Deriv 2")
 xlabel!("λ (Å)")
 ylabel!("f(λ), f'(λ), f''(λ), all standardized")
 title!("Template spectrum for chunk " * string(chunkid) )
 xlims!(spectral_orders_matrix.λ[idx[floor(Int,1+0.2*length(idx))]],spectral_orders_matrix.λ[idx[floor(Int,1+0.35*length(idx))]])
=#

line_width_50_default = 7.9e3
  lsf_width = 3.0e3
  max_mask_scale_factor = 4.0
  max_bc = RvSpectMLBase.max_bc
  #max_bc = RvSpectMLBase.max_bc_earth_rotation
  #if isfile(joinpath("data","solar_line_list_espresso.csv"))
   # line_list_espresso = CSV.read(joinpath("data","solar_line_list_espresso.csv"), DataFrame)
   # dont_need_to!(pipeline_plan,:clean_line_list_tellurics)
  #else
     #max_orders = min_order(NEID2D()):max_order(NEID2D())
     #good_orders = orders_to_use_default(NEID2D())
     #order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=max_orders, recalc=true )
     #order_list_timeseries = extract_orders(all_spectra,pipeline_plan,  orders_to_use=good_orders,  remove_bad_chunks=false, recalc=true )

     line_list_filename = "G2.espresso.mas"
     linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",line_list_filename)

     line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity,
        Δv_to_avoid_tellurics = 2*max_bc+5*line_width_50_default+max_mask_scale_factor*default_ccf_mask_v_width(NEID1D()), orders_to_use=good_orders, recalc=true, verbose=true)

     #CSV.write(joinpath("data","solar_line_list_espresso.csv"), line_list_espresso)
  #end
  #outputs["line_list_espresso"] = line_list_espresso
  line_width_50 = line_width_50_default

alt_sun = map(i->calc_solar_alt(all_spectra[i].metadata[:bjd]),1:length(all_spectra))
  obs_idx_max_alt_sun = argmax(alt_sun)
  time_max_alt_sun = chunk_list_timeseries.times[obs_idx_max_alt_sun]
  order_snrs = [RvSpectMLBase.calc_snr(chunk_list_timeseries[obsid].data[c]) for c in 1:length(chunk_list_timeseries[1].data),
                     obsid in 1:length(chunk_list_timeseries) ]
  order_weights =  (order_snrs[:,obs_idx_max_alt_sun]./order_snrs).^2
  for obsid in 1:length(chunk_list_timeseries)
     for ch in 1:length(chunk_list_timeseries[1].data)
         chunk_list_timeseries[obsid].data[ch].flux .*= order_weights[ch,obsid] #./mean_order_mean_flux[ch]
         chunk_list_timeseries[obsid].data[ch].var .*= (order_weights[ch,obsid]).^2 # ./mean_order_mean_flux[ch]).^2
     end
  end

msf = 0.62; fwtf = 1.5  # one pixel because tophat not using NEID's default v width
 #msf = 1.24; fwtf = 1.5   # two pixels
  println("# mask_scale_factor = ", msf, ", frac_width_to_fit =  ", fwtf, ". ")
  dont_need_to!(pipeline_plan,:extract_orders)
  ((ccfs_espresso, ccf_vars_espresso), v_grid) = ccf_total(chunk_list_timeseries, line_list_espresso, pipeline_plan,
    #mask_type=:gaussian,
    mask_scale_factor=msf, range_no_mask_change=5*line_width_50, ccf_mid_velocity=ccf_mid_velocity, v_step=155,
    v_max=max(5*line_width_50,2*max_bc), allow_nans=true, calc_ccf_var=true, recalc=true)

fwtf = 1.5; println("# mask_scale_factor = ", msf, ", frac_width_to_fit =  ", fwtf, ". ")
  alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
  rvs_ccf_espresso = RvSpectML.calc_rvs_from_ccf_total(ccfs_espresso, ccf_vars_espresso, pipeline_plan, v_grid=v_grid, times = chunk_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
  σ_rvs_espresso = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))

need_to!(pipeline_plan,:fit_lines)
 if need_to(pipeline_plan,:fit_lines)
    if verbose println("# Performing fresh search for lines in template spectra.")  end
    cl = ChunkList(map(grid->ChunkOfSpectrum(spectral_orders_matrix.λ,f_mean, var_mean, grid), spectral_orders_matrix.chunk_map), ones(Int64,length(spectral_orders_matrix.chunk_map)))
    # We're done with the spectral_orders_matrix, so we can release the memory now
    #spectral_orders_matrix = nothing
    #GC.gc()
    need_to!(pipeline_plan,:template)
    #lines_in_template_logy = RvSpectML.LineFinder.find_lines_in_chunklist(cl, plan=RvSpectML.LineFinder.LineFinderPlan(line_width=line_width_50,min_deriv2=0.5, use_logλ=true, use_logflux=true), verbose=true)  # TODO: Automate threshold for finding a line
    #line_width_50_default = 7.9e3
    #line_width_50 = line_width_50_default
    line_finder_plan = RvSpectML.LineFinder.LineFinderPlan(line_width=line_width_50,min_deriv2=0.5, use_logλ=true, use_logflux=false)
    lines_in_template = RvSpectML.LineFinder.find_lines_in_chunklist(cl, plan=line_finder_plan, verbose=true)  # TODO: Automate threshold for finding a line

    if verbose println("# Finding above lines in all spectra.")  end
    rvs_ccf = rvs_ccf_espresso
    @time fits_to_lines = RvSpectML.LineFinder.fit_all_lines_in_chunklist_timeseries(chunk_list_timeseries, lines_in_template, plan=line_finder_plan, rvs=rvs_ccf )

    #if save_data(pipeline_plan,:fit_lines)
       using CSV
       CSV.write(joinpath(output_dir, date_str * "_linefinder_lines.csv"), lines_in_template )
       CSV.write(joinpath(output_dir, date_str * "_linefinder_line_fits.csv"), fits_to_lines )
       #CSV.write(joinpath(output_dir,target_subdir * "_linefinder_line_fits_clean.csv"), lines_to_try )
    #end
    dont_need_to!(pipeline_plan,:fit_lines);
 end

function calc_lin_fit_coeff(times::AA1, y::AA2)  where { T1<:Real, T2<:Real, AA1<:AbstractArray{T1,1}, AA2<:AbstractArray{T2,1} }
    x = [ ones(length(times)) times ]
    linfit = (x'*x)\x'*y
 end

function calc_quad_fit_coeff(times::AA1, y::AA2)  where { T1<:Real, T2<:Real, AA1<:AbstractArray{T1,1}, AA2<:AbstractArray{T2,1} }
   x = [ ones(length(times)) times  times.^2 ]
   quadfit = (x'*x)\x'*y
end

function calc_ddt_coeff(times::AA1, y::AA2)  where { T1<:Real, T2<:Real, AA1<:AbstractArray{T1,1}, AA2<:AbstractArray{T2,1} }
   linfit = calc_lin_fit_coeff(times,y)
   slope = linfit[2]
end

using Printf

function select_line_fits_with_good_depth_width_slope(line_fits_df::DataFrame, quantile_threshold::Real; verbose::Bool = false, write_csv::Bool = false )
    conv_frac = line_fits_df |> @groupby(_.line_id) |> @map( { line_id=first(_.line_id), frac_converged=mean(_.fit_converged) } ) |> DataFrame

    fit_distrib = line_fits_df |> @filter(_.fit_converged == 1) |> @groupby(_.line_id) |>
             @map( { median_a=median(_.fit_a), median_b=median(_.fit_b), median_depth=median(_.fit_depth), median_σ²=median(_.fit_σ²), median_λc=median(_.fit_λc),
                    std_a=std(_.fit_a), std_b=std(_.fit_b), std_depth=std(_.fit_depth), std_σ²=std(_.fit_σ²), std_λc=std(_.fit_λc),
#                    Δt = df_files_use.Δt[_.obs_idx],
#                    ddt_depth=calc_ddt_coeff(df_files_use.Δt,_.fit_depth),
#                    ddt_σ²=calc_ddt_coeff(Δt,_.fit_σ²),
#                    ddt_λc=calc_ddt_coeff(Δt,_.fit_λc),
                    line_id=first(_.line_id) }) |>
             @join(conv_frac, _.line_id, _.line_id, {line_id=_.line_id, median_a=_.median_a, median_b=_.median_b, median_depth=_.median_depth, median_σ²=_.median_σ², median_λc=_.median_λc,
                    std_a=_.std_a, std_b=_.std_b, std_depth=_.std_depth, std_σ²=_.std_σ², std_λc=_.std_λc,
#                    ddt_depth=_.ddt_depth, ddt_σ²=_.ddt_σ², ddt_λc = _.ddt_λc,
                    frac_converged=__.frac_converged } ) |>
             @filter(_.frac_converged > 0.5 ) |> DataFrame

    std_depth_treshold = quantile(fit_distrib.std_depth,quantile_threshold)
    #median_σ_width_treshold = quantile(sqrt.(fit_distrib.median_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,1-quantile_threshold)
    median_σ_width_treshold_lo = 2000
    #median_σ_width_treshold_lo = 100
    #median_σ_width_treshold = line_width_50/3
    std_σ_width_treshold = quantile(sqrt.(fit_distrib.std_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,quantile_threshold)
    #std_depth_treshold = quantile(fit_distrib.std_depth,quantile_threshold)
    #median_σ_width_treshold_lo = quantile(sqrt.(fit_distrib.median_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,(1-quantile_threshold)/2)
    median_σ_width_treshold_hi = quantile(sqrt.(fit_distrib.median_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,1-(1-quantile_threshold)/2)
    std_b_treshold = quantile(fit_distrib.std_b,quantile_threshold)
    std_a_treshold = quantile(fit_distrib.std_a,quantile_threshold)
    good_lines_alt = fit_distrib |>
       @filter( 0.05 <= _.median_depth <= 1.0 ) |>
       #@filter( sqrt.(_.median_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps >= median_σ_width_treshold_lo ) |>
       @filter( sqrt.(_.std_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps <= std_σ_width_treshold ) |>
       #@filter(  median_σ_width_treshold_lo <= sqrt.(_.median_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps <= median_σ_width_treshold_hi ) |>
       @filter( _.std_depth <= std_depth_treshold ) |>
       @filter( _.std_b < std_b_treshold) |>
       DataFrame
    if verbose
       println("# Found ", size(good_lines_alt,1), " good lines (std_depth_width_slope), rejected ", size(fit_distrib,1)-size(good_lines_alt,1), " lines.")
    end
    if write_csv
       val_str = Printf.@sprintf("%1.2f",quantile_threshold)
       CSV.write(joinpath(output_dir,date_str * "_good_lines_fit_quant=" * val_str * ".csv"), good_lines_alt )
    end
    return good_lines_alt
 end

function calc_rvs_from_line_by_line_fits_alt(fits_to_lines::DataFrame, threshold_alt::Real)
    good_lines_alt = select_line_fits_with_good_depth_width_slope(fits_to_lines,threshold_alt,write_csv=true)
    df2 = fits_to_lines |> @join(good_lines_alt, _.line_id, _.line_id, {obs_id=_.obs_idx, line_id=_.line_id, Δv=(_.fit_λc.-__.median_λc).*RvSpectML.speed_of_light_mps./__.median_λc, weight=(min(__.median_λc/__.std_λc, RvSpectML.speed_of_light_mps/750))^2
    #= (__.median_λc/(RvSpectML.speed_of_light_mps*__.std_λc))^2 =# }) |>
             @groupby(_.obs_id) |> @map({obs_id=first(_.obs_id), Δv= sum(_.Δv.*_.weight)./sum(_.weight) }) |> DataFrame
    rms_alt = std(df2.Δv)
 end

good_lines_alt = select_line_fits_with_good_depth_width_slope(fits_to_lines,0.8) #0.8)
df2 = fits_to_lines |> @join(good_lines_alt, _.line_id, _.line_id, {obs_id=_.obs_idx, line_id=_.line_id, Δv=(_.fit_λc.-__.median_λc).*RvSpectML.speed_of_light_mps./__.median_λc, weight=(min(__.median_λc/__.std_λc, RvSpectML.speed_of_light_mps/750))^2 }) |>
          @groupby(_.obs_id) |> @map({obs_id=first(_.obs_id), Δv= sum(_.Δv.*_.weight)./sum(_.weight) }) |> DataFrame

using Plots
#thresholds = [ 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95]
thresholds = [ 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0 ]
 #thresholds = [0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0 ]
 #thresholds = [0.75, 0.8, 0.9, 0.95, 1.0 ]
  plt = plot()
  #out_vec_calc_rvs_from_line_by_line_fits_stdv = map(t->calc_rvs_from_line_by_line_fits_stdv(fits_to_lines,t),thresholds)
  #scatter!(plt,thresholds,out_vec_calc_rvs_from_line_by_line_fits_stdv,color=1,label="σ_λ")
  #plot!(plt,thresholds,out_vec_calc_rvs_from_line_by_line_fits_stdv,color=1,label=:none)
  out_vec_calc_rvs_from_line_by_line_fits_alt = map(t->calc_rvs_from_line_by_line_fits_alt(fits_to_lines,t),thresholds)
  scatter!(plt,thresholds,out_vec_calc_rvs_from_line_by_line_fits_alt,color=2,label="σ_other")
  plot!(plt,thresholds,out_vec_calc_rvs_from_line_by_line_fits_alt,color=2,label=:none)
  xlabel!("Threshold for accepting lines")
  ylabel!("RMS v (m/s)")
  savefig("rms_rv_vs_line_accept_threshold" * date_str * ".png")
  #display(plt)

#hpf_mask_data = readdlm(joinpath(pwd(),"data","HPF_TelMask_interpolated.txt"))


text:
rtol = 1e-2

# sunrise at equator
date = Dates.DateTime(2020, 2, 20, 6, 11, 0)
lon, lat = [FT(0.0), FT(0.0)]
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test sza ≈ π/2 rtol=rtol

# solar noon at equator
date = Dates.DateTime(2020, 2, 20, 12, 14, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test azi ≈ 3π/2 rtol=rtol

# sunset at equator
date = Dates.DateTime(2020, 2, 20, 18, 17, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test sza ≈ π/2 rtol=rtol

# solar noon at equator, eot correction = false
date = Dates.DateTime(2020, 2, 20, 12, 0, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set; eot_correction=false)
@test azi ≈ 3π/2 rtol=rtol

# sunset at equator, eot correction = false
date = Dates.DateTime(2020, 2, 20, 18, 0, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set; eot_correction=false)
@test sza ≈ π/2 rtol=rtol

## Test Polar Night
# polar night NH 1
date = Dates.DateTime(2020, 12, 20, 11, 0, 0)
lon, lat = [FT(0.0), FT(80.0)]
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test sza > π/2

# polar night NH 2
date = Dates.DateTime(2020, 12, 20, 23, 0, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test sza > π/2

# polar night SH 1
date = Dates.DateTime(2020, 6, 20, 11, 0, 0)
lon, lat = [FT(0.0), FT(-80.0)]
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test sza > π/2

# polar night SH 2
date = Dates.DateTime(2020, 6, 20, 23, 0, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test sza > π/2

## Test Distance
date = Dates.DateTime(2000, 3, 22, 0, 0, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
@test d ≈ orbit_semimaj(param_set) rtol=rtol


text:
using BinDeps
import BinDeps: PackageManager, can_use, package_available, available_version,
    libdir, generate_steps, LibraryDependency, provider, provides, pkg_name
type Mingw64 <: PackageManager
  package :: AbstractString
  Mingw64(package :: AbstractString) = new("mingw-w64-x86_64-" * package)
end

can_use(::Type{Mingw64}) = OS_NAME == :Windows
package_available(m :: Mingw64) = can_use(Mingw64) && exists(m.package)
available_version(m :: Mingw64) = get(info(m.package), "Version")
pkg_name(m :: Mingw64) = m.package
libdir(m :: Mingw64, dep) = joinpath(Pkg.dir("MSys2"), "deps", "msys64", "mingw64", "bin")

function generate_steps(dep :: LibraryDependency, m :: Mingw64, opts)
  if get(opts,:force_rebuild,false)
    error("Will not force yum to rebuild dependency \"$(dep.name)\".\n"*
        "Please make any necessary adjustments manually (This might just be a version upgrade)")
  end
  () -> install(m.package)
end


text:
"""
REINFORCE
---------

Williams, 1992

This implementation follows the pseudo code given in
Sutton Barto, "RL an Intro", Chapter 13 Monte Carlo Policy Gradient.

The policy parameters θ as well as the feature vector ϕ are
column vectors. Here we use a states x actions matrix for simpler access.
Also, the feature vector ϕ is a one-hot vector in our case. Therefore we
do not need to keep a vector but access the entries by the
state/action pair [s_t, a_t]

See learner/REINFORCElearner.jl for the core implementation.

Note
----
get_all_puzzles(), get_all_episodes_in_puzzle(), get_data_by_puzzle_episode()
are in rl_data_toolsforalgos.jl
"""
function getLL_REINFORCE(data::Array{Int,2}, alpha::Float64,
                        gamma::Float64, temp::Float64; do_debug_print = false)
    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_REINFORCE supports only one subject")
    end
    if nr_actions!=2; error("REINFORCE supports two actions only") end
    all_puzzles = get_all_puzzles(data)

    LL_unbiased = 0.
    action_counter = 0
    probs_unbiased = zeros(Float64, nr_samples)
    Δpolicy_record = zeros(Float64, nr_samples)

    PolicyParametersImbalance_record = zeros(Float64, nr_samples)
    PolicyImbalance_record = zeros(Float64, nr_samples)
    PolicyParametersImbalanceAllstates_record = Array{Array{Float64,1}, 1}(undef, nr_samples)
    [PolicyParametersImbalanceAllstates_record[i] = zeros(Float64, nr_states) for i in 1:nr_samples]
    PolicyParametersDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)
    PolicyDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)

    modulatingfactor = 1.
    learner = []
    for i_puzzle in all_puzzles
        # Puzzle = independent block of the experiment.
        # The code was made for general purposes, in this experiment we
        # have only ONE puzzle.
        learner = REINFORCE(nr_states = nr_states,
                        nr_actions = nr_actions,
                        α = alpha,
                        γ = gamma,
                        do_debug_print = do_debug_print)

        all_epi_in_puzzle = get_all_episodes_in_puzzle(data, i_puzzle)
        for i_episode in all_epi_in_puzzle
            episode_data = get_data_by_puzzle_episode(data, i_puzzle, i_episode)
            path_length = size(episode_data)[1]  # = T
            if do_debug_print
                println("i_puzzle:$i_puzzle | i_episode:$i_episode | path_length:$path_length")
            end
            reward = 0. # we have a single reward at goal state. If the episode was
            # not completed, this code still correctly works (updates are 0)
            for t in 1:path_length
                action_counter += 1
                state = episode_data[t, Int(c_state)]
                action = episode_data[t, Int(c_action)]
                reward = episode_data[t, Int(c_reward)]
                next_state = episode_data[t, Int(c_nextState)]
                done = reward>0 ? true : false

                p_actions = getSoftmaxProbs([learner.θ[state,1],learner.θ[state,2]], temp, [0., 0.])
                LL_unbiased += log(p_actions[action])
                probs_unbiased[action_counter] = copy(p_actions[action])

                PolicyParametersImbalance_record[action_counter] = abs(learner.θ[state, 1] - learner.θ[state, 2])
                PolicyImbalance_record[action_counter] = abs(p_actions[1] - p_actions[2])
                PolicyParametersImbalanceAllstates_record[action_counter] = abs.(learner.θ[:, 1] - learner.θ[:, 2])
                PolicyParametersDiffChosenVsUnchosen_record[action_counter] = learner.θ[state, action] - learner.θ[state, mod(action,2)+1]
                PolicyDiffChosenVsUnchosen_record[action_counter] = p_actions[action] - p_actions[mod(action,2)+1]

                θ_before = deepcopy(learner.θ)

                update!(learner, state, action, reward, p_actions[action], modulatingfactor)
                Δpolicy_record[action_counter] = sum(abs.(learner.θ .- θ_before))

                if do_debug_print
                    println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | p(s,a):$(round(p_action,digits=2))")
                end
            end # path_length (in single episode)
        end # episode
    end # puzzzle
  #---- Without bias
  LL = LL_unbiased
  trial_by_trial_p_sa = (probs_unbiased)

  algo = "REINFORCE"
  param_v = [alpha, gamma, temp]
  param_str = "alpha=$(round(alpha,digits=4)), gamma=$(round(gamma,digits=4)), temp=$(round(temp,digits=4))"

  results_dict = Dict{String,Any}("final_Theta" => learner.θ,
                                "DeltaPolicy" => Δpolicy_record,
                                "PolicyParametersImbalance" => PolicyParametersImbalance_record,
                                "PolicyImbalance" => PolicyImbalance_record,
                                "PolicyParametersImbalanceAllstates" => PolicyParametersImbalanceAllstates_record,
                                "PolicyParametersDiffChosenVsUnchosen" => PolicyParametersDiffChosenVsUnchosen_record,
                                "PolicyDiffChosenVsUnchosen" => PolicyDiffChosenVsUnchosen_record)

  signals = Step_by_step_signals(
    data_stats.all_subject_ids[1],
    [], # no RPE
    [], # no SPE
    trial_by_trial_p_sa,
    LL,
    data_stats,
    algo,
    param_v,
    param_str,
    99, #biased_action_idx, dummy value
    results_dict)
  return LL, signals

end


text:
import StaticArrays:
	SArray,
    MArray

################################################################################
# Static simplex where vertices are represented by some type of abstract vector.
################################################################################

struct Simplex{D, T} <: AbstractImmutableSimplex{D, T}
    vertices::Vector{AbstractVector{T}}
end

"""
    Simplex(pts::Vector{AbstractVector})

Construct a simplex from a vector of vectors.
"""
function Simplex(pts::Vector{Vector{T}}) where {T}
    if !(length(pts) == length(pts[1]) + 1)
        err = """ The input cannot be converted to a simplex.
        Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices.
        """
        throw(DomainError(pts, err))
    end
    dim = length(pts[1])
    Simplex{dim, T}([pts[i] for i = 1:length(pts)])
end

"""
    Simplex(pts::Vector{SVector{D, T}})

Construct a simplex from a vector of static vectors.
"""
function Simplex(pts::Vector{SArray{Tuple{D},T,1,D}}) where {D, T}
    if !(length(pts) == length(pts[1]) + 1)
        err = """ The input cannot be converted to a simplex.
        Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices.
        """
        throw(DomainError(pts, err))
    end
    dim = length(pts[1])
    Simplex{dim, T}([pts[i] for i = 1:length(pts)])
end

"""
    Simplex(pts::Vector{MVector{D, T}})

Construct a simplex from a vector of mutable static vectors.
"""
function Simplex(pts::Vector{MArray{Tuple{D},T,1,D}}) where {D, T}
    if !(length(pts) == length(pts[1]) + 1)
        err = """ The input cannot be converted to a simplex.
        Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices.
        """
        throw(DomainError(pts, err))
    end
    dim = length(pts[1])
    Simplex{dim, T}([pts[i] for i = 1:length(pts)])
end

"""
    Simplex(pts::AbstractArray{T, 2}) where T

Construct a simplex from an array of size `(dim + 1)-by-(dim)` or
`(dim)-by-(dim + 1)` (faster).
"""
function Simplex(pts::AbstractArray{T, 2}) where {T}
    s = size(pts)

    if (maximum(s) - minimum(s)) > 1
        err = """ The input cannot be converted to a simplex.
            size(pts) must be (dim, dim + 1) or (dim + 1, dim).
        """
        throw(DomainError(pts, err))
    end

    # Assuming rows of the input are vertices
    if s[1] > s[2]
        dim = s[2]
        return Simplex{dim, T}([pts[i, :] for i = 1:maximum(s)])
    end

    # Assuming columns of the input are vertices
    if s[2] > s[1]
        dim = s[1]
        return Simplex{dim, T}([pts[:, i] for i = 1:maximum(s)])
    end
end

export Simplex


text:
"""
OBJECTIVE
=========
This code anayzes the single-player game of snakes and ladders as a
discrete-time Markov chain with a single absorbing state, and calculates
[1] the expected number of turns to complete a game,
[2] the probability that a player ever moves to square j, given that she started
    in square j,
[3] the least number of moves (turns) to complete the game,
[4] the most likely path (or, sequence of board moves) to complete the game, and
[5] the probability that a player completes the game in k moves, k = 1, 2, ...

CHOICE OF BOARD
===============
[1] https://en.wikipedia.org/wiki/File:Cnl03.jpg
Input files: numberofsquares.tx, snakes.txt, ladders.txt

[2] https://www.etsy.com/listing/764625917/snakes-ladders-vintage-game-board-png
Input files: numberofsquares.txt, snakes_vintage.txt, ladders_vintage.txt

RESOURCES
=========
Online forums:
[1] https://math.stackexchange.com/questions/691494/expected-number-of-steps-between-states-in-a-markov-chain
[2] https://www.quora.com/What-is-the-expected-number-of-moves-required-to-finish-a-snakes-and-ladders-game
[3] https://en.wikipedia.org/wiki/Discrete_phase-type_distribution
[4] https://math.stackexchange.com/questions/1695191/markov-chain-snakes-and-ladders
[5] https://stackoverflow.com/questions/38894477/find-the-path-with-the-max-likelihood-between-two-vertices-in-markov-model
[6] https://mathematica.stackexchange.com/questions/19457/finding-the-likeliest-path-in-a-markov-process
[7] https://math.stackexchange.com/questions/3281928/what-is-the-distribution-of-time-to-absorption-for-an-absorbing-markov-chain

Lecture notes, videos, and book chapters:
[8] https://ocw.mit.edu/resources/res-6-012-introduction-to-probability-spring-2018/part-iii-random-processes/expected-time-to-absorption/
[9] http://www2.imm.dtu.dk/courses/02407/lectnotes/ftf.pdf
[10] Introduction to Probability Models, 11th Ed., Sheldon M. Ross, Section 4.6
[11] Introductionto Probability, 2nd Ed., Charles M. Grinstead and
     J. Laurie Snell, Chapter 11
     https://math.dartmouth.edu/~prob/prob/prob.pdf

AUTHOR
======
Kalyani Nagaraj
September 18, 2019

KEYWORDS
========
Snakes and ladders, Markov chain, discrete phase-type distribution, Dijkstra's
shortest path algorithm
"""

# 1. INCLUDE PACKAGES
# -------------------
using LinearAlgebra, Graphs
using DelimitedFiles

# 2. READ BOARD GAME DATA FROM FILE, CREATE A DICTIONARY
# ------------------------------------------------------
num_squares = readdlm("Data/numberofsquares.txt", Int64);
num_states = num_squares[1] + 1;
# First column of s  : locations of snake heads
# Second column of s : locations of snake tails
s = readdlm("Data/snakes_punitive.txt", Int64);
# First column of l  : location of bottom of ladder
# Second column of l : location of top of ladder
l = readdlm("Data/ladders_punitive.txt", Int64);
# Create a dictionary whose keys are path sources and whose corresponding
# values are path destination implied by the snakes and ladders
snakes  = [ [ s[i,1] s[i,2] ] for i in 1:size(s,1)]
ladders = [ [ l[i,1] l[i,2] ] for i in 1:size(l,1)]
dictJumps = Dict([la[1] => la[2] for la in ladders])
merge!(dictJumps, Dict([sn[1] => sn[2] for sn in snakes]))

# 3. GENERATE THE STATE TRANSITION MATRIX P USING BOARD GAME DATA
# ---------------------------------------------------------------
P = zeros(num_states, num_states)  # Initialize a 101 X 101 1-step transition
                                   # probability matrix
                                   # The game begins from square 0
for k = 1:(num_states-1)-5   # Loops from 1 to 95
   for j = k+1:k+6           # Fill in the transition probs. for a fair die
      P[k,j] = 1/6
   end
end
for k = (num_states-1)-4:num_states    # Loops from 96 to 101
                                       # Player must roll an exact number to
                                       # reach square 100
   for j = k+1:num_states
      P[k,j] = 1/6
   end
   P[k,k] = (k + 6 - num_states)/6  # If player overshoots square 100, she remains
                                    # in square k with probability (k+6-101)/6
end
# Now, adjust the transition matrix to account for jumps in position due snakes
# or ladders
for k in keys(dictJumps)
   P[:, dictJumps[k]+1] += P[:, k+1]
end
# Notice that P[i, j] = 0 if either square i or square j has a snake's head or the
# bottom of a ladder. So we can safely delete columns and rows corresponding to
# the start of "jumps", without affecting our analysis.
# This matrix truncation (naturally) reduces the size of the P,
# which in turn will make the matrix inversion operation go faster.
P = P[setdiff(collect(1:num_states), collect(keys(dictJumps)) .+ 1),
  setdiff(collect(1:num_states), collect(keys(dictJumps)) .+ 1)]

# 4. CALULATE THE EXPECTED NUMBER OF TURN TO COMPLETE THE GAME
# ------------------------------------------------------------
# P is a one-step transition matrix that corresponds to a discrete-time
# Markov chain with a single absorbing state.
# P = [T t;
#      0 1]
# Remove the last row and column to get the non-singular matrix T.
# The expected number of moves to complete the game (i.e., absorption times),
# when starting from squares 0 to 99 (minus the squares that were removed while
# truncating P), are given by E[Time] = inv(I - P) * vec(1).
T  = P[1:end-1, 1:end-1]
N  = size(T, 1)
IN = Matrix{Float64}(I, N, N)  # IN is the identity matrix of size N
M  = IN - T
S  = inv(M)  # Matrix of mean times s_ij the Markov chain is in state j given
             # that it starts in state i
A = S*ones(N) # Vector of mean absorption times a_i when starting in state i
println("Mean number of moves to complete the game = ", A[1], ".")

# 5. CALCULATE THE PROBABILITIES f_ij THAT THE MARKOV CHAIN EVER TRANSITIONS TO
# STATE j, GIVEN THAT ITS STARTS IN STATE i.
# -----------------------------------------------------------------------------
# f_ij = (s_ij - delta_ij) / s_jj,
# where delta_ij = 1 if i=j, 0 otherwise.
# Note that the matrix F is made up of probabilities only for the transient
# states.
F = (S .- IN) * diagm(0 => (1 ./ diag(S)))

# 6. FIND THE SMALLEST NUMBER OF MOVES TO COMPLETE THE GAME BY SOLVING A
# SHORTEST PATH PROBLEM (APPLY DIJKSTRA'S ALGORITHM)
# ----------------------------------------------------------------------
# 7. FIND THE MOST LIKELY PATH (AND THE LIKELIHOOD) TO COMPLETE THE GAME BY
# SOLVING ANOTHER SHORTEST PATH PROBLEM WITH EDGE WEIGHTS = -log(P[i, j])
# -------------------------------------------------------------------------
vnames = setdiff(collect(0:num_states-1), collect(keys(dictJumps))) # graph vertices
g = inclist(vnames, is_directed=true)
                            # Initialize a simple directed graph with vertex list
                            # `vnames` and no edges
                            # Note that vnames does not include the bottoms of
                            # ladders and the heads of snakes
AdjMat = Int.(sign.(P))     # Construct the adjacency matrix of graph g.
ne = sum(AdjMat)            # Count the number of edges from the adj. matrix.
dists = zeros(ne)           # Intialize each edge weight to 0.
k = 0
for i = 1:size(P,1)         # Loop to add an edge whenever trans. prob P_ij>0
   for j = 1:size(P,2)
      if P[i,j] > 0         # add edge from vnames[i] to vnames[j]
         add_edge!(g, vnames[i], vnames[j])
         global k = k + 1
         dists[k] = P[i,j] < 1 ? -log(P[i,j]) : 0.0
                                   # weight of the egde = - log(P_ij)
      end
   end
end
r1 = dijkstra_shortest_paths(g, 0) # Find smallest number of moves to go from
                                   # square 0 to squares 0, 1, 2, ..., 100
                                   # All edge weights equal 1
r2 = dijkstra_shortest_paths(g, dists, 0) # Find the most likely path from
                                          # square 0 to squares 0, 1, 2, ... 100
                                          # when edge weights = - log(P_ij)
# Enumerate the shortest paths from square 0 to all other vertices:
enpath_1 = enumerate_paths(vertices(g), r1.parent_indices)
enpath_2 = enumerate_paths(vertices(g), r2.parent_indices)
# When all edge weights equal one:
println("\nThe game can be completed in a minimum of ", r1.dists[end], " moves.")
println("One such set of moves is ", enpath_1[end], ".\n")
# When edge weights equal negative log of one-step transition probabilities:
println("\nThe most likely path from square 0 to square ", num_states-1, " is ", enpath_2[end], ".")
println("This path requires ", length(enpath_2[end])-1,
        " moves and is achieved with probility ", exp(-r2.dists[end]), ".")

# 8. COMPUTE THE PMF AND CDF OF RANDOM VARIABLE N, TIME TO COMPLETE THE GAME
# --------------------------------------------------------------------------
# Define the initial distribution tau
tau = zeros(size(T,1)); tau[1] = 1;
T0 = P[1:end-1, end]
pmf_phasetypedist = zeros(300);
cdf_phasetypedist = zeros(300);
for n=1:300
    pmf_phasetypedist[n] = transpose(tau) * T^(n-1) * T0;
    cdf_phasetypedist[n] = 1 - transpose(tau) * T^n * ones(size(T,1))
end
using Plots
cdfplot = Plots.plot(0:300, [0;cdf_phasetypedist], linetype=:steppre, title = "Discrete Phase-Type Distribution", label = "", lw = 1)
xlabel!("Number of moves to complete game")
ylabel!("Cummulative density function")
display(cdfplot)
savefig("Plots/cdf_punitive.pdf")
pmfplot = Plots.scatter(1:300, pmf_phasetypedist, title = "Discrete Phase-Type Distribution", label = "", lw = 1)
xlabel!("Number of moves to complete game")
ylabel!("Probability mass function")
display(pmfplot)
savefig("Plots/pmf_punitive.pdf")

max_stars_repo_name: kalyaninagaraj/Snakes-And-Ladders

max_issues_repo_name: kalyaninagaraj/snakesandladders

max_forks_repo_name: kalyaninagaraj/snakesandladders


text:
using Images, ImageView
img = load("sample-images/cats-3061372_640.jpg")
imshow(img[1:1, 1:1, :])


text:
@testset "simple workflow" begin
    env = StateTransformedEnv(CartPoleEnv{Float32}();state_mapping=deepcopy)
    policy = RandomPolicy(action_space(env))
    N_EPISODE = 10_000
    hook = TotalRewardPerEpisode()
    run(policy, env, StopAfterEpisode(N_EPISODE), hook)

    @test isapprox(sum(hook[]) / N_EPISODE, 21; atol=2)
end

@testset "multi agent" begin
    # https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/issues/393
    rps = RockPaperScissorsEnv() |> SequentialEnv
    ma_policy = MultiAgentManager(
        (
            NamedPolicy(p => RandomPolicy())
            for p in players(rps)
        )...
    )

    run(ma_policy, rps, StopAfterEpisode(10))
end


text:
# This is an example script demonstrating how PARDISO works on a small,
# sparse, real symmetric matrix. It computes the m solutions X to the
# collection of linear systems
#
#    A * X = B
#
# using the PARDISO solver, where A is a symmetric n x n matrix, B is an
# n x m matrix, and X is another n x m matrix.
using Pardiso
using SparseArrays
using Random
using Printf
using Test

# Script parameters.
# -----------------
verbose = false
n       = 4  # The number of equations.
m       = 3  # The number of right-hand sides.

A = sparse([ 1. 0 -2  3
             0  5  1  2
            -2  1  4 -7
             3  2 -7  5 ])

# Generate a random collection of right-hand sides.
B = rand(n,m)

# Initialize the PARDISO internal data structures.
# ps = PardisoSolver()
ps = MKLPardisoSolver()

if verbose
    set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
end

# If we want, we could just solve the system right now.
# Pardiso.jl will automatically detect the correct matrix type,
# solve the system and free the data
X1 = solve(ps, A, B)

# We also show how to do this in incremental steps.

# ps = PardisoSolver()
ps = MKLPardisoSolver()

# First set the matrix type to handle general real symmetric matrices
set_matrixtype!(ps, Pardiso.REAL_SYM_INDEF)

# Initialize the default settings with the current matrix type
pardisoinit(ps)

# Remember that we pass in a CSC matrix to Pardiso, so need
# to set the transpose iparm option.
fix_iparm!(ps, :N)

# Get the correct matrix to be sent into the pardiso function.
# :N for normal matrix, :T for transpose, :C for conjugate
A_pardiso = get_matrix(ps, A, :N)

# Analyze the matrix and compute a symbolic factorization.
set_phase!(ps, Pardiso.ANALYSIS)
set_perm!(ps, randperm(n))
pardiso(ps, A_pardiso, B)
@printf("The factors have %d nonzero entries.\n", get_iparm(ps, 18))

# Compute the numeric factorization.
set_phase!(ps, Pardiso.NUM_FACT)
pardiso(ps, A_pardiso, B)
@printf("The matrix has %d positive and %d negative eigenvalues.\n",
        get_iparm(ps, 22), get_iparm(ps, 23))

# Compute the solutions X using the symbolic factorization.
set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
X = similar(B) # Solution is stored in X
pardiso(ps, X, A_pardiso, B)
@printf("PARDISO performed %d iterative refinement steps.\n", get_iparm(ps, 7))

# Compute the residuals.
R = maximum(abs.(A*X - B))
@printf("The maximum residual for the solution X is %0.3g.\n", R)
@test R < 1e-10

# Free the PARDISO data structures.
set_phase!(ps, Pardiso.RELEASE_ALL)
pardiso(ps)


text:
struct RMP2{T} <: AbstractMPWavefunction# where T <: AbstractFloat
    CorrelationEnergy::T
    T2::_T2 where _T2 <: Fermi.AbstractTensor
end

function RMP2()
    refwfn = Fermi.HartreeFock.RHF()
    RMP2(refwfn)
end
function RMP2(ref::Fermi.HartreeFock.RHF)
    precision_selector = Dict{Any,Any}("single" => Float32,
                                       "double" => Float64)
    prec = precision_selector[Fermi.CurrentOptions["precision"]]
    RMP2{prec}(ref)
end

function RMP2{T}(ref::Fermi.HartreeFock.RHF) where T <: AbstractFloat
    # convert some options -> singletons
    alg_selector = Dict{Any,Any}("direct" => Direct(),
                                 "df" => DF(),
                                 "conv" => Conventional())
    alg = alg_selector[Fermi.CurrentOptions["mp2_type"]]
    RMP2{T}(ref,alg)
end

include("DF-RMP2.jl")
include("ConventionalRMP2.jl")
#"""
#This is the most generic constructor, so we'll put the docstring here
#"""
#function RMP2{T}(ref::Fermi.HartreeFock.RHF,alg::T2) where { T <: AbstractFloat,
#                                                                             T2 <: MP2Algorithm }
#end

# Concrete implementations of RMP2
#include("DirectRMP2.jl")


text:
using Instruments
using Test

# write your own tests here
@test 1 == 1
rm = ResourceManager()
insts = Vector{String}()
try
    push!(insts, Instruments.viFindRsrc(rm, "?*::INSTR")...)
catch LoadError
    viClose(rm)
end


text:
using Documenter, SigFigs

makedocs(;
    modules=[SigFigs],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/caseykneale/SigFigs.jl/blob/{commit}{path}#L{line}",
    sitename="SigFigs.jl",
    authors="Casey Kneale",
    assets=String[],
)

deploydocs(;
    repo="github.com/caseykneale/SigFigs.jl",
)


text:
module Immerse

using GtkUtilities, Colors, Reexport, Compat, REPL, Graphics
@reexport using Gadfly
import Gtk   # because both Gadfly and Gtk define draw
import Measures
using Compose, Cairo

export
    Figure,
    closefig,
    closeall,
    figure,
    gcf,
    scf,
    getproperty,
    setproperty!,
    hit,
    lasso_initialize

# Stuff for Gadfly/Compose
include("compose.jl")
# using .ImmerseCompose
include("display_gadfly.jl")
# using .DisplayGadfly

# Generic (?)
include("hit_test.jl")
# using .HitTest
include("select.jl")

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

function __init__()
    if !isprecompiling()
        Compose.link_cairo()
    end

    #set white background for default Theme
    t = Theme()
    t.panel_fill = colorant"white"
    Gadfly.push_theme(t)

    pushdisplay(_display)
end

end # module


text:
using LikelihoodProfiler

include("./cases_func.jl")

#include("../src/params_intervals.jl")
# using ParametersIdentification

res1 = params_intervals(
    [3., 4, 1.1, 10.],
    3,
    9.,
    f_3p_1im,
    logscale_all = true,
    method = :ONE_PASS
)

res2 = params_intervals(
    [3., 4, 1.1, 10.],
    1,
    9.,
    f_3p_1im,
    logscale_all = false,
    method = :D2D_PLE
)


text:
function GillFieldIntegrator(arg0::Field, arg1::RealFieldElement)
    return GillFieldIntegrator((Field, RealFieldElement), arg0, arg1)
end

function get_a(obj::GillFieldIntegrator)
    return jcall(obj, "getA", Vector{Vector{RealFieldElement}}, ())
end

function get_b(obj::GillFieldIntegrator)
    return jcall(obj, "getB", Vector{RealFieldElement}, ())
end

function get_c(obj::GillFieldIntegrator)
    return jcall(obj, "getC", Vector{RealFieldElement}, ())
end


text:
const ϵ = 1e-8

mutable struct TruncatedADAGrad
    eta::Float64
    tau::Float64
    n::Int

    iters::IdDict
    acc::IdDict
end

function TruncatedADAGrad(η = 0.1, τ = 1.0, n = 100)
    TruncatedADAGrad(η, τ, n, IdDict(), IdDict())
end

function apply!(o::TruncatedADAGrad, x, Δ)
    η = o.eta
    τ = o.tau

    g² = get!(
        o.acc,
        x,
        [fill(0.0, size(x)) for j = 1:o.n]
    )::Array{typeof(Tracker.data(x)), 1}
    i = get!(o.iters, x, 1)::Int

    # Example: suppose i = 12 and o.n = 10
    idx = mod(i - 1, o.n) + 1 # => idx = 2

    # set the current
    @inbounds @. g²[idx] = Δ^2 # => g²[2] = Δ^2 where Δ is the (o.n + 2)-th Δ

    # TODO: make more efficient and stable
    s = sum(g²)

    # increment
    o.iters[x] += 1

    # TODO: increment (but "truncate")
    # o.iters[x] = i > o.n ? o.n + mod(i, o.n) : i + 1

    @. Δ *= η / (τ + sqrt(s) + ϵ)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b425dc3c-89aa-11eb-29f0-271d5d515fb2
using IterTools: product

# ╔═╡ e80d9c13-4601-46b8-a6fe-6fc4cfddd31d
using Yao

# ╔═╡ 2f7f49f4-89b3-11eb-3bc3-ebc27d833a26
using LinearAlgebra: I

# ╔═╡ 1a3b87dd-b03e-44c9-8cac-906da4327396
md"""
The concepts to be introduced in this notebook, such as the Ising model, simulated annealing, and the transverse Ising model, play an important role in today's quantum algorithms and quantum computing paradigms, including quantum annealing, the quantum approximate optimization algorithm, and quantum-enhanced sampling. Here we give some insight on how these physical building blocks work.

# The Ising model

We would like to make a connection between the computational hardness of a problem and how difficult it is to solve a corresponding physical system. The Ising model is the most basic model to do this. It is an extensively studied model and one of the most basic examples to teach statistical mechanics and phase transitions, but we only require an elementary understanding of it.

Imagine that you have two magnets fixed on the same axis.
"""

# ╔═╡ b3d03fd8-89a6-11eb-31c3-d77d45414e64
begin
	import PlutoUI
	PlutoUI.LocalResource("../figures/two_magnets.svg", :width => 200, :alt => "Two magnets")
end

# ╔═╡ 31807b64-89a7-11eb-1957-9fd9c5e60f41
md"""
They will naturally anti-align: one will have north pole facing up, the second the south pole facing up. We can think of them as two binary variables, $\sigma_1$ and $\sigma_2$. Say, if the north pole is facing up, we assign the value +1 to the variable, and -1 otherwise. To abstract away from magnets, in general, we call these variables spins. So in the optimal configuration, their product is -1:

`` \sigma_1\sigma_2=-1 ``

We can think of this as the energy of the system: the lowest energy is called the ground state energy. Note that there are two physical configurations corresponding to this optimum: $\sigma_1=+1, \sigma_2=-1$, and $\sigma_1=-1, \sigma_2=+1$.

If we keep adding more magnets to this system, we can sum up their pairwise interaction to get the total energy. The total energy of the system is called the Hamiltonian, and we will denote it by $H$. So if we have $N$ magnets arranged along a straight line, we have

`` H=\sum_{i=1}^{N-1} \sigma_i \sigma_{i+1}``.

We did a simplification here: we assumed that remote magnets do not interact with each other (e.g. there is no such term as $\sigma_i\sigma_{i+2}$. In general, the interactions modeled depend on the layout of the spins and assumptions about the physical model: there will be some graph describing the connectivity of interactions. To reflect this, we write

`` H=\sum_{<i,j>} \sigma_i \sigma_{j}``,

where $<i,j>$ typically means nearest neighbours, but it is up to us to declare what nearest neighbours mean.

Now imagine that the distance is not the same between each pair. In other words, some pairs interact more than others. We can express this by adding a parameter that describes the interaction strength:

`` H=-\sum_{<i,j>} J_{ij} \sigma_i \sigma_{j}``,

where $J_{ij}$ is a real number. We added a negative sign to the Hamiltonian: this is by convention. If the spins are antiferromagnetic, that is, they behave as we would expect from magnets, then all $J_{ij}$ values would be negative. That cancels out the negative sign of the sum, so we still expect that each product $\sigma_i\sigma_j$ would give you -1 in the optimum configuration.

The model is fairly complicated by this point. Imagine that you have many spins and not all of them behave like magnets (that is, $J_{ij}$ can take both negative and positive values for different pairs). Nature still wants to find the lowest energy configuration, though. Let's take a look at how we would do it in code. Let's calculate the energy of spins on a line, given some couplings and a spin configuration:
"""

# ╔═╡ ca4a2698-7d5f-4480-9e6e-6bc4dd3f8a3a
calculate_energy(J, σ) = -sum(J_ij*σ[i]*σ[i+1] for (i,J_ij) in enumerate(J))

# ╔═╡ 0aff9cab-8f0d-42ee-ab52-cc96b0a05895
md"""
Let's give it a fixed set of couplings and a spin configuration on three sites:
"""

# ╔═╡ b84763a1-996a-4d50-ab7d-8126d758ce4b
begin
	J = [1.0, -1.0]
	σ = [+1, -1, +1]
end

# ╔═╡ d6241dd0-0540-40c9-8d7d-8554c182a1c0
md"""
The energy of this is
"""

# ╔═╡ 7972486c-eb23-461d-8998-ccb5474049f5
calculate_energy(J, σ)

# ╔═╡ f2133911-8380-47a1-823f-49f5425f5f15
md"""
Is this the ground state? How do we know? We are interested in the minimum, but we cannot use some gradient-based method to find it, since the variables are binary, plus the optimization landscape is nonconvex. So the easiest choice is an exhaustive search of all possibilities:
"""

# ╔═╡ bc02c9c9-3fc1-486e-b837-e6e3113f4069
let
	t = r =[]
	for σ in collect(product([+1,-1], [+1, -1], [+1, -1]))
	    t = [t;calculate_energy(J, σ)]
		r = [r; [i for i in σ]]
	end
	(t,r) = hcat(t), reshape(r, 8,3)
end

# ╔═╡ 6b89bdf8-7c27-4a87-8375-3e80aa4a2bdd
md"""
We see that -2 is the optimum, with two optimal configurations, but we had to enumerate all possibilities to figure this out. For this particular case, there are more clever ways to find the best solution, but in the general case, this is not the case.

To get to the general case, we need one more component, an external field. Imagine that you add a large magnet below each and every one of our magnets, creating an external magnetic field for each site. If this field is strong enough, it can override the pairwise interaction and flip the magnets. We model this by adding a linear term to the Hamiltonian:

`` H=-\sum_{<i,j>} J_{ij} \sigma_i \sigma_{j} - \sum_i h_i \sigma_i``,

where $h_i$ is the strength of the external field. This is the full description of the *classical Ising model*. The Hamiltonian describes the energy, but in computer science language, it means it expresses the objective function we want to minimize. The corresponding computer science problem is called quadratic unconstrained binary optimization (QUBO), where the only difference is that the variables take values in $\{0, 1\}$, but that is only a constant shift. QUBOs are NP-hard in general, that is, we are not aware of an efficient polynomial time algorithm to solve any given QUBO. So the generic strategy is the exhaustive search we did above, which takes exponentially many steps in the number of sites (variables).

As we mentioned, nature seeks the minimum energy configuration. So how does computational hardness map to physical difficulty? Imagine that the energy difference between the ground state and the next lowest energy state (also called the first excited state) is small, but the energetic cost of going from one to the other is high. A cartoon picture of this is the following:
"""

# ╔═╡ 41dd7fe8-89ac-11eb-3b00-a5df2682ae75
PlutoUI.LocalResource("../figures/energy_landscape.svg", :width => 300, :alt => "Energy Landscape")

# ╔═╡ 6e79b61e-89ac-11eb-238e-d5a2f77d0008
md"""
If we start from a random configuration, we might get stuck in the local optimum denoted by the green spot. This is what happens in metals if they are cooled down too quickly: the crystal lattice will have imperfections and the metal will not have the desired properties. A process called *annealing* helps in metallurgy: by increasing the temperature, the chance of overcoming the potential barrier increases and the crystal structure can reconfigure itself. If the barrier is high and the energy difference is small between the ground state and the first excited state, the probability of this happening drops. This is what it means that the problem is difficult to do in a physical system.

Annealing inspired a heuristic algorithm called *simulated annealing*. This defines a temperature to be able to hop out of local minima. The temperature is lowered over time to find the actual minimum. Simulated annealing has many implementations. Here we'll use the one implemented in dimod to solve our problem above:
"""

# ╔═╡ 54f1adc8-ef80-4801-86f5-835b1d499336
import dimod

# ╔═╡ 3adf3d6f-1d41-464c-a3d1-8558478600f3
md"""
The simulated annealing solver requires us to define the couplings as a dictionary between spins, and we must also pass the external field values as a dictionary. The latter is all zeros for us.
"""

# ╔═╡ 4405ace5-a383-4810-81b4-cae1fcb2222c
J = {(0, 1): 1.0, (1, 2): -1.0}
h = {0:0, 1:0, 2:0}

# ╔═╡ 5efc75f4-3de1-4b0e-b9f4-87e74cda9103
md"""
We instantiate an Ising model:
"""

# ╔═╡ 0ce517fd-1989-400a-bfe7-77ce043a7ac2
model = dimod.BinaryQuadraticModel(h, J, 0.0, dimod.SPIN)

# ╔═╡ 536f1777-09e3-48de-82a5-0058dc921b95
md"""
Finally, we create a simulated annealing sampler that pulls out potentially optimal solutions, and we read out 10 possible solutions:
"""

# ╔═╡ ee4c6fd4-7082-4a74-bcd5-9a101e6fb1e3
sampler = dimod.SimulatedAnnealingSampler()
response = sampler.sample(model, num_reads=10)

# ╔═╡ b9f62b85-8648-4c96-b571-69fe67ab4ccd
md"""
You can see that this configuration is actually easy, since you get the optimal solution -2 most of the time:
"""

# ╔═╡ cedb40a8-91b8-40e3-a37d-2e2da894dddc
[solution.energy for solution in response.data()].count(-2)

# ╔═╡ c0a1d5c7-7e97-4eb9-b423-543132776439
md"""
Simulated annealing is a classical heuristic algorithm. Quantum annealing uses physical effects to find the global optimum of an Ising model: it uses thermal effects just like annealing in metallurgy, but it also uses quantum effects like tunneling to overcome potential barriers.

The Ising model also plays an important role in quantum-enhanced sampling, but that idea requires a better understanding of the role of temperature, which we will revisit in a subsequent notebook.
"""

# ╔═╡ f0398f32-79da-4c74-b3f3-40f47d4a18f5
md"""
# The transverse-field Ising model
"""

# ╔═╡ 55ff9df7-574b-4963-992b-7905e2c24331
md"""
We discussed the Hamiltonian of the classical Ising model. We can write the same Hamiltonian in a quantum mechanical form. In quantum mechanics, the Hamiltonian is not a function of variables, but of operators. We will simulate what it means in a quantum circuit.
"""

# ╔═╡ 51f0a505-6608-435d-b28c-7ab2d2398e0a
md"""
The operator that replicates the effect of what we have seen in the classical case is the Pauli-Z matrix, defined as $\begin{bmatrix}1 & 0\\ 0& -1\end{bmatrix}$. Let's see what it does on the elements of the computational basis:
"""

# ╔═╡ 5cc78f38-5061-4212-b9dc-a2d740caf65f
let
	q = zero_state(1)
	state(q |> put(1, 1=>Z))
end

# ╔═╡ 529cba9b-cc11-4e20-b9a5-2eb432236ff2
md"""
This is nothing but the $|0\rangle$ state. In other words, it does not do anything to $|0\rangle$, which can also be thought of as multiplying it by +1. Let's try it on $|1\rangle$:
"""

# ╔═╡ d311efdc-8ea4-4579-950b-114d17d1b6ea
let
	q = zero_state(1)
	state(q |> chain(1, put(1=>X), put(1=>Z)))
end

# ╔═╡ a5c79e07-f541-4910-9bea-7553b7d4374c
md"""
We get $-|1\rangle$, which means it adds a minus sign to it. This way we have the +1, -1 values, just the same way as in the classical formalism. If we write $\sigma^Z_i$ for the operator $Z$ at a site $i$, the quantum mechanical Hamiltonian of the classical Ising model reads as

`` H=-\sum_{<i,j>} J_{ij} \sigma^Z_i \sigma^Z_{j} - \sum_i h_i \sigma^Z_i``.

Technically speaking, we should put a hat on $H$ and on all of the $\sigma^Z_i$ to indicate that they are operators, and not numbers or variables, but we omit this for notational simplicity.

The expectation value $<H>$ of the Hamiltonian is the energy of the system, and the corresponding quantum state $|\psi\rangle$ is the configuration of that energy level. We can create the quantum mechanical version of calculating the energy, matching the function we defined above for the classical mechanical variant:
"""

# ╔═╡ cb9c2854-2050-4767-9f2c-9499062f5539
calculate_energy_expectation(state, hamiltonian) = Float64(real(conj(state'ᵀ) * hamiltonian * state))

# ╔═╡ 39a8fa28-4750-496b-8f80-feb34f0fa876
md"""
It is a bit tricky to define the Hamiltonian with the $\sigma^Z_i$ operators, since saying that it acts on site $i$ means that it acts trivially on all other sites. So, for instance, for two sites, if we act on site one, the actual operator is $\sigma^Z\otimes I$, and acting on site two, we have $I \otimes \sigma^Z$. The above function to calculate the energy takes numpy arrays, so we manually define $\sigma^Z$ and calculate the energy of the Hamiltonian $H=-\sigma^Z_1\sigma^Z_2  - 0.5 (\sigma^Z_1 + \sigma^Z_2)$ on the state $|00\rangle$.
"""

# ╔═╡ 7b1acb49-4ea6-4b7e-9982-4c64a8cdd54f
begin
	PauliZ = [1 0; 0 -1]
	IZ = kron(Matrix(I, 2, 2), PauliZ)
	ZI = kron(PauliZ, Matrix(I, 2, 2))
	ZZ = kron(PauliZ, PauliZ)
	H = -ZZ + -0.5*(ZI+IZ)
	ψ = kron([1, 0], [1, 0])
	calculate_energy_expectation(ψ, H)
end

# ╔═╡ 132f9367-fabe-477e-a3f2-63d65b1cc1ae
md"""
This Hamiltonian commutes, which means all of its operators are commutative, which is a clear sign of nothing much quantum going on.

To make this a quantum Ising model, we need to add a term that does not commute with the rest of the terms. A transverse field is such, which is an on-site interaction just like the external field. Its effect is described by the Pauli-X operator (the NOT gate), which we will denote by $\sigma^X_i$ for a site $i$. It is very easy to see that the Pauli-Z and the Pauli-X do not commute:
"""

# ╔═╡ de97610f-33a2-4cae-a5b5-b3a0f56fa3dc
let
	q = zero_state(1)
	md"Pauli-X, then Pauli-Z:", state(q |> chain(1, put(1=>X), put(1=>Z)))
end

# ╔═╡ bc5307be-89b3-11eb-0083-11121b2d595e
let
	q = zero_state(1)
	md"Pauli-Z, then Pauli-X:", state(q |> chain(1, put(1=>Z), put(1=>X)))
end

# ╔═╡ c1920c45-c30b-40fb-9abc-58c5e176d891
md"""
There is a clear sign difference.

There are many other ways of making the Ising Hamiltonian noncommuting, but adding the onsite Pauli-X operations leads to the *transverse field Ising model*. Its Hamiltonian reads as

`` H=-\sum_{<i,j>} J_{ij} \sigma^Z_i \sigma^Z_{j} - \sum_i h_i \sigma^Z_i - \sum_i g_i \sigma^X_i``.

The transverse field Ising model is critically important to explain how quantum annealing works because by adding the $\sigma^X$ part to the Hamiltonian it becomes possible to exploit quantum effects like tunnelling. It is also important for understanding the quantum approximation optimization algorithms, since it was inspired by quantum annealing.
"""

# ╔═╡ Cell order:
# ╟─1a3b87dd-b03e-44c9-8cac-906da4327396
# ╟─b3d03fd8-89a6-11eb-31c3-d77d45414e64
# ╟─31807b64-89a7-11eb-1957-9fd9c5e60f41
# ╠═ca4a2698-7d5f-4480-9e6e-6bc4dd3f8a3a
# ╟─0aff9cab-8f0d-42ee-ab52-cc96b0a05895
# ╠═b84763a1-996a-4d50-ab7d-8126d758ce4b
# ╟─d6241dd0-0540-40c9-8d7d-8554c182a1c0
# ╠═7972486c-eb23-461d-8998-ccb5474049f5
# ╟─f2133911-8380-47a1-823f-49f5425f5f15
# ╠═b425dc3c-89aa-11eb-29f0-271d5d515fb2
# ╠═bc02c9c9-3fc1-486e-b837-e6e3113f4069
# ╟─6b89bdf8-7c27-4a87-8375-3e80aa4a2bdd
# ╟─41dd7fe8-89ac-11eb-3b00-a5df2682ae75
# ╟─6e79b61e-89ac-11eb-238e-d5a2f77d0008
# ╠═54f1adc8-ef80-4801-86f5-835b1d499336
# ╟─3adf3d6f-1d41-464c-a3d1-8558478600f3
# ╠═4405ace5-a383-4810-81b4-cae1fcb2222c
# ╟─5efc75f4-3de1-4b0e-b9f4-87e74cda9103
# ╠═0ce517fd-1989-400a-bfe7-77ce043a7ac2
# ╟─536f1777-09e3-48de-82a5-0058dc921b95
# ╠═ee4c6fd4-7082-4a74-bcd5-9a101e6fb1e3
# ╟─b9f62b85-8648-4c96-b571-69fe67ab4ccd
# ╠═cedb40a8-91b8-40e3-a37d-2e2da894dddc
# ╟─c0a1d5c7-7e97-4eb9-b423-543132776439
# ╟─f0398f32-79da-4c74-b3f3-40f47d4a18f5
# ╟─55ff9df7-574b-4963-992b-7905e2c24331
# ╠═e80d9c13-4601-46b8-a6fe-6fc4cfddd31d
# ╟─51f0a505-6608-435d-b28c-7ab2d2398e0a
# ╠═5cc78f38-5061-4212-b9dc-a2d740caf65f
# ╟─529cba9b-cc11-4e20-b9a5-2eb432236ff2
# ╠═d311efdc-8ea4-4579-950b-114d17d1b6ea
# ╟─a5c79e07-f541-4910-9bea-7553b7d4374c
# ╠═cb9c2854-2050-4767-9f2c-9499062f5539
# ╟─39a8fa28-4750-496b-8f80-feb34f0fa876
# ╠═2f7f49f4-89b3-11eb-3bc3-ebc27d833a26
# ╠═7b1acb49-4ea6-4b7e-9982-4c64a8cdd54f
# ╟─132f9367-fabe-477e-a3f2-63d65b1cc1ae
# ╠═de97610f-33a2-4cae-a5b5-b3a0f56fa3dc
# ╠═bc5307be-89b3-11eb-0083-11121b2d595e
# ╟─c1920c45-c30b-40fb-9abc-58c5e176d891

max_stars_repo_name: VarLad/qml-mooc

max_issues_repo_name: VarLad/qml-mooc

max_forks_repo_name: VarLad/qml-mooc


text:
using DecisionMakingUtils
using OnlineStats
using Test

@testset "ZeroOneNormalization Test" begin
    k = 3
    T = Float32
    low = -collect(T, 1:k)
    high = collect(T, 1:k)
    ranges = hcat(low, high)
    n = ZeroOneNormalization(ranges)
    x = zeros(T, k)
    buff = zero(x)
    y = n(x)
    @test typeof(y) == typeof(x)
    @test size(y) == size(x)

    @test all(y .≥ 0.0)
    @test all(y .≤ 1.0)
    y = n(low)
    @test all(y .≥ 0.0)
    @test all(y .≤ 1.0)
    y = n(high)
    @test all(y .≥ 0.0)
    @test all(y .≤ 1.0)
    n(buff, high)
    @test all(isapprox.(buff, y))

    @. low -= 1
    @. high += 1
    n(buff, low)
    @test all(buff .< 0.0)
    n(buff, low, fit=true)
    @test all(buff .≥ 0.0)
    n(buff, high)
    @test all(buff .> 1.0)
    n(buff, high, fit=true)
    @test all(buff .≤ 1.0)
end

@testset "PosNegNormalization Test" begin
    k = 3
    T = Float32
    low = -collect(T, 1:k)
    high = collect(T, 1:k)
    ranges = hcat(low, high)
    n = PosNegNormalization(ranges)

    x = zeros(T, k)
    buff = zero(x)
    y = n(x)
    @test typeof(y) == typeof(x)
    @test size(y) == size(x)
    @test all(y .≥ -1.0)
    @test all(y .≤ 1.0)
    y = n(low)

    @test all(y .≥ -1.0)
    @test all(y .≤ 1.0)
    y = n(high)
    @test all(y .≥ -1.0)
    @test all(y .≤ 1.0)
    n(buff, high)
    @test all(isapprox.(buff, y))

    @. low -= 1
    @. high += 1
    n(buff, low)
    @test all(buff .< -1.0)
    n(buff, low, fit=true)
    @test all(buff .≥ -1.0)
    n(buff, high)
    @test all(buff .> 1.0)
    n(buff, high, fit=true)
    @test all(buff .≤ 1.0)
end

@testset "GaussianNormalization Test" begin
    k = 3
    T = Float32
    x1 = collect(T, 1:k)
    x2 = collect(T, 1:k) .* 2
    x3 = collect(T, 1:k) .* 3

    n = GaussianNormalization(k)

    x = zeros(T, k)
    buff = zero(x)
    y = n(x)
    @test typeof(y) == typeof(x)
    @test size(y) == size(x)

    @test all(isapprox.(y, 0.0))
    @. x = 1
    y = n(x)
    @test all(isapprox.(y, 1.0))

    n(buff, x)
    @test all(isapprox.(buff, y))

    v = k*KahanVariance()
    x1 = collect(T, 1:k)
    x2 = collect(T, 1:k) .* 2
    x3 = collect(T, 1:k) .* 3

    fit!(v, x1)
    @. y = (x1 - mean.(v.stats)) / std.(v.stats)
    n(buff, x1, fit=true)
    @test all(isapprox(buff, y))

    fit!(v, x2)
    @. y = (x2 - mean.(v.stats)) / std.(v.stats)
    n(buff, x2, fit=true)
    @test all(isapprox(buff, y))

    fit!(v, x3)
    @. y = (x3 - mean(v.stats)) / std(v.stats)
    n(buff, x3, fit=true)
    @test all(isapprox(buff, y))
end


text:
using Test, HusimiFunctions

r0 = SVector(0.5, 0.5)
k0 = SVector(0.5, 0.25)

using Random
Random.seed!(5)
r = [rand(SVector{2}) for i in 1:10000]
ψ = [rand(ComplexF64) for i in 1:10000]

Q1 = husimi(ψ, r, r0, k0, 2.0)
@test Q1/1e6 ≈ 1.882235 atol=1e-6
Q2 = husimi(ψ, r, r0, k0, 2.0, 5.0)
@test Q2/1e6 ≈ 1.881620 atol = 1e-6

# # non magnetic
# @btime husimi(ψ, r, r0, k0, 2.0; thread = false)
# @btime husimi(ψ, r, r0, k0, 2.0; thread = true)
#
# # non magnetic
# @btime husimi(ψ, r, r0, k0, 2.0, 0.5; thread = false)
# @btime husimi(ψ, r, r0, k0, 2.0, 0.5; thread = true)


text:
#=
This file implements the algorithm described in:
"A Schur-Parlett Algorithm for Computing Matrix Functions"
(Davies and Higham, 2003)
(with some modifications to handle real matrices more efficiently).
=#

#=
blockpattern groups eigenvalues in sets, implementing Algorithm 4.1 of the paper.
The function returns S and p. S is the pattern, where S[i] = s means that
vals[i] has been assigned to set s. p is the number of sets identified.
=#
delta = 0.1
function blockpattern(vals::Vector{Complex128}, schurtype::Type)::Tuple{Vector{Int64}, Int64}
	unassigned = -1
	S = fill(unassigned, length(vals))
	p = 0

	# assign vals[i] to set s.
	function assign(i::Int64, s::Int64)
		lambda = vals[i]
		for k = i:length(vals)
			if vals[k] == lambda
				S[k] = s
			end
		end
	end
	function mergesets(s::Int64, t::Int64)
		if s == t
			return
		end
		smin, smax = minmax(s, t)
		for k = 1:length(vals)
			if S[k] == smax
				S[k] = smin
			elseif S[k] > smax
				S[k] -= 1
			end
		end
		p -= 1
	end

	for i = 1:length(vals)
		if S[i] == unassigned
			p += 1
			assign(i, p)
		end
		for j = i+1:length(vals)
			if abs(vals[i]-vals[j]) <= delta
				if S[j] == unassigned
					assign(j, S[i])
				else
					mergesets(S[i], S[j])
				end
			end
		end
	end

	if schurtype <: Real
		# complex conjugate eigenvalues can't be separated in the
		# real Schur factorization, so we merge the sets involving them:
		i = 1
		while i < length(vals)
			if imag(vals[i]) != 0
				mergesets(S[i], S[i+1])
				i += 2
			else
				i += 1
			end
		end
	end
	return S, p
end

# analogous to ordschur/trsen:
function ordvec!(v::Vector{T}, select::B) where {T, B<:Union{Vector{Bool}, BitVector}}
	ilst = 1
	for ifst = 1:length(v)
		if select[ifst]
			if ifst != ilst
				v[ilst], v[ilst+1:ifst] = v[ifst], v[ilst:ifst-1]
			end
			ilst += 1
		end
	end
end

#=
reorder reorders the Schur decomposition (T, Q, vals) according to pattern S,
using the swapping strategy described in Algorithm 4.2.
The returned vector, blockend, contains the indices at which each block ends
(after the reordering).
=#
function reorder!(T::Matrix{N}, Q::Matrix{N}, vals::Vector{Complex{R}}, S::Vector{Int64}, p::Int64)::Vector{Int64} where {
	R<:Real, N<:Union{R, Complex{R}}}

	# for each set, calculate its mean position in S:
	pos = zeros(Float64, p)
	cou = zeros(Int64, p)
	for i = 1:length(S)
		set = S[i]
		pos[set] += i
		cou[set] += 1
	end
	pos ./= cou

	blockend = zeros(Int64, p)
	for set = 1:p
		numordered = (set == 1) ? 0 : blockend[set-1]
		minset = indmin(pos)
		select = [i <= numordered || S[i] == minset for i = 1:length(S)]
		ordschur!(T, Q, select)
		ordvec!(vals, select)
		ordvec!(S, select)
		blockend[set] = count(select)
		pos[minset] = Inf
	end
	return blockend
end

#=
taylorf computes f(T) using Taylor as described in Algorithm 2.6.
=#
function taylorf(f::Func, T::Mat, shift::N)::Mat where {
	Func, N<:Number, Mat<:Union{Matrix{N}, UpperTriangular{N, Matrix{N}}}}

	maxiter = 300
	lookahead = 10
	tay = f(shift + Taylor1(N, 20+lookahead))
	M = T - shift*Mat(eye(T))
	P = M
	F = tay.coeffs[1]*Mat(eye(T))
	for k = 1:maxiter
		needorder = k + lookahead
		@assert needorder <= tay.order + 1
		if needorder > tay.order
			tay = f(shift + Taylor1(N, min(maxiter+lookahead, tay.order*2)))
		end

		Term = tay.coeffs[k+1] * P
		F += Term
		normP = norm(P, Inf)
		P *= M
		small = eps()*norm(F, Inf)
		if norm(Term, Inf) <= small
			#=
			The termination condition is loosely inspired to the one in the paper.
			We check that the next lookahead terms in the series are small,
			estimating ||M^(k+r)|| with max(||M^k||, ||M^(k+1)||).
			Works well in practice, at least for exp, log and sqrt.
			=#
			delta = maximum(abs.(tay.coeffs[k+2:k+1+lookahead]))
			if delta*max(normP, norm(P, Inf)) <= small
				return F
			end
		end
	end
	error("Taylor did not converge.")
end

#=
blockf computes f(T) where T is an "atomic" block.
=#
function blockf(f::Func, T::Matrix{N}, vals::Vector{Complex{R}})::Matrix{N} where {
	Func, R<:Real, N<:Union{R, Complex{R}}}

	n = size(T, 1)
	if n == 1
		return f.(T)
	end
	if N <: Complex
		# T is complex and triangular with one cluster of eigenvalues.
		return taylorf(f, UpperTriangular(T), mean(vals))
	end
	if all(imag(vals) .== 0)
		# T is real and triangular with one cluster of real eigenvalues.
		return taylorf(f, UpperTriangular(T), mean(real(vals)))
	end
	if any(abs.(imag(vals)) .<= delta/2)
		# T is real and quasi-triangular with one cluster of complex
		# conjugate eigenvalues (and possibly some real eigenvalue).
		return taylorf(f, T, mean(real(vals)))
	end
	if n == 2
		#=
		T is real and 2x2 with two well-separated complex conjugate eigenvalues.
		We use the Hermite interpolating polynomial obtained with the
		Lagrange-Hermite formula, simplified assuming f(conj(x)) == conj(f(x)).
		=#
		v1, v2 = vals[1], vals[2]
		@assert v1 == conj(v2)
		psi1 = f(v1)/(v1 - v2)
		return 2.0*real(psi1*(T - [v2 0; 0 v2]))
	end
	#=
	T is real with two well-separated, complex conjugate clusters of eigenvalues.
	We use the complex Schur factorization to further divide T.
	Doing complex Schur at the block level is more efficient than using complex
	arithmetic for the whole matrix.
	=#
	U, Z, vals = schur(Matrix{Complex{R}}(T))::Tuple{Matrix{Complex{R}}, Matrix{Complex{R}}, Vector{Complex{R}}}
	select = imag(vals) .> 0
	ordschur!(U, Z, select)
	ordvec!(vals, select)
	@assert count(select) == n÷2
	F = Z*recf(f, U, vals, [n÷2, n])*Z'
	realF = real(F)
	if norm(imag(F), Inf) > sqrt(eps())*norm(realF, Inf)
		warn("T is real but f(T) has non-negligible imaginary part.")
		# note: we need f(T) to be real, because trsyl can't handle
		# quasi-triangular complex matrices.
	end
	return realF
end

#=
recf computes f(T) using the Parlett recurrence (T is block-triangular).
The paper suggests to do the recurrence iteratively (Algorithm 5.1),
we do it recursively as this version performs better.
=#
function recf(f::Func, T::Matrix{N}, vals::Vector{Complex{R}}, blockend::Vector{Int64})::Matrix{N} where {
	Func, R<:Real, N<:Union{R, Complex{R}}}

	@assert length(blockend) > 0
	if length(blockend) == 1
		return blockf(f, T, vals)
	end

	# split T in 2x2 superblocks of size ~n/2:
	n = size(T, 1)
	b = indmin(abs.(blockend .- n/2))
	bend = blockend[b]
	T11, T12 = T[1:bend, 1:bend], view(T, 1:bend, bend+1:n)
	T21, T22 = view(T, bend+1:n, 1:bend), T[bend+1:n, bend+1:n]

	F11 = recf(f, T11, vals[1:bend], blockend[1:b])
	F22 = recf(f, T22, vals[bend+1:n], blockend[b+1:end] .- bend)

	# T11*F12 - F12*T22 = F11*T12 - T12*F22
	F12, scale = LAPACK.trsyl!('N', 'N', T11, T22, F11*T12 - T12*F22, -1)

	return [F11 F12/scale; zeros(T21) F22]
end

#=
Computes f(A) using the Schur-Parlett algorithm.
When A is a real matrix, computation will be done mostly in real arithmetic
and the algorithm will assume f(conj(x)) == conj(f(x)).
=#
function schurparlett(f::Func, A::Matrix{N})::Matrix{N} where {Func, N<:Union{Float64, Complex128}}
	T, Q, vals = schur(A)
	return schurparlett(f, T, Q, vals)
end

#=
Analogous to schurparlett(f, A), but A is provided as its Schur decomposition.
=#
function schurparlett(f::Func, T::Matrix{Float64}, Q::Matrix{Float64}, vals::Vector{Float64})::Matrix{Float64} where {Func}
	return schurparlett(f, T, Q, Vector{Complex128}(vals))
end
function schurparlett(f::Func, T::Matrix{N}, Q::Matrix{N}, vals::Vector{Complex128})::Matrix{N} where {
	Func, N<:Union{Float64, Complex128}}

	if size(T, 1) != size(T, 2)
		throw(DimensionMismatch("T must be square"))
	end
	if size(Q) != size(T) || length(vals) != size(T, 1)
		throw(DimensionMismatch("T, Q, vals dimension mismatch"))
	end
	if size(T, 1) == 0
		error("T is empty")
	end

	d = diag(T)
	D = diagm(d)
	if norm(T-D, Inf) <= 1000*eps()*norm(D, Inf) # T is diagonal
		return Q*diagm(f.(d))*Q'
	end

	S, p = blockpattern(vals, N)
	blockend = reorder!(T, Q, vals, S, p)
	return Q*recf(f, T, vals, blockend)*Q'
end


text:
#Basic function to use

"""
    absolute_error(x, x̂)

return the Absolute error of x and x̂
"""
absolute_error(x, x̂) = abs(x .- x̂)

"""
    relative_error(x, x̂)

The relative error between X1 and X2
"""
relative_error(x, x̂) = abs((x .- x̂) ./ x)

"""
    ξ(x₀, x)

return the number ξ(x) bwteen x₀ and x.
"""
@inline ξ(x₀, x) = rand(1)[1] * (x - x₀)/x + x₀

"""
    NthDrivative(f::Function, x::Real, n::Int)

return the nth derivative value at x
"""
@inline function NthDerivative(f::Function, x::Real, n::Int)
    if n == 0
        return f(x)
    else
        return ForwardDiff.derivative(x->NthDerivative(f, x, n-1), x)
    end
end

"""
    TaylorPolynomials(f::Function, x::Real, x₀::Real, n::Int)

return the value of nth Taylor Ploynomial at x.
"""
@inline function TaylorPolynomials(f::Function, x::Real, x₀::Real, n::Int)
    result = []
    for i in 0:n
        push!(result, NthDerivative(f, x₀, i)*(x - x₀)^i/factorial(i))
    end
    return sum(result)

end


text: # Part 2: another PDE


text:
using Dates

struct Price
    pipsize::Rational{Int64}
    pips::Int64
end

struct Tick
    time::DateTime
    ask::Price
    bid::Price
    volume::Int64
end

struct Position
    symbol::String
    costbasis::Float64
    size::Int64
end

struct Account
    cash::Price
    positions::Set{Position}
end

struct Indicator
    inputs::Vector{Union{Float64,Nothing}}
    input::String
    value::Function
end

struct Engine
    acct::Account
    ts::Vector{Tick}
    index::Int64
    indicators::Dict{String,Indicator}
end