text:
module Extraction
using Logging
using LightGraphs, MetaGraphs

import JSON

# TODO: Create functions for each
# Cause - is concept label, Effect - definition label
# within a RelationMention with labels "Definition, Entity"

"""    definitiongraph(dir::String, namefunc)

read a directory of json files and ingest all the Automates matches into a metagraph.

dir is the directory containing the json files and
namefunc is a function to convert objects into vertex names.

See also: sequentialnamer
"""
function definitiongraph(dir::String, namefunc)
    files = [joinpath(dir, f) for f in readdir(dir)]
    @info "Processing files: $files"
    return definitiongraph(files, namefunc)
end

"""    term(record)

get the text of a definition's term
"""
function term(record)
    return record["arguments"]["cause"][1]["text"]
end

"""    definition(record)

get the text of a definition's body
"""
function definition(record)
    return record["arguments"]["effect"][1]["text"]
end

"""    sequentialnamer(prefix="")

A closure used to get sequential vertex names from a stream of strings.
"""
function sequentialnamer(prefix="")
    i = 0
    function genname(s)
        i += 1
        return "$prefix$i"
    end
    return genname
end

function definitiongraph(files::Vector{String}, namefunc)
    graph = DiGraph()
    changes_recorded = []
    v_count = 0
    for def_file in files
        @info "Reading Definitions from: $def_file"
        json_res = JSON.parsefile(def_file)
        mentions = json_res["mentions"]

        for val in mentions

            if haskey(val,"type") && val["type"] == "RelationMention"
                @debug "Parsed JSON record" json=val
                v_count += 1
                cause = term(val)
                effect = definition(val)
                @info "Found Definition" term=cause definition=effect
                @debug("Add new vertex for Concept to Graph")
                add_vertex!(graph)
                push!(changes_recorded, Dict(:v_id=>v_count,
                                             :type=> "AddVertex",
                                             :values => Dict(:name=> string("Term_", namefunc(effect)),
                                                             :id => v_count,
                                                             :text => cause)))
                v_count += 1
                add_vertex!(graph)
                push!(changes_recorded, Dict(:v_id=>v_count,
                                             :type=> "AddVertex",
                                             :values => Dict(:name=> string("Def_",namefunc(effect)),
                                                             :id => v_count,
                                                             :text => effect)))
                @debug("Adding edge from Concept to Definition")
                add_edge!(graph, v_count-1, v_count)
                push!(changes_recorded, Dict(:edge=> Edge(v_count-1,v_count),
                                             :type=> "AddEdge",
                                             :values=> Dict(:name=> "is defined as")))
            end
        end
    end

    # Create MetaGraph after digraph is created using changes record
    metagraph = MetaDiGraph(graph)
    for change in changes_recorded
        if change[:type] == "AddEdge"
            set_props!(metagraph, change[:edge], change[:values])
            @info("Adding Edge",edge=change[:edge], props(metagraph, change[:edge])...)
        end
        if change[:type] == "AddVertex"
            set_props!(metagraph, change[:v_id], change[:values])
            @info("Adding Vertex",vertex=change[:v_id], props(metagraph, change[:v_id])...)
        end
    end
    set_indexing_prop!(metagraph, :name)
    return metagraph
end

end


text:
module Errors

using JuLIP: Atoms, energy, forces, virial
using IPFitting: LsqDB, Dat, configtype, err_weighthook, eval_obs,
                      observation, hasobservation, observations,
                      vec_obs, tfor_observations, truncate_string
using IPFitting.DB: matrows
using FileIO, Printf, Base.Threads

using Statistics: quantile
using LinearAlgebra: norm
using ProgressMeter, DataFrames, PrettyTables
using Random: shuffle!

export add_fits!, add_fits_serial!, rmse, mae, rmse_table, rmse_table_wg, mae_table,
       lsqerrors, splittraintest

# ------------------------- USER INTERFACE FUNCTIONS -------------------------

rmse(data::Vector{Dat}; fitkey="fit") = _fiterrs(data; fitkey=fitkey, p=2)
mae(data::Vector{Dat}; fitkey="fit") = _fiterrs(data; fitkey=fitkey, p=1)

rmse(D::Dict) = D["rmse"], D["relrmse"]
mae(D::Dict) = D["mae"], D["relmae"]

rmse_table(data::AbstractVector{Dat}) = rmse_table(rmse(data)...)
mae_table( data::AbstractVector{Dat}) =  mae_table( mae(data)...)

rmse_table(D::Dict) = rmse_table(rmse(D)...)

rmse_table_wg(D::Dict, weights) = rmse_table_wg(rmse(D)..., weights)

mae_table(D::Dict) =  mae_table( mae(D)...)

rmse_table(errs::Dict, errs_rel::Dict; configtypes=:) =
   _err_table(errs, errs_rel, "RMSE", configtypes)

rmse_table_wg(errs::Dict, errs_rel::Dict, weights; configtypes=:) =
      _err_table_wg(errs, errs_rel, "RMSE WEIGHTED TABLE", weights, configtypes)

mae_table(errs::Dict, errs_rel::Dict) = _err_table(errs, errs_rel, "MAE")

"""
`add_fits(IP, data::Vector{Dat}) -> nothing`

for each `d in data`, compute the observations with the IP (fitted to this
data) and store in `d.info["fit"][...]`

The key `"fit"` can be replaced with any string, through the kwarg
`fitkey`, e.g.,
```julia
data = ...
add_fits!(PIP4Benv, data; fitkey = "PIP4Benv")
add_fits!(GAP, data; fitkey = "GAP2010")
```
"""
function add_fits!(IP, configs::Vector{Dat}; fitkey = "fit")
   # create the nec essary dictionaries
   for cfg in configs
      cfg.info[fitkey] = Dict{String, Vector{Float64}}()
   end
   tfor_observations( configs,
         (n, okey, cfg, lck) -> begin
            obs = vec_obs(okey, eval_obs(okey, IP, cfg))
            lock(lck)
            cfg.info[fitkey][okey] = obs
            unlock(lck)
         end;
         msg="Add Fit info to configs")

   return nothing
end

function add_fits_serial!(IP, configs::Vector{Dat}; fitkey = "fit")
   # create the nec essary dictionaries
   for cfg in configs
      cfg.info[fitkey] = Dict{String, Vector{Float64}}()
   end
   for (okey, cfg, _) in observations(configs)
      obs = vec_obs(okey, eval_obs(okey, IP, cfg))
      cfg.info[fitkey][okey] = obs
   end
   return nothing
end

function _relerr(relerrs, ct, ot)
   if haskey(relerrs[ct], ot)
      if relerrs[ct][ot] isa Number
         return 100*relerrs[ct][ot]
      end
   end
   return NaN
end

function _err(errs, ct, ot)
   if  haskey(errs[ct], ot)
      if errs[ct][ot] isa Number
         return errs[ct][ot]
      end
   end
   return NaN
end

# function _err_table_wg(errs, relerrs, title, weights, configtypes=:)
#    lentitle = length(title)
#    print("┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n")
#    print("┃ "); printstyled(title; bold=true);
#    print(repeat(' ', 70-3-lentitle)); print("┃\n")
#    print("┣━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━┫\n")
#    print("┃ config type  ┃     E [meV]     │     F [eV/A]    │     V [meV]     ┃\n")
#    print("┠──────────────╂────────┬────────┼────────┬────────┼────────┬────────┨\n")
#    s_set = ""
#    for ct in keys(errs)  # ct in configtypes
#       if (configtypes != :)
#          if !(ct in configtypes)
#             continue
#          end
#       end
#
#       if haskey(weights, ct)
#          w = weights[ct]
#       else
#          w = weights["default"]
#       end
#
#       s = @sprintf("┃ %12s ┃ %6.2f ┊ %5.3f%% │ %6.3f ┊ %5.2f%% │ %6.1f ┊ %5.2f%% ┃\n",
#          truncate_string(ct, 12),
#          _err(errs, ct, "E")*1000 * w["E"], _relerr(relerrs, ct, "E") * w["E"],
#          _err(errs, ct, "F") * w["F"], _relerr(relerrs, ct, "F") * w["F"],
#          _err(errs, ct, "V")*1000 * w["V"], _relerr(relerrs, ct, "V") * w["V"])
#       if ct == "set"
#          s_set = s
#       else
#          print(s)
#       end
#    end
#    if s_set != ""
#       print("┠──────────────╂────────┼────────┼────────┼────────┼────────┼────────┨\n")
#       print(s_set)
#    end
#       print("┗━━━━━━━━━━━━━━┻━━━━━━━━┷━━━━━━━━┷━━━━━━━━┷━━━━━━━━┷━━━━━━━━┷━━━━━━━━┛\n")
# end

function _err_table(errs, relerrs, title, configtypes=:)
   #@show relerrs
   lentitle = length(title)
   datf = DataFrame()
   datf[!, Symbol("cf")] = Any[]
   datf[!, Symbol("E")]= Any[]
   datf[!, Symbol("F")]= Any[]
   datf[!, Symbol("Fr")]= Any[]
   datf[!, Symbol("V")]= Any[]

   #need to add relative errors!

   for ct in sort(collect(keys(errs)))
      if ct != "set"
         push!(datf, [ct, _err(errs, ct, "E")*1000, _err(errs, ct, "F"), _err(relerrs, ct, "F")*100, _err(errs, ct, "V")*1000])
      end
   end

   push!(datf, ["set", _err(errs, "set", "E")*1000, _err(errs, "set", "F"), _err(relerrs, "set", "F")*100, _err(errs, "set", "V")*1000])

   header = ["config type" "E [meV]" "F [eV/A]" "F [%]" "V[meV]"]

   pretty_table(datf, header, formatters = ft_printf("%5.3f"), body_hlines = [length(datf[!, 1])-1])

   # print("┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓\n")
   # print("┃ "); printstyled(title; bold=true);
   # print(repeat(' ', 70-3-lentitle)); print("┃\n")
   # print("┣━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━━━━━┫\n")
   # print("┃ config type  ┃     E [meV]     │     F [eV/A]    │     V [meV]     ┃\n")
   # print("┠──────────────╂────────┬────────┼────────┬────────┼────────┬────────┨\n")
   # s_set = ""
   # for ct in keys(errs)  # ct in configtypes
   #    if (configtypes != :)
   #       if !(ct in configtypes)
   #          continue
   #       end
   #    end
   #
   #    s = @sprintf("┃ %12s ┃ %6.2f ┊ %5.3f%% │ %6.3f ┊ %5.2f%% │ %6.1f ┊ %5.2f%% ┃\n",
   #       truncate_string(ct, 12),
   #       _err(errs, ct, "E")*1000, _relerr(relerrs, ct, "E"),
   #       _err(errs, ct, "F"), _relerr(relerrs, ct, "F"),
   #       _err(errs, ct, "V")*1000, _relerr(relerrs, ct, "V"))
   #    if ct == "set"
   #       s_set = s
   #    else
   #       print(s)
   #    end
   # end
   # if s_set != ""
   #    print("┠──────────────╂────────┼────────┼────────┼────────┼────────┼────────┨\n")
   #    print(s_set)
   # end
   #    print("┗━━━━━━━━━━━━━━┻━━━━━━━━┷━━━━━━━━┷━━━━━━━━┷━━━━━━━━┷━━━━━━━━┷━━━━━━━━┛\n")
end

# ---------------------------------------------------------------------
# ------------------------- PRIVATE FUNCTIONS -------------------------
# ---------------------------------------------------------------------

"""
initialize an error dictionary
"""
function errdict(configtypes)
   D = Dict{String, Dict{String, Float64}}()
   for ct in configtypes
      D[ct] = Dict{String, Float64}()
   end
   D["set"] = Dict{String, Float64}()
   return D
end

function _fiterrs(configs::Vector{Dat}; fitkey="fit", p=2)
   configtypes = unique(configtype.(configs))
   err = errdict(configtypes)
   relerr = errdict(configtypes)
   nrm = errdict(configtypes)
   len = errdict(configtypes)

   # assemble the errors
   for (okey, dat, _) in observations(configs)
      ct = configtype(dat)
      if !haskey(err[ct], okey)
         len[ct][okey] = 0.0
         nrm[ct][okey] = 0.0
         err[ct][okey] = 0.0
      end
      if !haskey(err["set"], okey)
         len["set"][okey] = 0.0
         nrm["set"][okey] = 0.0
         err["set"][okey] = 0.0
      end
      w = err_weighthook(okey, dat)  # the default weight for this observation
      exval = dat.D[okey] * w
      fitval = dat.info[fitkey][okey] * w
      len[ct][okey] += length(exval)
      len["set"][okey] += length(exval)
      nrm[ct][okey] += norm(exval, p)^p
      nrm["set"][okey] += norm(exval, p)^p
      err[ct][okey] += norm(exval-fitval, p)^p
      err["set"][okey] += norm(exval-fitval, p)^p
   end

   # normalise correctly
   for ct in keys(err)
      for okey in keys(err[ct])
         err[ct][okey] = (err[ct][okey]/len[ct][okey])^(1/p)
         nrm[ct][okey] = (nrm[ct][okey]/len[ct][okey])^(1/p)
         relerr[ct][okey] = err[ct][okey] / nrm[ct][okey]
      end
   end

   return err, relerr
end

# TODO: the next fucntion is to be retired!!

@noinline function lsqerrors(db, c, Ibasis; cfgtypes = Colon(), Vref = nothing,
                                            Icfg=Colon())

   cfgtypes isa Colon ?  cfgtypes = configtypes(db) :
                         cfgtypes = collect(cfgtypes)

   # create the dict for the fit errors
   errs = Dict("rmse" => errdict([cfgtypes; "set"]),
               "relrmse" => errdict([cfgtypes; "set"]),
               "mae"  => errdict([cfgtypes; "set"]),
               "maxe" => errdict([cfgtypes; "set"]),
               "nrm2" => errdict([cfgtypes; "set"]),
               "nrminf" => errdict([cfgtypes; "set"]),
               "nrm1" => errdict([cfgtypes; "set"]) )
   lengths = Dict{String, Dict{String, Int}}()
   for ct in [cfgtypes; "set"]
      lengths[ct] = Dict{String, Int}()
   end

   for (okey, dat, icfg) in observations(db, Icfg)

      ct = configtype(dat)
      if !(ct in cfgtypes)
         continue
      end
      if !haskey(errs["rmse"][ct], okey)
         errs["rmse"][ct][okey] = 0.0
         errs["nrm2"][ct][okey] = 0.0
         errs["mae"][ct][okey] = 0.0
         errs["nrm1"][ct][okey] = 0.0
         errs["maxe"][ct][okey] = 0.0
         errs["nrminf"][ct][okey] = 0.0
         lengths[ct][okey] = 0
      end
      if !haskey(errs["rmse"]["set"], okey)
         errs["rmse"]["set"][okey] = 0.0
         errs["nrm2"]["set"][okey] = 0.0
         errs["mae"]["set"][okey] = 0.0
         errs["nrm1"]["set"][okey] = 0.0
         errs["maxe"]["set"][okey] = 0.0
         errs["nrminf"]["set"][okey] = 0.0
         lengths["set"][okey] = 0
      end
      # assemble the observation vector
      y = Float64[]
      o = observation(dat, okey)
      # correct the observation if we are fitting from a reference potential
      if Vref != nothing
         o -= vec_obs(okey, eval_obs(okey, Vref, dat))
      end
      append!(y, o)
      # get the lsq block
      irows = matrows(dat, okey)
      block = db.Ψ[irows, Ibasis]
      # the error on this particular data-point
      e = block * c - y
      # compute the various errors for this ct/okey combination
      w = err_weighthook(okey, dat)
      errs["rmse"][ct][okey] += norm(w .* e)^2
      errs["nrm2"][ct][okey] += norm(w .* y)^2
      errs["mae"][ct][okey]  += norm(w .* e, 1)
      errs["nrm1"][ct][okey] += norm(w .* y, 1)
      errs["maxe"][ct][okey] = max(errs["maxe"][ct][okey], norm(e, Inf))
      errs["nrminf"][ct][okey] = max(errs["nrminf"][ct][okey], norm(y, Inf))
      lengths[ct][okey] += length(y)
      errs["rmse"]["set"][okey] += norm(w .* e)^2
      errs["nrm2"]["set"][okey] += norm(w .* y)^2
      errs["mae"]["set"][okey]  += norm(w .* e, 1)
      errs["nrm1"]["set"][okey] += norm(w .* y, 1)
      errs["maxe"]["set"][okey] = max(errs["maxe"]["set"][okey], norm(e, Inf))
      errs["nrminf"]["set"][okey] = max(errs["nrminf"]["set"][okey], norm(y, Inf))
      lengths["set"][okey] += length(y)
   end

   for ct in [cfgtypes; "set"]
      for okey in keys(errs["rmse"][ct])
         len = lengths[ct][okey]
         errs["rmse"][ct][okey] = sqrt( errs["rmse"][ct][okey] / len )
         errs["nrm2"][ct][okey] = sqrt( errs["nrm2"][ct][okey] / len )
         errs["mae"][ct][okey] = errs["mae"][ct][okey] / len
         errs["nrm1"][ct][okey] = errs["nrm1"][ct][okey] / len
         errs["relrmse"][ct][okey] = errs["rmse"][ct][okey] / errs["nrm2"][ct][okey]
      end
   end

   return errs
end

splittraintest(db::LsqDB, ptrain=0.8) = splittraintest(db.configs, ptrain)

function splittraintest(cfgs::AbstractVector{Dat}, ptrain=0.8)
   cfgtypes = unique(configtype.(cfgs))
   Itrain = Int[]
   Itest = Int[]
   for ct in cfgtypes
      # find all indices of configurations with the current configtype
      Ict = findall(configtype.(cfgs) .== ct)
      # randomly permute the indices
      shuffle!(Ict)
      # first 80% go into train, remaining 20% into test
      ntrain = ceil(Int, length(Ict) * ptrain)
      append!(Itrain, Ict[1:ntrain])
      append!(Itest, Ict[ntrain+1:end])
   end
   return Itrain, Itest
end

end

ext: jl


text:
using Dates: Millisecond
using TimeZones: transition_range

@testset "lt_local / lt_utc" begin
    t = Transition(DateTime(1977, 5, 26, 1), FixedTimeZone("PDT", -8 * 3600, 3600))
    ms = Millisecond(1)

    @testset "lt_local" begin
        local_dt = DateTime(1977, 5, 25, 18)

        @test TimeZones.lt_local(local_dt - ms, t)
        @test !TimeZones.lt_local(local_dt, t)
        @test !TimeZones.lt_local(local_dt + ms, t)

        @test !TimeZones.lt_local(t, local_dt - ms)
        @test !TimeZones.lt_local(t, local_dt)
        @test TimeZones.lt_local(t, local_dt + ms)
    end

    @testset "lt_utc" begin
        utc_dt = DateTime(1977, 5, 26, 1)

        @test TimeZones.lt_utc(utc_dt - ms, t)
        @test !TimeZones.lt_utc(utc_dt, t)
        @test !TimeZones.lt_utc(utc_dt + ms, t)

        @test !TimeZones.lt_utc(t, utc_dt - ms)
        @test !TimeZones.lt_utc(t, utc_dt)
        @test TimeZones.lt_utc(t, utc_dt + ms)
    end
end

@testset "transition_range" begin
    # TODO: Redundancy with test from ZonedDateTime
    @testset "multiple hour transitions" begin
        # Transitions changes that exceed an hour. Results in having two sequential
        # non-existent hour and two sequential ambiguous hours.
        tz = VariableTimeZone("Testing", [
            Transition(DateTime(1800,1,1), FixedTimeZone("TST",0,0)),
            Transition(DateTime(1950,4,1), FixedTimeZone("TDT",0,7200)),
            Transition(DateTime(1950,9,1), FixedTimeZone("TST",0,0)),
        ])

        ### Local ###

        # Initial transition
        @test transition_range(DateTime(1799, 12, 31, 23), tz, Local) == 1:0
        @test transition_range(DateTime(1800, 01, 01, 00), tz, Local) == 1:1

        # A "spring forward" where 2 hours are skipped.
        @test transition_range(DateTime(1950, 03, 31, 23), tz, Local) == 1:1
        @test transition_range(DateTime(1950, 04, 01, 00), tz, Local) == 2:1  # 00:00 TST/TDT
        @test transition_range(DateTime(1950, 04, 01, 01), tz, Local) == 2:1  # 01:00 TST/TDT
        @test transition_range(DateTime(1950, 04, 01, 02), tz, Local) == 2:2

        # A "fall back" where 2 hours are duplicated. Never appears to occur in reality.
        @test transition_range(DateTime(1950, 08, 31, 23), tz, Local) == 2:2
        @test transition_range(DateTime(1950, 09, 01, 00), tz, Local) == 2:3  # 00:00 TDT/TST
        @test transition_range(DateTime(1950, 09, 01, 01), tz, Local) == 2:3  # 01:00 TDT/TST
        @test transition_range(DateTime(1950, 09, 01, 02), tz, Local) == 3:3

        ### UTC ###

        # Initial transition
        @test transition_range(DateTime(1799, 12, 31, 23), tz, UTC) == 1:0
        @test transition_range(DateTime(1800, 01, 01, 00), tz, UTC) == 1:1

        # A "spring forward" where 2 hours are skipped.
        @test transition_range(DateTime(1950, 03, 31, 23), tz, UTC) == 1:1  # 23:00 TST
        @test transition_range(DateTime(1950, 04, 01, 00), tz, UTC) == 2:2  # 02:00 TDT

        # A "fall back" where 2 hours are duplicated. Never appears to occur in reality.
        @test transition_range(DateTime(1950, 08, 31, 21), tz, UTC) == 2:2
        @test transition_range(DateTime(1950, 08, 31, 22), tz, UTC) == 2:2  # 00:00 TDT
        @test transition_range(DateTime(1950, 08, 31, 23), tz, UTC) == 2:2  # 01:00 TDT
        @test transition_range(DateTime(1950, 09, 01, 00), tz, UTC) == 3:3  # 00:00 TST
        @test transition_range(DateTime(1950, 09, 02, 00), tz, UTC) == 3:3  # 01:00 TST
        @test transition_range(DateTime(1950, 09, 03, 00), tz, UTC) == 3:3
    end
end

# Contains both positive and negative UTC offsets and observes daylight saving time.
apia = first(compile("Pacific/Apia", tzdata["australasia"]))

ambiguous_pos = DateTime(2011,4,2,3)
non_existent_pos = DateTime(2011,9,24,3)
ambiguous_neg = DateTime(2012,4,1,3)
non_existent_neg = DateTime(2012,9,30,3)

@test_throws AmbiguousTimeError ZonedDateTime(ambiguous_pos, apia)
@test_throws NonExistentTimeError ZonedDateTime(non_existent_pos, apia)
@test_throws AmbiguousTimeError ZonedDateTime(ambiguous_neg, apia)
@test_throws NonExistentTimeError ZonedDateTime(non_existent_neg, apia)

@test isempty(TimeZones.shift_gap(ambiguous_pos, apia))
@test TimeZones.shift_gap(non_existent_pos, apia) == (
    ZonedDateTime(2011, 9, 24, 2, 59, 59, 999, apia),
    ZonedDateTime(2011, 9, 24, 4, apia),
)
@test isempty(TimeZones.shift_gap(ambiguous_neg, apia))
@test TimeZones.shift_gap(non_existent_neg, apia) == (
    ZonedDateTime(2012, 9, 30, 2, 59, 59, 999, apia),
    ZonedDateTime(2012, 9, 30, 4, apia),
)

# Valid local datetimes close to the non-existent hour should have no boundaries as are
# already valid.
@test isempty(TimeZones.shift_gap(non_existent_pos - Second(1), apia))
@test isempty(TimeZones.shift_gap(non_existent_pos + Hour(1), apia))
@test isempty(TimeZones.shift_gap(non_existent_neg - Second(1), apia))
@test isempty(TimeZones.shift_gap(non_existent_neg + Hour(1), apia))

# Create custom VariableTimeZones to test corner cases
zone = Dict{AbstractString,FixedTimeZone}()
zone["T+0"] = FixedTimeZone("T+0", 0)
zone["T+1"] = FixedTimeZone("T+1", 3600)
zone["T+2"] = FixedTimeZone("T+2", 7200)

# A time zone with a two hour gap
long = VariableTimeZone("Test/LongGap", [
    Transition(DateTime(1800,1,1,0), zone["T+1"])
    Transition(DateTime(1900,1,1,0), zone["T+0"])
    Transition(DateTime(1935,4,1,2), zone["T+2"])
])

# A time zone with an unnecessary transition that typically is hidden to the user
hidden = VariableTimeZone("Test/HiddenTransition", [
    Transition(DateTime(1800,1,1,0), zone["T+1"])
    Transition(DateTime(1900,1,1,0), zone["T+0"])
    Transition(DateTime(1935,4,1,2), zone["T+1"])  # The hidden transition
    Transition(DateTime(1935,4,1,2), zone["T+2"])
])

non_existent_1 = DateTime(1935,4,1,2)
non_existent_2 = DateTime(1935,4,1,3)

# Both "long" and "hidden" are identical for the following tests
for tz in (long, hidden)
    local tz
    boundaries = (
        ZonedDateTime(1935, 4, 1, 1, 59, 59, 999, tz),
        ZonedDateTime(1935, 4, 1, 4, tz),
    )

    @test_throws NonExistentTimeError ZonedDateTime(DateTime(0), tz)
    @test_throws NonExistentTimeError ZonedDateTime(non_existent_1, tz)
    @test_throws NonExistentTimeError ZonedDateTime(non_existent_2, tz)

    # Unhandled datetimes should not be treated as a gap
    @test isempty(TimeZones.shift_gap(DateTime(0), tz))

    @test TimeZones.shift_gap(non_existent_1, tz) == boundaries
    @test TimeZones.shift_gap(non_existent_2, tz) == boundaries
end

# Various DateTimes in Pacific/Apia
valid = DateTime(2013,1,1)
ambiguous = ambiguous_pos
non_existent = non_existent_pos

# first_valid/last_valid with a step
@test TimeZones.first_valid(valid, apia, Hour(1)) == ZonedDateTime(valid, apia)
@test TimeZones.last_valid(valid, apia, Hour(1)) == ZonedDateTime(valid, apia)

@test TimeZones.first_valid(non_existent, apia, Hour(1)) == ZonedDateTime(2011,9,24,4,apia)
@test TimeZones.last_valid(non_existent, apia, Hour(1)) == ZonedDateTime(2011,9,24,2,apia)

@test TimeZones.first_valid(ambiguous, apia, Hour(1)) == ZonedDateTime(ambiguous,apia,1)
@test TimeZones.last_valid(ambiguous, apia, Hour(1)) == ZonedDateTime(ambiguous,apia,2)

# first_valid/last_valid with no step
@test TimeZones.first_valid(valid, apia) == ZonedDateTime(valid, apia)
@test TimeZones.last_valid(valid, apia) == ZonedDateTime(valid, apia)

@test TimeZones.first_valid(non_existent, apia) == ZonedDateTime(2011,9,24,4,apia)
@test TimeZones.last_valid(non_existent, apia) == ZonedDateTime(2011,9,24,2,59,59,999,apia)

@test TimeZones.first_valid(ambiguous, apia) == ZonedDateTime(ambiguous,apia,1)
@test TimeZones.last_valid(ambiguous, apia) == ZonedDateTime(ambiguous,apia,2)


text:
module Constraint

import MutableArithmetics
const MA = MutableArithmetics

import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

import ComplexOptInterface
const COI = ComplexOptInterface

include("split_zero.jl")

end


text:
function _generate_backupsofsimplecontainergenerator_content_1(config::Config)
    simplecontainergenerator_installation_command = config.simplecontainergenerator_installation_command
    lines = String[
        "empty!(Base.DEPOT_PATH)",
        "pushfirst!(Base.DEPOT_PATH, \"/opt/simplecontainergenerator_containers/depot_backup_simplecontainergenerator\")",
        "import Pkg",
        "$(simplecontainergenerator_installation_command)",
    ]
    content = join(lines, "\n") * "\n"
    return content
end

function _generate_backupsofsimplecontainergenerator_content_2(config::Config)
    packagecompiler_installation_command = config.packagecompiler_installation_command
    lines = String[
        "empty!(Base.DEPOT_PATH)",
        "pushfirst!(Base.DEPOT_PATH, \"/opt/simplecontainergenerator_containers/depot_backup_packagecompiler\")",
        "import Pkg",
        "$(packagecompiler_installation_command)",
    ]
    content = join(lines, "\n") * "\n"
    return content
end


text:
abstract type AbstractModel end

transform_param(model, param) = 10. .^ (-6 .+ 8 .* param)


text:
using Luxor

"""
    drawlogo(logo_width::Int)

Draw a `logo_width × logo_width`-sized square image logo for JuliaConstraints.
"""
function drawlogo(logo_width::Int)

    base_width = 1024
    s = logo_width/base_width

    origin() # Set origin to center of image
    scale(s) # Call to scale() must come after origin()

    ref_dist = 320

    ball_centers = [Point(-ref_dist,ref_dist), Point(0,-ref_dist), Point(ref_dist,ref_dist)]
    ball_radius = div(ref_dist,2)
    ball_colors = ["brown3", "forestgreen", "mediumorchid3"]

    link_small_thickness = div(ref_dist,16)*s # line thickness is not scaled http://juliagraphics.github.io/Luxor.jl/stable/transforms/#Scaling-of-lines
    link_large_length = div(ref_dist, 5)
    link_large_width = div(ref_dist, 7)
    link_num = 3 # number of large chains in each direction
    link_color = "royalblue"

    @layer begin # Draw triangle contour for a "small chain links" effect
        setcolor(link_color)
        setline(link_small_thickness)
        poly(ball_centers, :stroke, close=true)
    end

    @layer begin # Draw rounded rectangles for the "large chain links"
        setcolor(link_color)
        angles = (-atan(2), atan(2), 0)
        for i in 1:3
            pt1 = ball_centers[mod1(i, 3)]
            pt2 = ball_centers[mod1(i + 1, 3)]
            for j in 1:link_num
                link_large_shape = squircle(O, link_large_length, link_large_width, rt=0.5, vertices=true)
                link_large = poly(link_large_shape, :transparent)
                polyrotate!(link_large, angles[i])
                polymove!(link_large, O, between(pt1, pt2, j/(link_num+1)))
                setline(link_small_thickness)
                poly(link_large, :fill)
            end
        end
    end

    @layer begin # Draw "chain" balls
        for (ball_center, ball_color) in zip(ball_centers, ball_colors)
            setcolor(ball_color)
            circle(ball_center, ball_radius, :fill)
        end
    end
end

for logo_width in (1024, 512, 128, 32) # large, medium, small and favicon sizes
    Drawing(logo_width, logo_width, joinpath("logo_images", "jc-logo_solid_"*string(logo_width)*"x"*string(logo_width)*".png"))
    drawlogo(logo_width)
    finish()
end

max_stars_repo_name: JuliaConstraints/jc-logo

max_issues_repo_name: JuliaConstraints/jc-logo

max_forks_repo_name: JuliaConstraints/jc-logo


text:
"""
    PrescribedHeatFluxRegion(q::Real,b::Body)

Set up a prescribed heat flux (per unit area) `q` in a region described by the
interior of body `b`.
"""
struct PrescribedHeatFluxRegion{BT}
    q :: Float64
    body :: BT
end

struct PrescribedHeatFlux{GT,MT,BT}
    q :: GT
    mask :: MT
    body :: BT
    cache :: GT
end

PrescribedHeatFlux(p::PrescribedHeatFluxRegion,u::GridData,g::PhysicalGrid) =
            (o = similar(u); fill!(o,1); PrescribedHeatFlux(p.q*o,Mask(p.body,g,u),p.body,o))

(e::PrescribedHeatFlux)() = e.mask(e.cache,e.q)

"""
    PrescribedHeatModelRegion(hc::Real,Tc::Real,b::Body)

Set up a prescribed heat flux model `hc*(Tc - T)` when acting on data `T` in a region
described by the interior of body `b`.
"""
struct PrescribedHeatModelRegion{BT}
    hc :: Float64
    Tc :: Float64
    body :: BT
end

struct PrescribedHeatModel{GT,MT,BT}
    hc :: Float64
    Tc :: GT
    mask :: MT
    body :: BT
    cache1 :: GT
    cache2 :: GT
end

PrescribedHeatModel(p::PrescribedHeatModelRegion,u::GridData,g::PhysicalGrid) =
            (o = similar(u); fill!(o,1); PrescribedHeatModel(p.hc,p.Tc*o,Mask(p.body,g,u),p.body,o,similar(o)))

(c::PrescribedHeatModel{GT})(T::GT) where {GT} = (c.cache1 .= c.hc.*(c.Tc.-T); c.mask(c.cache2,c.cache1))


text:
# Cross Polar Cap Potential (CPCP) calculation at z=2.0 cut plane.
#
# Part 2: Doing surface flux integral and line integral to obtain CPCP.
#
# Hongyang Zhou, hyzhou@umich.edu 03/22/2020

using DelimitedFiles, Glob, PyPlot, Statistics, LinearAlgebra
using NumericalIntegration, ConcaveHull, Random

const Rg = 2634000.                  # [m], radius of Ganymede
const e  = 1.60217662e-19            # [C], electron charge
const Δx = 1/32                      # [Rg],grid resolution
const Δ = Δx*Rg                      # [m]
const Δ²= Δ^2                        # [m^2]
const CoefP = 1e3 * 1e-9 * Rg * 1e-3 # [km/s]*[nT]*[Rg] --> [kV]
const CoefJ = 1e-6 /e * 1e-6* 1e-3  # [μA/m^2]/[#/cc] --> [km/s]
const ϵ = Float32(1e-5) # for comparing floating point values

"Integrate total magnetic flux over the closed surface."
function integrate_surface(data; DoPlot=false, verbose=false)

   bz= @view data[:,7] # [nT]
   x = @view data[:,11]
   y = @view data[:,12]

   # I can also use sort instead.
   bzLow = quantile(bz,0.01)
   bzReasonableRange = quantile(bz,0.95) - quantile(bz,0.05)

   ϕ = 0.0f0
   for z in bz
      bzLow - z > bzReasonableRange && continue # skip suspicious values
      ϕ += z
   end

   if verbose
      b̄z = mean(bz)
      Random.seed!(1234)
      xC = [i + ϵ * rand(Float32) for i in x]
      yC = [i + ϵ * rand(Float32) for i in y]
      v = [[xC[i], yC[i]] for i in 1:length(xC)]
      hull = concave_hull(v, 100) # ConcaveHull
      hull_area = area(hull)

      xB = [p[1] for p in hull.vertices]
      yB = [p[2] for p in hull.vertices]

      if DoPlot
         figure(figsize=(4.5,7.5))
         plt.rc("font", family = "serif", size = 14)
         scatter(xC, yC, s=1, c=bz)
         colorbar()
         plot(xB, yB, "r")

         xlabel(L"x\, [R_G]", size = 14)
         ylabel(L"y\, [R_G]", size = 14)
         title(L"B_z\, [nT]", size = 14)
         xlim(-2.,2.)
         #axis("scaled")
         #tight_layout()
      end

      @info "mean bz = $b̄z"
      @info "area = $hull_area"
      @info "estimate: $(b̄z*hull_area)"
      @info "compute: $(ϕ/32/32)"
   end

   ϕ *= Δ²*1e-9 # [T*m^2]
end

"Not to be confused with the built-in extrema function. Faster."
function minmax(x::AbstractArray)
   a = b = first(x)
   @inbounds @simd for i in eachindex(x)
      if x[i] > b
         b = x[i]
      elseif x[i] < a
         a = x[i]
      end
   end
   a, b
end

"Integrate electric field over the x=xMax line and return the difference."
function integrate_line(data; verbose=false)

   xMax = maximum(@view data[:,11])
   index_ = [i for i in 1:size(data)[1] if isapprox(data[i,11], xMax, atol=ϵ)]

   nP = length(index_)

   ρ = data[index_,1]
   ux= data[index_,2]
   uy= data[index_,3]
   uz= data[index_,4]
   bx= data[index_,5]
   by= data[index_,6]
   bz= data[index_,7]
   jx= data[index_,8]
   jy= data[index_,9]
   jz= data[index_,10]
   y = data[index_,12]

   # Get E field: E = -uxB
   Ex = Vector{Float64}(undef,nP)
   Ey = Vector{Float64}(undef,nP)
   Ez = Vector{Float64}(undef,nP)

   # Calculate electron bulk velocity in Hall MHD
   # Transform the hall velocity into [km/s]
   uex = @. ux - jx/ρ * CoefJ
   uey = @. uy - jy/ρ * CoefJ
   uez = @. uz - jz/ρ * CoefJ

   for i = 1:nP
      Ex[i], Ey[i], Ez[i] =
      -cross([uex[i], uey[i], uez[i]], [bx[i], by[i], bz[i]])
   end

   ϕy = cumul_integrate(y, Ey)
   ϕymin, ϕymax = minmax(ϕy)

   potential_line = (ϕymax - ϕymin) * CoefP # [kV]

   verbose && @info potential_line, y[end]-y[1]

   return potential_line, y[end]-y[1]
end

"""
	get_E_integral(filenames; ftype="hall")

Obtain the upstream electric field integral from surface flux integral and
middle line integration.
"""
function get_E_integral(filenames; ftype="hall")

   nfile = length(filenames)
   ϕ = zeros(Float32, nfile)
   potential_M     = Vector{Float32}(undef, nfile)
   potential_Mhalf = Vector{Float32}(undef, nfile-1)
   yLength         = Vector{Float32}(undef, nfile)

   for i = 1:nfile
      @info i
      fieldline = readdlm(filenames[i], ',', Float32, header=true)
      data = fieldline[1]

      ϕ[i] = integrate_surface(data)
      potential_M[i], yLength[i] = integrate_line(data,verbose=true)
   end

   potential_T = diff(ϕ) .* 1e-3 # [kV]

   # Get the potential drop in the middle line at half times
   for i in 1:nfile-1
      potential_Mhalf[i] = 0.5*(potential_M[i] + potential_M[i+1])
   end

   potential_U = potential_T .+ potential_Mhalf # Upstream potential, [kV]

   if ftype == "hall"
      t = deleteat!(collect(0:1200), [355, 457, 505, 568] .+ 1)
   else # PIC
      t = deleteat!(collect(300:1500), [567, 591, 610, 680] .- 300 .+ 1)
   end

   t_half = zeros(Float32, length(t)-1)

   # Get half times
   for i in 1:length(t)-1
      t_half[i] = 0.5*(t[i] + t[i+1])
   end

   outname = joinpath(dirname(filenames[1]),"CPCP_test.txt")
   isfile(outname) && rm(outname) # Clear the previous existing file

   for i = 1:length(t_half)
      open(outname, "a") do io
         writedlm(io, zip(t_half[i], round(potential_U[i],digits=4),
         round(potential_T[i],digits=4), round(potential_Mhalf[i],digits=4)))
      end
   end

   #@info mean(yLength)
end

function plot_surface(filenames)

   for i = 1:length(filenames)
      fig, ax = subplots(1,1)
      fieldline = readdlm(filenames[i], ',', Float32, header=true)
      data = fieldline[1]
      bz= @view data[:,7] # [nT]
      x = @view data[:,11]
      y = @view data[:,12]
      s = ax.scatter(x, y, s=1, c=bz)
      ax.set_xlim(-2.5,2.0)
      ax.set_ylim(-3.1,3.1)
      axis("scaled")
      fig.colorbar(s)
      title("t=$i")
      savefig("figure/ux"*lpad(i,4,"0")*".png")
      close(fig)
   end

end

"Upstream boundary line integral from closed region point data."
function upstream_boundary_integral(data; DoPlot=false)

   x = @view data[:,11]
   y = @view data[:,12]

   xMax = maximum(x)
   yMin, yMax = minmax(y)

   # Get boundary points through concave hull
   Random.seed!(1234)
   xC = [i + ϵ * rand(Float32) for i in x]
   yC = [i + ϵ * rand(Float32) for i in y]
   v = [[xC[i], yC[i]] for i in 1:length(xC)]
   hull = concave_hull(v, 100) # ConcaveHull
   xH = [p[1] for p in hull.vertices]
   yH = [p[2] for p in hull.vertices]

   # Find boundary indexes
   boundary_index_ = Int32[]
   for iH = axes(xH,1)
      for i = axes(x,1)
         if isapprox(xH[iH], x[i], atol=2ϵ) && isapprox(yH[iH], y[i], atol=2ϵ)
            push!(boundary_index_, i)
            break
         end
      end
   end

   xB = data[boundary_index_,11]
   yB = data[boundary_index_,12]

   line_index_ = @. ((xB == xMax) & (yB !== yMin) & (yB !== yMax))

   up_index_ = .!line_index_
   bc_seq_ = sortperm(yB[up_index_])
   ρ  = data[boundary_index_,1][up_index_][bc_seq_]
   ux = data[boundary_index_,2][up_index_][bc_seq_]
   uy = data[boundary_index_,3][up_index_][bc_seq_]
   uz = data[boundary_index_,4][up_index_][bc_seq_]
   bx = data[boundary_index_,5][up_index_][bc_seq_]
   by = data[boundary_index_,6][up_index_][bc_seq_]
   bz = data[boundary_index_,7][up_index_][bc_seq_]
   jx = data[boundary_index_,8][up_index_][bc_seq_]
   jy = data[boundary_index_,9][up_index_][bc_seq_]
   jz = data[boundary_index_,10][up_index_][bc_seq_]
   xUp= xB[up_index_][bc_seq_]
   yUp= yB[up_index_][bc_seq_]

   if DoPlot
      fig, ax = plt.subplots(1,1,figsize=(6.0,7.5))
      plt.rc("font", family = "serif", size = 14)
      scatter(x, y, s=1, c=data[:,7])
      colorbar()
      plot(xB, yB, "C0", linewidth=3)
      plot([xB[end],xB[1]], [yB[end],yB[1]], "C0", linewidth=3)
      plot(xUp, yUp, "r", linewidth=3)
      #scatter(xB, yB, s=5, c="C0")
      #scatter(xUp, yUp, s=5, c="r")
      plot(xUp[1], yUp[1], "k^",markersize=8)
      plot(xUp[end], yUp[end],"kv",markersize=8)
      ax.minorticks_on()
      xlabel(L"x\, [R_G]", size = 14)
      ylabel(L"y\, [R_G]", size = 14)
      title(L"B_z\, [nT]", size = 14)
      xlim(-2.,2.)
      ax.annotate("U", xy=(-1, 2.03), xytext=(-1.5, 2.5),
         arrowprops=Dict("facecolor"=>"black", "shrink"=>0.02))
      ax.annotate("M", xy=(1.3, 0.1), xytext=(1.7, 0.),
         arrowprops=Dict("facecolor"=>"black", "shrink"=>0.02))
      ax.annotate("A", xy=(1.0, 0.0), xytext=(0., 0.), color="w",size=20)
      #=
      figure(figsize=(2.7,4.8))
      scatter(xB, yB, s=5, c="C0")
      scatter(xUp, yUp, s=5, c="r")
      axis("scaled")
      xlabel(L"x\, [R_G]")
      ylabel(L"y\, [R_G]")
      #title("Boundary points")
      tight_layout()
      =#
   end

   nP = length(xUp)

   # Get E field: E = -uxB
   Ex = Vector{Float64}(undef,nP)
   Ey = Vector{Float64}(undef,nP)
   Ez = Vector{Float64}(undef,nP)

   # Calculate electron bulk velocity in Hall MHD
   # Transform the hall velocity into [km/s]
   uex = @. ux - jx/ρ * CoefJ
   uey = @. uy - jy/ρ * CoefJ
   uez = @. uz - jz/ρ * CoefJ

   for i = 1:nP
      Ex[i], Ey[i], Ez[i] =
      -cross([uex[i], uey[i], uez[i]], [bx[i], by[i], bz[i]])
   end

   ϕup = cumul_integrate(xUp, Ex) .+ cumul_integrate(yUp, Ey)

   line_index_ = (xB .== xMax)
   line_seq_ = sortperm(yB[line_index_],rev=true)
   yL = yB[line_index_][line_seq_]
   xL = xB[line_index_][line_seq_]

   ρL  = data[boundary_index_,1][line_index_][line_seq_]
   uxL = data[boundary_index_,2][line_index_][line_seq_]
   uyL = data[boundary_index_,3][line_index_][line_seq_]
   uzL = data[boundary_index_,4][line_index_][line_seq_]
   bxL = data[boundary_index_,5][line_index_][line_seq_]
   byL = data[boundary_index_,6][line_index_][line_seq_]
   bzL = data[boundary_index_,7][line_index_][line_seq_]
   jxL = data[boundary_index_,8][line_index_][line_seq_]
   jyL = data[boundary_index_,9][line_index_][line_seq_]
   jzL = data[boundary_index_,10][line_index_][line_seq_]

   nL = length(xL)

   # Get E field: E = -uxB
   ExL = Vector{Float64}(undef,nL)
   EyL = Vector{Float64}(undef,nL)
   EzL = Vector{Float64}(undef,nL)

   # Calculate electron bulk velocity in Hall MHD
   # Transform the hall velocity into [km/s]
   uexL = @. uxL - jxL/ρL * CoefJ
   ueyL = @. uyL - jyL/ρL * CoefJ
   uezL = @. uzL - jzL/ρL * CoefJ

   for i = 1:nL
      ExL[i], EyL[i], EzL[i] =
      -cross([uexL[i], ueyL[i], uezL[i]], [bxL[i], byL[i], bzL[i]])
   end

   ϕL = cumul_integrate(yL, EyL)

   if DoPlot
      EIntUp = ϕup.*CoefP
      EIntL = ϕL.*CoefP .+ EIntUp[end]
      figure(figsize=(4.5,3.0))
      scatter(1,EIntUp[1],s=30,c="k",marker="^")
      scatter(nP,EIntUp[end],s=30,c="k",marker="v")
      plot(1:nP,EIntUp,"r")
      plot(nP:nP+1,[EIntUp[end],EIntL[1]],"r")
      plot(nP+1:nP+nL,EIntL, "C0")
      xlabel("Points along curve clockwise")
      ylabel("Electric field integral [kV]")
      tight_layout()
   end

   ϕmin, ϕmax = minmax(ϕup)

   potential_line = (ϕmax - ϕmin) * CoefP # [kV]

   if DoPlot
      @info "Boundary curve potential drop = ", potential_line
      @info "Middle line potential drop = ", (maximum(ϕL) - minimum(ϕL))*CoefP
   end
   return potential_line
end

## Potential drop from surface integration and line integration

#filenames = glob("surface_value*txt","dataSurf_Hall")
#get_E_integral(filenames, ftype="hall")

filenames = glob("surface_value*txt","dataSurf_PIC")
#get_E_integral(filenames, ftype="pic")

## Demonstration for one snapshot

i = 281#952
fieldline = readdlm(filenames[i], ',', Float32, header=true)
data = fieldline[1]
#phi1 = integrate_surface(data,DoPlot=true,verbose=true)
#potential_line1, yL1 = integrate_line(data,verbose=true)
upstream_boundary_integral(data,DoPlot=true)

## Upstream boundary curve integration from closed points data
#=
ftype = "pic"
outname = joinpath("dataBoundary_PIC","CPCP_test.txt")
isfile(outname) && rm(outname) # Clear the previous existing file

if ftype == "hall"
   t = deleteat!(collect(0:1200), [355, 457, 505, 568] .+ 1)
else # PIC
   t = deleteat!(collect(300:1500), [567, 591, 610, 680] .- 300 .+ 1)
end

potential_line = zeros(length(filenames))

for i in 1:length(filenames)
   @info i
   filename = filenames[i]
   fieldline = readdlm(filename, ',', Float32, header=true)
   data = fieldline[1]

   potential_line[i] = upstream_boundary_integral(data,DoPlot=false)

   open(outname, "a") do io
      writedlm(io, [t[i] round(potential_line[i],digits=4)])
   end
end
=#

##
# Next step is time-series plotting.


text:
using PyCall
using SymPy
using REPL

tosub = Dict(
    '1' => '₁',
    '2' => '₂',
    '3' => '₃',
    '4' => '₄',
    '5' => '₅',
    '6' => '₆',
    '7' => '₇',
    '8' => '₈',
    '9' => '₉',
    '0' => '₀',
)

"""
Generate Sym object its names are
x₁,x₂,...,xₙ
and return them as Sym[]
Example
julia> generate_symbol("x",3)
Sym[x₁, x₂, x₃]
"""
function generate_symbol(v::String, n::Int)
    data::Vector{Sym} = []
    for i = 1:n
        seqdigit = [tosub[c] for c in string(i)]
        v_s = "$v" * join(seqdigit)
        eval(Meta.parse("$(v_s)=Sym(\"$(v_s)\")"))
        push!(data, eval(Meta.parse("$(v_s)")))
    end
    return data
end

"""
This func cant use n more than 9
"""
function _generate_symbol(v, n)
    data = Sym[]
    for i = 1:n
        target = "$(v)\\_$i"
        _sub = REPL.REPLCompletions.completions(target, length(target))[1][1]
        sub = REPL.REPLCompletions.completion_text(_sub)
        v_s = v * sub
        eval(Meta.parse("$(v_s)=Sym(\"$(v_s)\")"))
        push!(data, eval(Meta.parse("$(v_s)")))
    end
    return data
end

n = 5
xs = generate_symbol("x", n)
ys = generate_symbol("y", n)

# Uncomment if you can switch randomized numerical data
#xs, ys=rand(n), rand(n)
#

@vars a b

model(x) = a * x + b
E = 1 / 2 * sum(@. (model(xs) - ys)^2)

function bruteforce()
    ∂E∂a, ∂E∂b = diff(E, a), diff(E, b)
    solution = solve([∂E∂a, ∂E∂b], (a, b))
    solution
end

function theoretical()
    denom = n * sum(xs .^ 2) - sum(xs)^2
    a_num = n * sum(xs .* ys) - sum(xs) * sum(ys)
    â = a_num / denom
    b_num = sum(xs .^ 2) * sum(ys) - sum(xs) * sum(xs .* ys)
    b̂ = b_num / denom
    return Dict(a => â, b => b̂)
end

function uselinsolve1()
    ∂E∂a, ∂E∂b = diff(E, a), diff(E, b)
    solution = linsolve((∂E∂a, ∂E∂b), (a, b))
    return solution
end

function uselinsolve2()
    A = [
        diff(E, a, 2) diff(diff(E, a), b)
        diff(diff(E, b), a) diff(E, b, 2)
    ] .|> simplify
    v = -[diff(E, a)(a => 0, b => 0), diff(E, b)(a => 0, b => 0)]
    solution = linsolve((A, v), a, b)
    return solution
end

py"""
def decompose(r):
    return [[eq for eq in eqs] for eqs in r.args]
"""

r0 = theoretical()
r1 = bruteforce()
r2 = py"decompose"(uselinsolve1())
r3 = py"decompose"(uselinsolve2())

@show r0[a]
@show r1[a]
@show r2[1] # a
@show r3[1] # a

@show r0[b]
@show r1[b]
@show r2[2] # b
@show r3[2] # b

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
include("backwards_activation_functions.jl")

"""
Partial derivates of the components of linear forward function using the linear output (δZ) and caches of these components (cache).
"""
function linear_backward(δZ, cache)
    # unpack cache
    A_prev, W, b = cache
    m = size(A_prev, 2)

    # partial derivates of each component
    δW = δZ * (A_prev') / m
    δb = sum(δZ, dims = 2) / m
    δA_prev = (W') * δZ

    @assert (size(δA_prev)[1] == size(A_prev)[1])
    @assert (size(δW) == size(W))
    @assert (size(δb)[1] == size(b)[1])

    return δW, δb, δA_prev
end # function

"""
Unpack the linear activated caches (cache) and compute their derivates from the applied activation function.
"""
function linear_activation_backwards(δA, cache; activation_function="relu")
    @assert activation_function ∈ ("sigmoid", "relu", "softmax", "tanh", "swish")

    linear_cache, cache_activation = cache

    if (activation_function == "relu")

        δZ = relu_backwards(δA, cache_activation)

    elseif (activation_function == "sigmoid")

        δZ = sigmoid_backwards(δA, cache_activation)

    elseif (activation_function ==  "tanh")

        δZ = tanh_backwards(δA, cache_activation)

    elseif (activation_function == "softmax")

        δZ = softmax_backwards(δA, cache_activation)

    elseif (activation_function == "swish")

        δZ = swish_backwards(δA, cache_activation)

    end # if

    δW, δb, δA_prev = linear_backward(δZ, linear_cache)

    return δW, δb, δA_prev
end # function

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
w11 = [0.1, 6.8, -3.5, 2.0, 4.1, 3.1, -0.8, 0.9, 5.0, 3.9];
w12 = [1.1, 7.1, -4.1, 2.7, 2.8, 5.0, -1.3, 1.2, 6.4, 4.0];
w21 = [7.1, -1.4, 4.5, 6.3, 4.2, 1.4, 2.4, 2.5, 8.4, 4.1];
w22 = [4.2, -4.3, 0.0, 1.6, 1.9, -3.2, -4.0, -6.1, 3.7, -2.2];
w31 = [-3.0, 0.5, 2.9, -0.1, -4.0, -1.3, -3.4, -4.1, -5.1, 1.9];
w32 = [-2.9, 8.7, 2.1, 5.2, 2.2, 3.7, 6.2, 3.4, 1.6, 5.1];
w41 = [2.0, 8.9, 4.2, 8.5, 6.7, 0.5, 5.3, 8.7, 7.1, 8.0] * -1;
w42 = [8.4, -0.2, 7.7, 3.2, 4.0, 9.2, 6.7, 6.4, 9.7, 6.3] * -1;
const SAMPLE = [w11 w12 w21 w22 w31 w32 w41 w42];

function normalize_data(positive, negative)
    " normalization adds additional one column to data "
    normalize = data -> [data ones(size(data)[1])]

    positive = normalize(positive)
    negative = normalize(negative) * -1
    assert(size(positive) == size(negative))
    data = transpose([positive; negative])

    d, n = size(data)
    weight = zeros(d)

    data, weight
end

function pseudo_inverse(positive, negative)
    data, weight = normalize_data(positive, negative)
    d, n = size(data)

    weight = pinv(data * transpose(data)) * data * ones(n)
end

" least mean square method to computer weight vector "
function lms(positive, negative; eta=1, margin=1, theta=0.01)
    gradient(y) = (1 - dot(y, weight)) / dot(y, y) * y
    loss(weight) = sum([1 - dot(data[:, i], weight) for i = 1:n])

    data, weight = normalize_data(positive, negative)
    d, n = size(data)
    weight = [0.08276694629366121,0.12681539660924512,0.4783531495651045]

    function compute_delta(data, step, weight, eta, margin)
        i = (step % n + 1)
        for i = 1:n
            delta = eta / i * gradient(data[:, i])
            if dot(delta, delta) >= theta^2
                return delta
            end
        end
        nothing
    end

    step, delta = 1, compute_delta(data, 1, weight, eta, margin)
    while delta != nothing
        println("$step,$(loss(weight))")

        weight += delta
        step += 1
        delta = compute_delta(data, step, weight, eta, margin)
    end

    weight, step
end

function single_sample_relax(positive, negative; eta=1, margin=1, epsilon=0.00001)
    data, weight = normalize_data(positive, negative)
    d, n = size(data)

    function find_error(data, weight, eta, margin)
        for i = 1:n
            if dot(weight, data[:, i]) <= margin - epsilon
                return data[:, i]
            end
        end
        nothing
    end

    step, y = 1, find_error(data, weight, eta, margin)
    while y != nothing
        gradient = (margin - dot(weight, y)) / dot(y, y) * y
        weight += eta * gradient

        step += 1
        y = find_error(data, weight, eta, margin)
    end

    weight, step
end

function relax_procedure(positive, negative; eta=1, margin=1, epsilon=0.00001)
    data, weight = normalize_data(positive, negative)
    d, n = size(data)

    println("data = $data")

    function compute_error(data, weight, eta, margin)
        error_samples = [y for y in filter([data[:, i] for i = 1:n]) do y
            dot(y, weight) <= margin - epsilon
        end]
    end

    step, error_samples = 1, compute_error(data, weight, eta, margin)
    while length(error_samples) > 0
        println("== $(length(error_samples)) ==========================================")
        println("error_samples = $error_samples")
        println("weight = $weight")

        gradient = sum(error_samples) do y
            #println(y);
            #println(margin - dot(weight, y) / dot(y, y))
            (margin - dot(weight, y)) / dot(y, y) * y
        end

        weight += eta * gradient
        println("new weight $weight\t<= $gradient")

        step += 1
        error_samples = compute_error(data, weight, eta, margin)
    end

    weight, step
end

function batch_perception(positive, negative)
    data, weight = normalize_data(positive, negative)
    d, n = size(data)

    function compute_error(data, weight)
        error_samples = [y for y in filter([data[:, i] for i = 1:n]) do y
            dot(weight, y) <= 0
        end]
    end

    step, error_samples = 1, compute_error(data, weight)

    while length(error_samples) > 0
        gradient = sum(error_samples)
        weight += gradient
        #println("new weight $weight\t<= $gradient")

        step += 1
        error_samples = compute_error(data, weight)
    end

    weight, step
end

"""
Write a program to implement the “batch perception” algorithm
(see page 44 or 45 in PPT).

(a). Starting with a = 0, apply your program to the training data from 1 and 2.
Note that the number of iterations required for convergence.

(b). Apply your program to the training data from 3 and 2.
Again, note that the number of iterations required for convergence.

(c). Explain the difference between the iterations required in the above two cases.

写一个基于松驰准则的感知函数学习的代码，并采用 w1 和 w2 的数据进行训练。对程序运算过程和结果进行分析。

写一个最小平方准则的线性判别函数学习的代码，并分别采用 w1 和 w2、w2 和 w3、w3 和 w4 的数据进行训练。
"""
function hw2_1()
    weight, step = batch_perception(SAMPLE[:, 1:2], SAMPLE[:, 3:4])
    println("batch perception w1 and w2: iteration steps: $step, final weights: $weight")
    weight, step = batch_perception(SAMPLE[:, 3:4], SAMPLE[:, 5:6])
    println("batch perception w2 and w3: iteration steps: $step, final weights: $weight")

    weight, step = single_sample_relax(SAMPLE[:, 1:2], SAMPLE[:, 3:4])
    println("single sample relax w1 and w2: iteration steps: $step, final weights: $weight")

    weight = pseudo_inverse(SAMPLE[:, 1:2], SAMPLE[:, 3:4])
    println("pseudo_inverse MSE w1 and w2: final weights: $weight")
    weight = pseudo_inverse(SAMPLE[:, 3:4], SAMPLE[:, 5:6])
    println("pseudo_inverse MSE w2 and w3: final weights: $weight")
    weight = pseudo_inverse(SAMPLE[:, 5:6], SAMPLE[:, 7:8])
    println("pseudo_inverse MSE w3 and w4: final weights: $weight")

end

if !isinteractive()
    hw2_1()
end

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
module Scruff
include("../MultiInterface.jl/src/MultiInterface.jl")  ## Scruff.MultiInterface module

include("utils.jl")            ## Scruff.Utils module
include("core.jl")             ## Scruff module
include("operators.jl")        ## Scruff.Operators module
include("sfuncs.jl")           ## Scruff.SFuncs module
include("models.jl")           ## Scruff.Models module
include("runtime.jl")          ## Scruff module
include("runtime_utils.jl")    ## Scruff.RTUtils module
include("algorithms.jl")       ## Scruff.Algorithms module
include("logger.jl")           ## Scruff module

end


text:
abstract type Dof <: Map end

# """
#     abstract type Dof <: Map

# Abstract type representing a degree of freedom (DOF), a basis of DOFs, and related objects.
# These different cases are distinguished by the return type obtained when evaluating the `Dof`
# object on a `Field` object. See function [`evaluate_dof!`](@ref) for more details.

# The following functions needs to be overloaded

# - [`dof_cache`](@ref)
# - [`evaluate_dof!`](@ref)

# The following functions can be overloaded optionally

# - [`dof_return_type`](@ref)

# The interface is tested with

# - [`test_dof`](@ref)

# In most of the cases it is not strictly needed that types that implement this interface
# inherit from `Dof`. However, we recommend to inherit from `Dof`, when possible.

# """
# abstract type Dof <: Map end

# """
#     return_cache(dof,field)

# Returns the cache needed to call `evaluate_dof!(cache,dof,field)`
# """
# function return_cache(dof::Dof,field)
#   @abstractmethod
# end

# """
#     evaluate_dof!(cache,dof,field)

# Evaluates the dof `dof` with the field `field`. It can return either an scalar value or
# an array of scalar values depending the case. The `cache` object is computed with function
# [`dof_cache`](@ref).

# When a mathematical dof is evaluated on a physical field, a scalar number is returned. If either
# the `Dof` object is a basis of DOFs, or the `Field` object is a basis of fields,
# or both objects are bases, then the returned object is an array of scalar numbers. The first
# dimensions in the resulting array are for the `Dof` object and the last ones for the `Field`
# object. E.g, a basis of `nd` DOFs evaluated at physical field returns a vector of `nd` entries.
# A basis of `nd` DOFs evaluated at a basis of `nf` fields returns a matrix of size `(nd,nf)`.
# """
# function evaluate!(cache,dof::Dof,field)
#   @abstractmethod
# end

# """
#     dof_return_type(dof,field)

# Returns the type for the value obtained with evaluating `dof` with `field`.

# It defaults to

#     typeof(evaluate_dof(dof,field))
# """
# function return_type(dof::Dof,field)
#   typeof(evaluate(dof,field))
# end

# Testers

"""
    test_dof(dof,field,v;cmp::Function=(==))

Test that the `Dof` interface is properly implemented
for object `dof`. It also checks if the object `dof`
when evaluated at the field `field` returns the same
value as `v`. Comparison is made with the `comp` function.
"""
function test_dof(dof::Dof,field,v;cmp::Function=(==))
  _test_dof(dof,field,v,cmp)
end

function test_dof_array(dof::AbstractArray{<:Dof},field,v;cmp::Function=(==))
  _test_dof(dof,field,v,cmp)
end

function _test_dof(dof,field,v,cmp)
  if isa(dof,Dof)
    test_map(v,dof,field;cmp=cmp)
  end
  r = evaluate(dof,field)
  @test cmp(r,v)
  @test typeof(r) == return_type(dof,field)
end

#struct DofEval <: Map end
#
#function return_cache(k::DofEval,dof,field)
#  return_cache(dof,field)
#end
#
#@inline function evaluate!(cache,k::DofEval,dof,field)
#  evaluate!(cache,dof,field)
#end
#
#function return_type(k::DofEval,dof,field)
#  return_type(dof,field)
#end


text:
function plot_flx_corrs(simid, Ds, ϵs, cg, is_glclim_stst)

    MEmode = "ME_AVZ_EXPECTED_AVUG_BOUNDED"

    # dat pool
    datpool = Dict()

    _push!(d, k, v) = push!(get!(d, k, typeof(v)[]), v)

    # Collect
    for D in Ds, (ϵi, ϵ) in enumerate(ϵs)
        simparams = (;D, ϵ, cg)

        ## ------------------------------------------------------
        # Dynamic
        simparams = (;D, ϵ, cg)
        status = get_status(simid, simparams)
        if (status != STST_SIM_STATUS)
            # @warn("Wrong status ", (D, ϵ, cg), status)
            continue
        end

        datfile = Dyn.simdat_file(simid, simparams)
        if !isfile(datfile)
            # @warn("File missing ", (D, ϵ, cg), datfile)
            continue
        end

        if !is_glclim_stst[simparams]
            @warn("Not glc-lim STST ", status)
            continue
        end

        S = ldat(datfile)
        PX = S.P
        V = S.V
        Vz = Vi(V, :z)
        Vug = Vi(V, :ug)
        Vuo = Vi(V, :uo)

        # Dyn dat
        z_avPX = sum(PX .* Vz)
        ug_avPX = sum(PX .* Vug)
        uo_avPX = sum(PX .* Vuo)
        cgD_X = S.cgD_X

        ## ------------------------------------------------------
        # LP
        net0 = S.V.net

        # fmm & fmM
        net = deepcopy(net0)
        Dyn.fix!(net, :z, D)
        Dyn.ub!(net, :ug, cgD_X)
        Lug, _ = Dyn.fva(net, :ug)
        Dyn.fix!(net, :ug, Lug)
        Dyn.reducebox!(net)
        Luo, Uuo = Dyn.fva(net, :uo)

        z_fba_fmm, ug_fba_fmm, uo_fba_fmm = D, Lug, Luo
        z_fba_fmM, ug_fba_fmM, uo_fba_fmM = D, Lug, Uuo

        # fMm & fMM
        net = deepcopy(net0)
        Dyn.fix!(net, :z, D)
        Dyn.ub!(net, :ug, cgD_X)
        _, Uug = Dyn.fva(net, :ug)
        Dyn.fix!(net, :ug, Uug)
        Dyn.reducebox!(net)
        Luo, Uuo = Dyn.fva(net, :uo)

        z_fba_fMm, ug_fba_fMm, uo_fba_fMm = D, Uug, Luo
        z_fba_fMM, ug_fba_fMM, uo_fba_fMM = D, Uug, Uuo

        # M00
        net = deepcopy(net0)
        Dyn.ub!(net, :ug, cgD_X)
        Dyn.reducebox!(net)
        _, Uz = Dyn.fva(net, :z)
        Dyn.fix!(net, :z, Uz)
        Dyn.reducebox!(net)
        Lug, Uug = Dyn.fva(net, :ug)
        if !isapprox(Lug, Uug)
            @warn("Lug != Uug", (D, ϵ, cg), Lug, Uug)
            continue
        end
        Dyn.fix!(net, :ug, Lug)
        Dyn.reducebox!(net)
        Luo, Uuo = Dyn.fva(net, :uo)
        !isapprox(Luo, Uuo) && @warn("Luo != Uuo", (D, ϵ, cg), Luo, Uuo)

        z_fba_M00, ug_fba_M00, uo_fba_M00 = Uz, Lug, Luo

        ## ------------------------------------------------------
        # MaxEnt
        status = Dyn.get_status(simid, simparams)
        if (status != STST_SIM_STATUS)
            @warn("Wrong status ", (D, ϵ, cg), status)
            continue
        end

        datfile = Dyn.maxent_file(simid, MEmode, simparams)
        if !isfile(datfile)
            @warn("File missing ", (D, ϵ, cg), datfile)
            continue
        end
        ME_Dat = ldat(datfile)
        PME = ME_Dat.PME

        # me dat
        z_avPME = sum(PME .* Vz)
        ug_avPME = sum(PME .* Vug)
        uo_avPME = sum(PME .* Vuo)

        if (abs(z_avPME - D) / D) > 0.05
            @warn("D and z not converged", status, z_avPME, D)
            continue
        end

        # collect
        for (plot_id, dyn_v, model_v) in [
            (:fba_fmm_z, z_avPX, z_fba_fmm),
            (:fba_fmm_ug, ug_avPX, ug_fba_fmm),
            (:fba_fmm_uo, uo_avPX, uo_fba_fmm),

            (:fba_fmM_z, z_avPX, z_fba_fmM),
            (:fba_fmM_ug, ug_avPX, ug_fba_fmM),
            (:fba_fmM_uo, uo_avPX, uo_fba_fmM),

            (:fba_fMM_z, z_avPX, z_fba_fMM),
            (:fba_fMM_ug, ug_avPX, ug_fba_fMM),
            (:fba_fMM_uo, uo_avPX, uo_fba_fMM),

            (:fba_fMm_z, z_avPX, z_fba_fMm),
            (:fba_fMm_ug, ug_avPX, ug_fba_fMm),
            (:fba_fMm_uo, uo_avPX, uo_fba_fMm),

            (:fba_M00_z, z_avPX, z_fba_M00),
            (:fba_M00_ug, ug_avPX, ug_fba_M00),
            (:fba_M00_uo, uo_avPX, uo_fba_M00),

        ]
            plt_dat = get!(datpool, plot_id, Dict())
            _push!(plt_dat, :dyn, dyn_v)
            _push!(plt_dat, :model, model_v)
            _push!(plt_dat, :ϵi, ϵi)
            _push!(plt_dat, :ϵ, ϵ)
        end

        for (plot_id, dyn_v, model_v) in [
            (:me_z, z_avPX, z_avPME)
            (:me_ug, ug_avPX, ug_avPME)
            (:me_uo, uo_avPX, uo_avPME)
        ]
            plt_dat = get!(datpool, plot_id, Dict())
            _push!(plt_dat, :dyn, dyn_v)
            _push!(plt_dat, :model, model_v)
            _push!(plt_dat, :ϵi, ϵi)
            _push!(plt_dat, :ϵ, ϵ)
        end

    end

    # # Test
    # datpool = lcache(ChE, :TEMP_CACHE) do
    #     datpool
    # end

    # sort
    _sort(f::Function, arr) = sort(arr; by = f)
    sdatpool = _sort(datpool) do elm
        model_label = split(string(elm), "_")
        replace!(model_label, "z" => "1z", "ug" => "2ug", "uo" => "3uo")
        join(model_label)
    end

    # plot
    ps = Plots.Plot[]
    for (plot_id, dat) in sdatpool
        @info("At", plot_id)

        model_label = split(string(plot_id), "_")
        replace!(model_label, "ug" => "(u_g)", "uo" => "(u_o)", "z" => "(z)", "fba" => "FBA", "me" => "ME")
        flux_label = last(model_label)
        model_label = join(model_label, " ")
        model_label = _textbf(model_label)
        dyn_label = _textbf("Dynamic", " ", flux_label)

        p = plot(;
            xlabel = dyn_label,
            ylabel = model_label
        )
        push!(ps, p)

        dyn_vs = dat[:dyn]
        model_vs = dat[:model]
        ϵs_ = dat[:ϵ]
        # ϵis_ = dat[:ϵi]

        max_ϵs_ = maximum(ϵs)
        ms0 = 6

        scatter!(p, dyn_vs, model_vs;
            label = "",
            color = :black,
            alpha = 0.65,
            m = :circle,
            ms = ms0 .+ (2 * ms0) .* (ϵs_ ./ max_ϵs_)
        )

        all_ = sort!([dyn_vs; model_vs])
        lims = (last(all_) - first(all_))
        lims = (first(all_) - (lims * 0.1), last(all_) + (lims * 0.1))
        plot!(p, all_, all_;
            label = "", ls = :dash, c = :black, lw = 3, xlim = lims, ylim = lims
        )

        # annotation
        # if plot_id in [:fba_M00_ug, :fba_fmm_ug, :fba_fmM_ug, :fba_fMm_ug, :fba_fMM_ug, :me_ug]
            # plot!(p; leftmargin = 18mm)
            # m_, M_ = sort(all_)[[1, end]]
            # text_ = split(string(plot_id), "_")
            # replace!(text_, "ug" => "(u_g)", "uo" => "(u_o)", "z" => "(z)", "fba" => "FBA", "me" => "ME")
            # text_ = _textbf(text_...)
            # annotate!(p, m_ - (M_ - m_) * 0.25, m_ + (M_ - m_) * 0.5, text(text_, 25, :black, rotation = 90))
        # end

        # # title
        # if plot_id == :fba_M00_ug
        #     plot!(p; title = _textbf("u_g"))
        # elseif plot_id == :fba_M00_uo
        #     plot!(p; title = _textbf("u_o"))
        # elseif plot_id == :fba_M00_z
        #     plot!(p; title = _textbf("z"))
        # else
        #     plot!(p; title = _textbf(""))
        # end
        # plot!(p; title = "$plot_id") # Test

    end

    return ps

    # return [lp_Mug_pug, lp_mug_pug, me_pug]

end


text:
using CMBLensing
using CMBLensing: basis, Basis, BasisTuple, @SVector, @SMatrix, RK4Solver

##

using Test
using SparseArrays
using LinearAlgebra
using Zygote

##

@testset "CMBLensing" begin

##

# basic printing sanity checks, which were super annoying to get right...
# see also: https://discourse.julialang.org/t/dispatching-on-the-result-of-unwrap-unionall-seems-weird/25677

@testset "Printing" begin

    # concrete types:
    for f in [FlatMap(rand(4,4)), FlatQUMap(rand(4,4),rand(4,4))]
        @test occursin("pixels",sprint(show, MIME("text/plain"), f))
        @test occursin("pixels",sprint(show, MIME("text/plain"), [f,f]))
    end

    for m in ((), (MIME("text/plain"),))
        # unionall types: (the presence of "where" indicates printing correctly
        # forwarded to the default behavior)
        @test occursin("where",sprint(show, m..., FieldTuple{<:Any,<:NamedTuple{(:Q,:U)}}))
        @test occursin("where",sprint(show, m..., FlatMap{<:Any,<:Any,<:Matrix{Real}}))
        @test occursin("where",sprint(show, m..., FΦTuple))
        @test occursin("where",sprint(show, m..., FlatQUMap))
    end

end
##

@testset "FieldTuples" begin

    f = FlatMap(rand(4,4))

    @testset "Constructors" begin

        # Enumerate the combinations of
        # * 1) no basis or names specified, 2) only basis specified, 3) only names specified, 4) basis and names specified
        # * 1) args/kwargs form 2) single Tuple/NamedTuple argument form

        for F in [
            FieldTuple,
            FieldTuple{QUMap},
            FieldTuple{QUMap,<:NamedTuple{(:Q,:U)}},
            FieldTuple{<:Basis,<:NamedTuple{(:Q,:U)}}
            ]

            @testset "F :: $F" begin
                @test (@inferred F(;Q=f, U=f)) isa F
                @test (@inferred F((Q=f, U=f))) isa F
                @test (@inferred F(f,f)) isa F
                @test (@inferred F((f,f))) isa F
            end

        end

    end

    @testset "Basis conversions" begin

        # basis conversions
        for f_basistuple in [FieldTuple(f, f), FieldTuple(A=f, B=f)] # named and unnamed
            @test basis(@inferred    Fourier(f_basistuple)) <: BasisTuple{Tuple{Fourier,Fourier}}
            @test basis(@inferred        Map(f_basistuple)) <: BasisTuple{Tuple{Map,Map}}
            @test basis(@inferred DerivBasis(f_basistuple)) <: BasisTuple{Tuple{Fourier,Fourier}}
            @test basis(@inferred BasisTuple{Tuple{Fourier,Fourier}}(f_basistuple)) <: BasisTuple{Tuple{Fourier,Fourier}}
        end

        f_concretebasis = FieldTuple{QUMap, <:NamedTuple{(:Q,:U)}}(f,f)
        @test basis(@inferred    Fourier(f_concretebasis)) <: QUFourier
        @test basis(@inferred        Map(f_concretebasis)) <: QUMap
        @test basis(@inferred DerivBasis(f_concretebasis)) <: QUFourier

    end

end

##

@testset "Flat Constructors" begin

    N = 2
    θpix = 3
    kwargs = (Nside=N, θpix=θpix)
    P = Flat(;kwargs...)
    Ix = rand(N,N)
    Il = rand(N÷2+1,N) + im*rand(N÷2+1,N)

    for (F,args) in [
            (FlatMap,        (Ix,)),
            (FlatFourier,    (Il,)),
            (FlatQUMap,      (Ix,Ix)),
            (FlatQUFourier,  (Il,Il)),
            (FlatEBMap,      (Ix,Ix)),
            (FlatEBFourier,  (Il,Il)),
            (FlatIQUMap,     (Ix,Ix,Ix)),
            (FlatIQUFourier, (Il,Il,Il)),
            (FlatIEBMap,     (Ix,Ix,Ix)),
            (FlatIEBFourier, (Il,Il,Il))
        ]
        @testset "f::$F" begin
            @test F(args...; kwargs...) isa F{P}
            @test (@inferred F{P}(args...)) isa F{P}
            @test (@inferred broadcast(real, (F{P}(args...)))) isa F{P}
            if eltype(args[1]) <: Complex
                @test (@inferred F{P}(map(real,args)...)) isa F{P}
            end
            @test real(eltype(@inferred F{P,Float32}(args...))) == Float32
        end
    end

end

##

@testset "Algebra" begin

    fs = ((B0,f0),(B2,f2),(Bt,ft)) = [
        (Fourier,   FlatMap(rand(4,4))),
        (EBFourier, FlatQUMap(rand(4,4),rand(4,4))), # named FieldTuple
        (Fourier,   FieldTuple(FlatMap(rand(4,4)),FlatMap(rand(4,4)))), # unnamed FieldTuple
        (BasisTuple{Tuple{Fourier,EBFourier}}, FlatIQUMap(rand(4,4),rand(4,4),rand(4,4))), # named nested FieldTuple,
        (BasisTuple{Tuple{Fourier,EBFourier}}, FieldTuple(FlatMap(rand(4,4)),FlatQUMap(rand(4,4),rand(4,4)))) # unnamed nested FieldTuple
    ]

    for (B,f) in fs

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

            local Ðf, Ðv, g, H

            @test similar(f) isa typeof(f)
            @test zero(f) isa typeof(f)
            @test similar(f,Float32) isa Field
            @test eltype(similar(f,Float32)) == Float32

            # broadcasting
            @test (@inferred f .+ f) isa typeof(f)
            @test (@inferred f .+ Float32.(f)) isa typeof(f)

            # promotion
            @test (@inferred f + B(f)) isa typeof(f)
            @test (@inferred f + B(Float32.(f))) isa typeof(f)

            # gradients
            @test (Ðf = @inferred ∇[1]*f) isa Field
            @test (∇[1]'*f ≈ -∇[1]*f)
            @test (-∇[1]'*f ≈ ∇[1]*f)
            @test (@inferred mul!(Ðf,∇[1],Ð(f))) isa Field
            @test (Ðv = @inferred ∇*f) isa FieldVector
            @test (@inferred mul!(Ðv,∇,Ð(f))) isa FieldVector
            @test ((g,H) = map(Ł, (@inferred gradhess(f)))) isa NamedTuple{<:Any, <:Tuple{FieldVector, FieldMatrix}}

            # Diagonal broadcasting
            @test (@inferred Diagonal(f) .* Diagonal(f) .* Diagonal(f)) isa typeof(Diagonal(f))

            # inverses
            @test (@inferred pinv(Diagonal(f))) isa Diagonal{<:Any,<:typeof(f)}
            @test_throws SingularException inv(Diagonal(0*f))

            # Field dot products
            D = Diagonal(f)
            @test (@inferred f' * f) isa Real
            @test (@inferred f' * B(f)) isa Real
            @test (@inferred f' * D * f) isa Real
            @test (@inferred sum(f)) isa Number

            if f isa FlatS0
                # FieldVector dot product
                @test (@inferred Diagonal.(g)' * g) isa typeof(g[1])
                @test (@inferred mul!(similar(g[1]), Diagonal.(g)', g)) isa typeof(g[1])

                # FieldMatrix-FieldVector product
                @test (@inferred Diagonal.(H) * g) isa FieldVector
                @test (@inferred Diagonal.(H) * Diagonal.(g)) isa FieldOrOpVector
                @test (@inferred mul!(Diagonal.(similar.(g)), Diagonal.(H), Diagonal.(g))) isa FieldOrOpVector

            end

            # Explicit vs. lazy DiagOp algebra
            @test (Diagonal(Ð(f)) + Diagonal(Ð(f))) isa DiagOp{<:Field{basis(Ð(f))}}
            @test (Diagonal(Ł(f)) + Diagonal(Ð(f))) isa LazyBinaryOp
            @test (Diagonal(Ł(f)) + Diagonal(Ł(f))) isa DiagOp{<:Field{basis(Ł(f))}}

            # tuple adjoints
            v = similar.(@SVector[f0,f0])
            @test (@inferred mul!(f0, tuple_adjoint(f), f)) isa Field{<:Any,S0}
            @test (@inferred mul!(v, tuple_adjoint(f), @SVector[f,f])) isa FieldVector{<:Field{<:Any,S0}}

        end

    end

    # mixed-spin
    @test (@inferred f0 .* f2) isa typeof(f2)
    @test (@inferred f0 .* ft) isa typeof(ft)

    # matrix type promotion
    @test_broken (@inferred FlatMap(rand(Float32,2,2)) + FlatMap(spzeros(Float64,2,2))) isa FlatMap{<:Any,Float64,<:Matrix}

end

##

@testset "FlatS02" begin

    ΣTT, ΣTE, ΣEE, ΣBB = [Diagonal(Fourier(FlatMap(rand(3,3)))) for i=1:4]
    L = FlatIEBCov(@SMatrix([ΣTT ΣTE; ΣTE ΣEE]), ΣBB)
    f = IEBFourier(FlatIEBMap(rand(3,3),rand(3,3),rand(3,3)))

    @test (sqrt(L) * @inferred(@inferred(sqrt(L)) * f)) ≈ (L * f)
    @test (L * @inferred(@inferred(pinv(L)) * f)) ≈ f
    @test @inferred(L * L) isa FlatIEBCov
    @test @inferred(L + L) isa FlatIEBCov
    @test L * Diagonal(f) isa FlatIEBCov
    @test Diagonal(f) * L isa FlatIEBCov
    @test_broken @inferred L * Diagonal(f)
    @test @inferred(diag(L)) isa DiagOp{<:FlatIEBFourier}
    @test @inferred(L + I) isa FlatIEBCov
    @test @inferred(2 * L) isa FlatIEBCov
    @test @inferred(similar(L)) isa FlatIEBCov
    @test (L .= 2L) isa FlatIEBCov

end

##

@testset "Gradients" begin

    @test (@inferred ∇[1] * FlatMap(rand(3,3), ∂mode=fourier∂)) isa FlatFourier
    @test (@inferred ∇[1] * FlatQUMap(rand(3,3), rand(3,3), ∂mode=fourier∂)) isa FlatQUFourier
    @test (@inferred ∇[1] * FlatIQUMap(rand(3,3), rand(3,3), rand(3,3), ∂mode=fourier∂)) isa FlatIQUFourier

    @test (@inferred ∇[1] * Fourier(FlatMap(rand(3,3), ∂mode=map∂))) isa FlatMap
    @test (@inferred ∇[1] * QUFourier(FlatQUMap(rand(3,3), rand(3,3), ∂mode=map∂))) isa FlatQUMap
    @test (@inferred ∇[1] * BasisTuple{Tuple{Fourier,QUFourier}}(FlatIQUMap(rand(3,3), rand(3,3), rand(3,3), ∂mode=map∂))) isa FlatIQUMap

end

##

@testset "Misc" begin

    f = FlatMap(rand(4,4))

    @test                  @inferred(MidPass(100,200) .* Diagonal(Fourier(f))) isa Diagonal
    @test_throws Exception           MidPass(100,200) .* Diagonal(        f)

end

##

@testset "ParamDependentOp" begin

    D = Diagonal(FlatMap(rand(4,4)))
    mem = similar(D)

    @test_throws ArgumentError ParamDependentOp((;x=1, y=1)->x*y*D)(mem) # passing memory to non-inplace op
    @test_throws ArgumentError ParamDependentOp((mem;x=1, y=1)->mem.=x*y*D,similar(D))(1) # passing wrong-type memory

    @test ParamDependentOp((;x=1, y=1)->x*y*D)() ≈ D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)(z=2) ≈ D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)(x=2) ≈ 2D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)((x=2,y=2)) ≈ 4D # tuple calling form
    @test ParamDependentOp((mem;x=1, y=1)->mem.=x*y*D,similar(D))(D) ≈ D # inplace
end
##

# make sure we can take type-stable gradients of scalar functions of our Fields
# (like the posterior)
@testset "Autodiff" begin

    f = FlatMap(rand(2,2))
    D = Diagonal(f)

    # the very basics
    grad1() = Zygote.gradient(function (θ)
        g = θ * f
        dot(g,g)
    end, 1)[1]
    # this one *is* inferred sometimes, can't figure out what combination of
    # versions though, so leaving it without @inferred for now
    @test (grad1()) ≈ 2*norm(f,2)^2

    grad2() = Zygote.gradient(function (θ)
        g = (θ * D * f)
        dot(g,g)
    end, 1)[1]
    @test_broken                  grad2()  ≈ 2*norm(f.^2,2)^2

    # derivatives through ParamDependentOps
    Dr = ParamDependentOp((;r=1)-> r * D)
    grad3() = Zygote.gradient(function (r)
        g = (Dr(r=r) * f)
        dot(g,g)
    end,1)[1]
    @test_broken           grad3()  ≈ 2*norm(f.^2,2)^2
    @test_broken @inferred(grad3()) ≈ 2*norm(f.^2,2)^2 # would be nice to get this inferred

    @test_broken @inferred(Zygote.gradient(r->logdet(Dr(r=r)), 3)[1]) ≈ 4/3

end

##

@testset "Lensing" begin

    local f,ϕ

    Cℓ = camb().unlensed_total
    nside = 128

    for T in (Float32, Float64)

        @testset "T :: $T" begin

            ϵ = sqrt(eps(T))
            Cϕ = Cℓ_to_Cov(Flat(Nside=nside), T, S0, Cℓ.ϕϕ)
            @test (ϕ = @inferred simulate(Cϕ)) isa FlatS0
            Lϕ = LenseFlow(ϕ)

            ## S0
            Cf = Cℓ_to_Cov(Flat(Nside=nside), T, S0, Cℓ.TT)
            @test (f = @inferred simulate(Cf)) isa FlatS0
            @test (@inferred Lϕ*f) isa FlatS0
            # adjoints
            f,g = simulate(Cf),simulate(Cf)
            @test f' * (Lϕ * g) ≈ (f' * Lϕ) * g
            # gradients
            δf, δϕ = simulate(Cf), simulate(Cϕ)
            @test (FΦTuple(δf,δϕ)'*(δf̃ϕ_δfϕ(Lϕ,Lϕ*f,f)'*FΦTuple(f,ϕ))) ≈ (f'*((LenseFlow(ϕ+ϵ*δϕ)*(f+ϵ*δf))-(LenseFlow(ϕ-ϵ*δϕ)*(f-ϵ*δf)))/(2ϵ)) rtol=1e-2

            # S2 lensing
            Cf = Cℓ_to_Cov(Flat(Nside=nside), T, S2, Cℓ.EE, Cℓ.BB)
            @test (f = @inferred simulate(Cf)) isa FlatS2
            @test (@inferred Lϕ*f) isa FlatS2
            # adjoints
            f,g = simulate(Cf),simulate(Cf)
            @test f' * (Lϕ * g) ≈ (f' * Lϕ) * g
            # gradients
            δf, δϕ = simulate(Cf), simulate(Cϕ)
            @test (FΦTuple(δf,δϕ)'*(δf̃ϕ_δfϕ(Lϕ,Lϕ*f,f)'*FΦTuple(f,ϕ))) ≈ (f'*((LenseFlow(ϕ+ϵ*δϕ)*(f+ϵ*δf))-(LenseFlow(ϕ-ϵ*δϕ)*(f-ϵ*δf)))/(2ϵ)) rtol=1e-2

        end

    end

end

##

@testset "Posterior" begin

    Cℓ = camb()
    L = LenseFlow{RK4Solver{7}}
    T = Float64

    for pol in (:I,:P)

        @testset "pol = $pol" begin

            @unpack f,f̃,ϕ,ds,ds₀ = load_sim_dataset(
                seed  = 0,
                Cℓ    = Cℓ,
                θpix  = 3,
                Nside = 128,
                T     = T,
                beamFWHM = 3,
                pol   = pol,
                L     = L,
                pixel_mask_kwargs = (edge_padding_deg=2,)
                );
            @unpack Cf,Cϕ,D = ds₀
            f° = L(ϕ)*D*f

            @test lnP(0,f,ϕ,ds) ≈ lnP(1,    f̃,  ϕ ,ds) rtol=1e-4
            @test lnP(0,f,ϕ,ds) ≈ lnP(:mix, f°, ϕ, ds) rtol=1e-4

            ε = sqrt(eps(T))
            δf,δϕ = simulate(Cf),simulate(Cϕ)

            @test FΦTuple(gradient((f, ϕ)->lnP(0,   f, ϕ,ds),f, ϕ))'*FΦTuple(δf,δϕ) ≈ (lnP(0,   f +ε*δf,ϕ+ε*δϕ,ds)-lnP(0,   f -ε*δf,ϕ-ε*δϕ,ds))/(2ε)  rtol=1e-2
            @test FΦTuple(gradient((f̃ ,ϕ)->lnP(1,   f̃, ϕ,ds),f̃, ϕ))'*FΦTuple(δf,δϕ) ≈ (lnP(1,   f̃ +ε*δf,ϕ+ε*δϕ,ds)-lnP(1,   f̃ -ε*δf,ϕ-ε*δϕ,ds))/(2ε)  rtol=1e-1
            @test FΦTuple(gradient((f°,ϕ)->lnP(:mix,f°,ϕ,ds),f°,ϕ))'*FΦTuple(δf,δϕ) ≈ (lnP(:mix,f°+ε*δf,ϕ+ε*δϕ,ds)-lnP(:mix,f°-ε*δf,ϕ-ε*δϕ,ds))/(2ε)  rtol=5e-2

        end

    end

end

##

end


text:
using PMDPs
using Test
using POMDPSimulators

using StaticArrays
using POMDPs
using MCTS, DiscreteValueIteration
using StatsBase
# using Plots
using POMDPSimulators
# using D3Trees
using POMDPPolicies
using POMDPLinter
using Random
using DataFrames
using POMDPSimulators

"""
Options
1) Rollout policy
2) Initialize N and Q
"""

mdp_vi = PMDPs.create_PMDP(PMDPe)
mdp_mc = PMDPs.create_PMDP(PMDPg)

policy = PMDPs.get_VI_policy(mdp_vi)
# planner = PMDPs.get_MCTS_planner(mdp_mc)

solver = MCTSSolver(n_iterations=5000,
                    depth=4,
                    exploration_constant=50.0,
                    enable_tree_vis=true)
planner = solve(solver, mdp_mc);

rng = MersenneTwister(1234)

PMDPs.compare_actions(mdp_vi, policy, planner; rng_seed=123, verbose=true)

mdp_mc = PMDPs.create_PMDPe10(PMDPg)
s0 = rand(rng, initialstate(mdp_mc))
solver = MCTSSolver(n_iterations=5000,
                    depth=1,
                    exploration_constant=50.0,
                    enable_tree_vis=true)
planner = solve(solver, mdp_mc);
action(planner, s0)

rng = MersenneTwister(123)
hr = HistoryRecorder(max_steps=1000, capture_exception=true, rng=rng)
h = simulate(hr, mdp_mc, planner)

# PMDPs.run_sim(mdp_mc, planner)
# s0 = rand(rng, initialstate(mdp_mc))

# function run_sim(mdp::PMDP, policy::Policy; rng_seed=1234)
#     rng = MersenneTwister(rng_seed)
#     hr = HistoryRecorder(max_steps=100, capture_exception=true, rng=rng)
#     h = simulate(hr, mdp, policy)
#     collect(eachstep(h, "s, a, r, user_budget"))
#     # sum(h[:r])
# end

# rng = MersenneTwister(12)
# hr = HistoryRecorder(max_steps=100, capture_exception=false, rng=rng)
# h = simulate(hr, mdp_mc, planner)
# collect(eachstep(h, "s, a, r, info"))

# @show run_sim(mdp_mc, policy; rng_seed = 1235)

# ch = run_sim(mdp_mc, planner; rng_seed = 1236)
# @show ch
# any(ch[end][:s].c .< 0)

# for i in 1:10000
#     ch = run_sim(mdp_vi, planner; rng_seed = i)
#     print(i, " ")
#     any(ch[end][:s].c .< 0) ? break : continue
# end
# action(planner, s0)

# params_mdp = Dict(pairs( (n_edges = 3, c_init = 2, demand = Float64[1.5,1.5,1.5], selling_horizon_end = [20,25,30], actions = 15:5:90)))
params_mdp = Dict(pairs( (n_edges = 2, c_init = 1, demand = Float64[1,1], selling_horizon_end = [20,25], actions = 15:5:30)))

mdp_vi = PMDPs.create_PMDP(PMDPe; params_mdp...);
policy = PMDPs.get_VI_policy(mdp_vi);

result, filepath = PMDPs.vi_policy(params_mdp, mdp_vi)


text:
import Markdown: latex, latexinline

function latex(io::IO, tex::Markdown.LaTeX)
    math_envs = ["align", "equation", "eqnarray"]
    use_dollars = !any([occursin("\\begin{$me", tex.formula) for me in math_envs])
    use_dollars && write(io, "\\[")
    write(io, string("\n", tex.formula, "\n"))
    use_dollars && write(io, "\\]\n")
end

#Remove comments that can occur inside a line
function latexinline(io, comment::WeaveMarkdown.Comment)
    write(io, "")
end

function latex(io::IO, comment::WeaveMarkdown.Comment)
    for line in split(comment.text, r"\r\n|\n")
        write(io, "% $line\n")
    end
end


text:
module Deno

using Pkg.Artifacts
using Base: SHA1

export deno_cmd

function load_artifact(artifact_name)
        artifact_info = artifact_meta(artifact_name, joinpath(@__DIR__, "..", "Artifacts.toml"))
        artifact_info === nothing && return nothing

        return artifact_path(SHA1(artifact_info["git-tree-sha1"]))
end

const deno_path = load_artifact("deno_exe")
const deno_exe_name = Sys.iswindows() ? "deno.exe" : "deno"
const deno_exe_path = deno_path === nothing ? deno_exe_name : joinpath(deno_path, deno_exe_name)

"""
Returns the full path to the deno executable
"""
function deno_cmd()
        return `$deno_exe_path`
end

end


text:
using Test

@testset "MathOptInterface" begin
    include("MOI_wrapper.jl")
end


text:
using DashBootstrapComponents, DashHtmlComponents

tabs = html_div([
    dbc_tabs(
        [
            dbc_tab(label = "Tab 1", tab_id = "tab-1"),
            dbc_tab(label = "Tab 2", tab_id = "tab-2"),
        ],
        id = "tabs",
        active_tab = "tab-1",
    ),
    html_div(id = "content"),
]);

callback!(app, Output("content", "children"), Input("tabs", "active_tab")) do at
    if at == "tab-1"
        return tab1_content
    elseif at == "tab-2"
        return tab2_content
    else
        return html_p("This shouldn't ever be displayed...")
    end
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
eval(quote
    using .DynamicPolynomials: @polyvar
end)

eval(load_kron_dynamicpolynomials())


text:
module AnalyzeRegistry

# Standard libraries
using Pkg, TOML
# Third-party packages
using FLoops # for the `@floops` macro
using MicroCollections # for `EmptyVector` and `SingletonVector`
using BangBang # for `append!!`

export general_registry, find_packages, analyze

struct Package
    name::String # name of the package
    repo::String # URL of the repository
    reachable::Bool # can the repository be cloned?
    docs::Bool # does it have documentation?
    runtests::Bool # does it have the test/runtests.jl file?
    github_actions::Bool # does it use GitHub Actions?
    travis::Bool # does it use Travis CI?
    appveyor::Bool # does it use AppVeyor?
    cirrus::Bool # does it use Cirrus CI?
    circle::Bool # does it use Circle CI?
    drone::Bool # does it use Drone CI?
    buildkite::Bool # does it use Buildkite?
    azure_pipelines::Bool # does it use Azure Pipelines?
    gitlab_pipeline::Bool # does it use Gitlab Pipeline?
end
function Package(name, repo;
                 reachable=false,
                 docs=false,
                 runtests=false,
                 github_actions=false,
                 travis=false,
                 appveyor=false,
                 cirrus=false,
                 circle=false,
                 drone=false,
                 buildkite=false,
                 azure_pipelines=false,
                 gitlab_pipeline=false,
                 )
    return Package(name, repo, reachable, docs, runtests, github_actions, travis,
                   appveyor, cirrus, circle, drone, buildkite, azure_pipelines, gitlab_pipeline)
end

function Base.show(io::IO, p::Package)
    body = """
        Package $(p.name):
          * repo: $(p.repo)
          * is reachable: $(p.reachable)
        """
    if p.reachable
        body *= """
              * has documentation: $(p.docs)
              * has tests: $(p.runtests)
            """
        ci_services = (p.github_actions => "GitHub Actions",
                       p.travis => "Travis",
                       p.appveyor => "AppVeyor",
                       p.cirrus => "Cirrus",
                       p.circle => "Circle",
                       p.drone => "Drone CI",
                       p.buildkite => "Buildkite",
                       p.azure_pipelines => "Azure Pipelines",
                       p.gitlab_pipeline => "GitLab Pipeline",
                       )
        if any(first.(ci_services))
            body *= "  * has continuous integration: true\n"
            for (k, v) in ci_services
                if k
                    body *= "    * $(v)\n"
                end
            end
        else
            body *= "  * has continuous integration: false\n"
        end
    end
    print(io, strip(body))
end

"""
    general_registry() -> String

Guess the path of the General registry.
"""
general_registry() =
    first([joinpath(d, "registries", "General") for d in Pkg.depots() if isfile(joinpath(d, "registries", "General", "Registry.toml"))])

"""
    find_packages(dir = general_registry()) -> Vector{String}

Find all packages in the given registry, the General registry by default.
Return a vector with the paths to the directories of each package in the
registry.
"""
function find_packages(dir = general_registry())
    # Get the list of packages in the registry by parsing the `Registry.toml`
    # file in the given directory.
    packages = TOML.parsefile(joinpath(general_registry(), "Registry.toml"))["packages"]
    # Get the directories of all packages.  Filter out JLL packages: they are
    # automatically generated and we know that they don't have testing nor
    # documentation.
    packages_dirs = [joinpath(dir, p["path"]) for (_, p) in packages if !endswith(p["name"], "_jll")]
end

"""
    analyze(dir::AbstractString) -> Package

Analyze the package whose entry in the registry is in the `dir` directory.
NOTE: the git repository of the package will be cloned, so you need Internet
access to use this package.

## Example
```julia
julia> analyze(joinpath(general_registry(), "B", "BinaryBuilder"))
Package BinaryBuilder:
  * repo: https://github.com/JuliaPackaging/BinaryBuilder.jl.git
  * is reachable: true
  * has documentation: true
  * has tests: true
  * has continuous integration: true
    * GitHub Actions
    * Azure Pipelines
```
"""
function analyze(dir::AbstractString)
    # Parse the `Package.toml` file in the given directory.
    toml = TOML.parsefile(joinpath(dir, "Package.toml"))
    name = toml["name"]::String
    repo = toml["repo"]::String
    # Clone the repository into a temporary directory
    package = mktempdir() do dest
        reachable = try
            # Clone only latest commit on the default branch.  Note: some
            # repositories aren't reachable because the author made them private
            # or deleted them.  In these cases git would ask for username and
            # password, provide it with fake values just to move on:
            # https://stackoverflow.com/a/65705346/2442087
            run(pipeline(`git clone -q --depth 1 --config credential.helper='!f() { echo -e "username=git\npassword="; }; f' $(repo) $(dest)`; stderr=devnull))
            true
        catch
            # The repository may be unreachable
            false
        end
        if reachable
            docs = isfile(joinpath(dest, "docs", "make.jl")) || isfile(joinpath(dest, "doc", "make.jl"))
            testpath = joinpath(dest, "test", "runtests.jl")
            runtests = isfile(testpath) && (length(readlines(testpath))>6 || length(readdir(joinpath(dest, "test")))>1)
            travis = isfile(joinpath(dest, ".travis.yml"))
            appveyor = isfile(joinpath(dest, "appveyor.yml"))
            cirrus = isfile(joinpath(dest, ".cirrus.yml"))
            circle = isfile(joinpath(dest, ".circleci", "config.yml"))
            drone = isfile(joinpath(dest, ".drone.yml"))
            azure_pipelines = isfile(joinpath(dest, "azure-pipelines.yml"))
            buildkite = isfile(joinpath(dest, ".buildkite", "pipeline.yml"))
            gitlab_pipeline = isfile(joinpath(dest, ".gitlab-ci.yml"))
            github_workflows = joinpath(dest, ".github", "workflows")
            if isdir(github_workflows)
                # Find all workflows
                files = readdir(github_workflows)
                # Exclude TagBot and CompatHelper
                filter(f -> lowercase(f) ∉ ("compathelper.yml", "tagbot.yml"), files)
                # Assume all other files are GitHub Actions for CI.  May not
                # _always_ be the case, but it's a good first-order approximation.
                github_actions = length(files) > 0
            else
                github_actions = false
            end
            Package(name, repo; reachable, docs, runtests, travis, appveyor, cirrus,
                    circle, drone, buildkite, azure_pipelines, gitlab_pipeline, github_actions)
        else
            Package(name, repo)
        end
    end
    return package
end

"""
    analyze(packages::AbstractVector{<:AbstractString}) -> Vector{Package}

Analyze all packages in the iterable `packages`, using threads.  Return a
`Vector{Package}`.
"""
function analyze(packages::AbstractVector{<:AbstractString})
    @floop for p in packages
        ps = SingletonVector((analyze(p),))
        @reduce(result = append!!(EmptyVector(), ps))
    end
    result
end

end # module


text:
include("generator_utils.jl")

using REPL
write_dict_to_file("repl_latex_symbols.jl",
                   "repl_latex_symbols",
                   REPL.REPLCompletions.latex_symbols)


text:
# Macros and functions to check if the BigNumber RNG is defined.
# Based on _os_test from osutils.jl
function bigrng_test(qm, ex, test)
    @assert qm == :?
    @assert isa(ex, Expr)
    @assert ex.head == :(:)
    @assert length(ex.args) == 2
    if test
        return esc(ex.args[1])
    else
        return esc(ex.args[2])
    end
end

macro bigrng(qm, ex)
    bigrng_test(qm, ex, isdefined(:BigRNG))
end

macro have_bigrng(ex)
    @bigrng? esc(ex) : nothing
end

macro no_bigrng(ex)
    @bigrng? nothing : esc(ex)
end


text:
function Base.convert(::Type{cusolverEigType_t}, typ::Int)
    if typ == 1
        CUSOLVER_EIG_TYPE_1
    elseif typ == 2
        CUSOLVER_EIG_TYPE_2
    elseif typ == 3
        CUSOLVER_EIG_TYPE_3
    else
        throw(ArgumentError("Unknown eigenvalue solver type $typ."))
    end
end

function Base.convert(::Type{cusolverEigMode_t}, jobz::Char)
    if jobz == 'N'
        CUSOLVER_EIG_MODE_NOVECTOR
    elseif jobz == 'V'
        CUSOLVER_EIG_MODE_VECTOR
    else
        throw(ArgumentError("Unknown eigenvalue solver mode $jobz."))
    end
end


text:
function nonsym_qstateclock_mpo(beta::Float64,q::Int)
    dat = zeros(ComplexF64,q,q,q,q);
    comega(d) = cos(2*pi*d/q);

    for i = 1:q, j = 1:q, k = 1:q, l = 1:q
        dat[i,j,k,l] = exp(beta*(comega(i-j) + comega(j-k) + comega(k-l) + comega(l-i)));
    end

    return PeriodicMPO(TensorMap(dat,ℂ^q*ℂ^q,ℂ^q*ℂ^q))
end


text:
using Test
using DFTK

include("testcases.jl")

# Silicon redHF (without xc) is a metal, so we add a bit of temperature to it

# TODO There is a lot of code duplication in this file ... once we have the ABINIT reference
#      stuff in place, this should be refactored.

function run_silicon_redHF(T; Ecut=5, test_tol=1e-6, n_ignored=0, grid_size=15, scf_tol=1e-6)
    # T + Vloc + Vnloc + Vhartree
    # These values were computed using ABINIT with the same kpoints as testcases.jl
    # and Ecut = 25
    ref_redHF = [
        [0.17899118507651615, 0.6327279881297371, 0.6327279881326648, 0.6327279881356039,
         0.706557757783828, 0.7065577577877139, 0.7065577577915956, 0.7397951816714727,
         0.8532089291297222, 0.8978914445971602],
        [0.23220003663858457, 0.42189409862896016, 0.5921574659414509, 0.5921574659446628,
         0.672858189872362, 0.7372271903827399, 0.7372271903861028, 0.8643640848936627,
         0.9011792204214196, 0.9011792204356576],
        [0.2517502116803524, 0.445206025448218, 0.5328870916963034, 0.532887091701182,
         0.6211365856991057, 0.661989858948651, 0.8863951918546257, 0.8863951918584175,
         0.973261179805555, 0.9771287508158364],
        [0.30685586314464863, 0.376375429632464, 0.4438764716222098, 0.5459065154292047,
         0.651122698647485, 0.8164293660861612, 0.8515978828421051, 0.8735213568005982,
         0.8807275612483988, 0.8886454931307763]
    ]
    ref_etot = -5.440593269861395

    n_bands = length(ref_redHF[1])
    fft_size = grid_size * ones(3)
    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.atnum, functional="lda", family="hgh"))
    model = model_DFT(Array{T}(silicon.lattice), [Si => silicon.positions], [];
                      temperature=0.05)
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    scfres = self_consistent_field(basis, tol=scf_tol, n_bands=n_bands)

    for ik in 1:length(silicon.kcoords)
        @test eltype(scfres.eigenvalues[ik]) == T
        @test eltype(scfres.ψ[ik]) == Complex{T}
        println(ik, "  ", abs.(ref_redHF[ik] - scfres.eigenvalues[ik][1:n_bands]))
    end
    for ik in 1:length(silicon.kcoords)
        # Ignore last few bands, because these eigenvalues are hardest to converge
        # and typically a bit random and unstable in the LOBPCG
        diff = abs.(ref_redHF[ik] - scfres.eigenvalues[ik][1:n_bands])
        @test maximum(diff[1:n_bands - n_ignored]) < test_tol
    end
end

function run_silicon_lda(T ;Ecut=5, test_tol=1e-6, n_ignored=0, grid_size=15, scf_tol=1e-6,
                         n_noconv_check=0)
    # These values were computed using ABINIT with the same kpoints as testcases.jl
    # and Ecut = 25
    ref_lda = [
        [-0.178566465714968, 0.261882541175914, 0.261882541178847, 0.261882541181782,
          0.354070367072414, 0.354070367076363, 0.354070367080310, 0.376871160884678],
        [-0.127794342370963, 0.064395861472044, 0.224958824747686, 0.224958824750934,
          0.321313617512188, 0.388442495007398, 0.388442495010722, 0.542078732298094],
        [-0.108449612789883, 0.077125812982728, 0.172380374761464, 0.172380374766260,
          0.283802499666810, 0.329872296009131, 0.525606867582028, 0.525606867585921],
        [-0.058089253154566, 0.012364292440522, 0.097350168867990, 0.183765652148129,
          0.314593174568090, 0.470869435132365, 0.496966579772700, 0.517009645871194],
    ]
    ref_etot = -7.911817522631488
    n_bands = length(ref_lda[1])

    fft_size = grid_size * ones(3)
    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.atnum, functional="lda", family="hgh"))
    model = model_DFT(Array{T}(silicon.lattice), [Si => silicon.positions], [:lda_x, :lda_c_vwn])
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    scfres = self_consistent_field(basis; n_bands=n_bands, tol=scf_tol,
                                   eigensolver=lobpcg_hyper, n_ep_extra=n_noconv_check)

    for ik in 1:length(silicon.kcoords)
        @test eltype(scfres.eigenvalues[ik]) == T
        @test eltype(scfres.ψ[ik]) == Complex{T}
    end
    for ik in 1:length(silicon.kcoords)
        # Ignore last few bands, because these eigenvalues are hardest to converge
        # and typically a bit random and unstable in the LOBPCG
        @test maximum(ref_lda[ik][1:n_bands - n_ignored] - scfres.eigenvalues[ik][1:n_bands - n_ignored]) < test_tol
    end

    energies = scfres.energies
    @test sum(values(energies)) ≈ ref_etot atol=test_tol
    @test eltype(sum(values(scfres.energies))) == T
end

function run_silicon_pbe(T ;Ecut=5, test_tol=1e-6, n_ignored=0, grid_size=15, scf_tol=1e-6)
    # These values were computed using ABINIT with the same kpoints as testcases.jl
    # and Ecut = 25
    ref_pbe = [
        [-0.181210259413818, 0.258840553222639, 0.258840553225549, 0.258840553228459,
          0.351692348652324, 0.351692348656259, 0.351692348660193, 0.380606400669216,
          0.540705881744348, 0.540705883460555],
        [-0.130553299114991, 0.062256443775155, 0.221871391287580, 0.221871391290802,
          0.322398722411882, 0.386194327436667, 0.386194327439986, 0.546859898649217,
          0.550571701390781, 0.550571701394327],
        [-0.111170738096744, 0.074494899973125, 0.169461730083372, 0.169461730088140,
          0.284305392082236, 0.330468937070505, 0.524509288492752, 0.524509288496625,
          0.616964090764029, 0.619623658242765],
        [-0.061054203629684, 0.009700769243041, 0.095769985640881, 0.180784778430457,
          0.315000287382235, 0.471042322838057, 0.495281775946584, 0.517469860611792,
          0.530124341745161, 0.539044739392045],
    ]
    ref_etot = -7.854477356672080
    n_bands = length(ref_pbe[1])

    fft_size = grid_size * ones(3)
    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.atnum, functional="pbe", family="hgh"))
    model = model_DFT(Array{T}(silicon.lattice), [Si => silicon.positions], [:gga_x_pbe, :gga_c_pbe])
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    scfres = self_consistent_field(basis; tol=scf_tol, n_bands=n_bands)

    for ik in 1:length(silicon.kcoords)
        @test eltype(scfres.eigenvalues[ik]) == T
        @test eltype(scfres.ψ[ik]) == Complex{T}
        println(ik, "  ", abs.(ref_pbe[ik] - scfres.eigenvalues[ik][1:n_bands]))
    end
    for ik in 1:length(silicon.kcoords)
        # Ignore last few bands, because these eigenvalues are hardest to converge
        # and typically a bit random and unstable in the LOBPCG
        diff = abs.(ref_pbe[ik] - scfres.eigenvalues[ik][1:n_bands])
        @test maximum(diff[1:n_bands - n_ignored]) < test_tol
    end

    energies = scfres.energies
    @test sum(values(energies)) ≈ ref_etot atol=test_tol
    @test eltype(sum(values(scfres.energies))) == T
end


text:
input = "input/shell_sort.nl"
manual = false #If it's true, all the steps on the analysis will require the user to press enter. Otherwise, it will be executed automatically

input_name = split(split(input,"/")[2],".")[1] #The name without the '.nl' and without the directory
tree_file = open("../outputs/tree/$input_name","w+")

if !isdefined(:already_executed_flag) #Only calls the includes if it's the first time you call include("main.jl")
	try
		info("Loading the needed packages. Since it's the first execution of the code for this workspace, it may take some seconds.")
		info("On the next time you call call include(\"main.jl\"), it will be faster since the packages will not be loaded again.")

		using CSV
		#using DataStructures
		using DataFrames
		using JSON
		using Match
	catch
		info("Installing missing packages. It may take some minutes depending on your internet connection.")
		Pkg.add("CSV")
		#Pkg.add("DataStructures")
		Pkg.add("DataFrames")
		Pkg.add("JSON")
		Pkg.add("Match")

		using CSV
		#using DataStructures
		using JSON
		using Match
	end
end
already_executed_flag =  true

include("tokens.jl")
include("data_structures.jl")
include("auxiliar_funcs.jl")
include("grammar.jl")

include("lexical.jl")
table = CSV.read("grammar_table.csv")

include("sinthatic.jl")

close(f)
close(tree_file)
close(g)

#Cria uma cópia da tabela LL1 na pasta especificações, onde o professor vai ver o que é pedido no arquivo
#Ocasionalmente esta tabela pode ser atualizada/melhorada. Por isso esse comando para sempre enviar ela para a pasta da avaliação
cp("grammar_table.csv","../especificações/tabela_LL1.csv",remove_destination=true)

#Converte a gramática da sintaxe usada na linguagem para a sintaxe vista em sala de aula
#convertGrammar()

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler