text:
using PyPlot
using Triangulation_diagrams

N = Float64[ 1  0  1  0  -1   0  -1
             1  1  0  0   0  -1  -1 ]

T = [ 1 4 4 3 4 5
      2 3 2 4 5 7
      3 2 5 6 6 6 ]

E = [ 6 3 1 2 5 7
      3 1 2 5 7 6 ]

tri = Triangulation(N, T, E, 4)

figure(1)
draw_triangles(tri)
enumerate_triangles(tri)
enumerate_vertices(tri)
mark_first_local_vertex(tri, 0.15)
enumerate_bdry_edges(tri, 0.1)
axis("off")
axis("equal")
savefig("ex1_triangulation.pdf")

N = Float64[ 4  2  0  1  3  5  2  4  6  6
             2  2  2  1  1  1  0  0  0  2 ]

T = [ 3  5  4  7  6  1  5   6  10  8
      4  2  7  8  5  2  6   9   1  9
      2  4  5  5  8  5  1  10   6  6 ]

E = [ 10  1  2  3  4  7  8   9
       1  2  3  4  7  8  9  10 ]

tri = Triangulation(N, T, E, 5)

figure(2)
draw_triangles(tri)
enumerate_triangles(tri)
enumerate_vertices(tri)
mark_first_local_vertex(tri, 0.3)
enumerate_bdry_edges(tri, 0.2)
axis("off")
axis("equal")
savefig("ex2_triangulation.pdf")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
using Test, QSimulator
using QSimulator: index
using LsqFit: curve_fit

############################ Ramsey ##########################

# ramsey evolution of a single transmon should produce sinusoidal oscillations at the qubit frequency
@testset "Ramsey" begin
    dims = 3
    qubit_freq = 5.0
    anharm = -0.2
    q0 = DuffingTransmon("q0", dims, DuffingSpec(qubit_freq, anharm))
    cqs = CompositeQSystem([q0])
    add_hamiltonian!(cqs, q0)
    times = collect(range(0,stop=1,length=201))
    ψ_init = (1/sqrt(2)) * ComplexF64[1; 1; 0]
    ψs = unitary_state(cqs, times, ψ_init)
    signal = Float64[real(ψ_init'*ψ) for ψ in ψs]
    expected = 0.5 .+ 0.5*cos.(2π*qubit_freq * times)
    @test isapprox(signal, expected, rtol=1e-3)

    # check that the propagator gives the same result
    us = unitary_propagator(cqs, times)
    ψs = [u*ψ_init for u in us]
    signal = Float64[real(ψ_init'*ψ) for ψ in ψs]
    @test isapprox(signal, expected, rtol=1e-4)
end

########################### Rabi ##########################

# rabi flops should produce sinusoidal oscillations between ground and excited state
@testset "Rabi" begin
    dims = 3
    qubit_freq = 5.0
    anharm = -0.2
    nutation_freq = 0.02
    q0 = DuffingTransmon("q0", dims, DuffingSpec(qubit_freq, anharm))
    cqs = CompositeQSystem([q0])
    add_hamiltonian!(cqs, q0)
    add_hamiltonian!(cqs, dipole_drive(q0, t -> nutation_freq*cos(2π*qubit_freq * t)), q0)
    ψ_init = ComplexF64[1; 0; 0]
    times = collect(range(0,stop=100,length=101))
    ψs = unitary_state(cqs, times, ψ_init)

    g_sim = [abs2(s[1]) for s in ψs]
    e_sim = [abs2(s[2]) for s in ψs]

    # we would normally expect a factor of 0.5 from lab -> rotating frame but by defining the qubit
    # drive in terms of X rather than 0.5 X undoes that
    g_expected = 0.5 .+ 0.5*cos.(2π*nutation_freq * times)
    @test isapprox(g_sim, g_expected, rtol=1e-2)
    @test isapprox(e_sim, 1 .- g_expected, rtol=1e-2)

    #rwa version
    cqs_rwa = CompositeQSystem([q0])
    add_hamiltonian!(cqs_rwa, q0)
    add_hamiltonian!(cqs_rwa, rwa_dipole(q0, t -> nutation_freq*cos(2π*qubit_freq * t)), q0)
    ψs_rwa = unitary_state(cqs_rwa, times, ψ_init)

    g_sim_rwa = [abs2(s[1]) for s in ψs_rwa]
    e_sim_rwa = [abs2(s[2]) for s in ψs_rwa]

    @test isapprox(g_sim_rwa, g_expected, rtol=1e-2)
    @test isapprox(e_sim_rwa, 1 .- g_expected, rtol=1e-2)

    # check that the propagator gives the same result
    # use floquet_propagator to speed it up since the Hamiltonian is periodic
    floquet_prop = floquet_propagator(unitary_propagator, 1/qubit_freq)
    us = floquet_prop(cqs, times)
    ψs = [u*ψ_init for u in us]
    g_sim = [abs2(s[1]) for s in ψs]
    e_sim = [abs2(s[2]) for s in ψs]
    @test isapprox(g_sim, g_expected; rtol=1e-2)
    @test isapprox(e_sim, 1 .- g_expected; rtol=1e-2)
end

########################### Parmetric Flux Drive ##########################

# parametric flux drive should produce flops between resonant states parameters from S. A.
# Caldwell, et. al (2018). Parametrically Activated Entangling Gates Using Transmon Qubits. Physical
# Review Applied, 10(3). https://doi.org/10.1103/PhysRevApplied.10.034050

@testset "parametric flux drive" begin
    # create Hamiltonian
    dims = [3,3]
    duffing = DuffingSpec(3.94015, -0.1807)
    transmon = TransmonSpec(0.172, 12.71, 3.69)
    q0 = DuffingTransmon("q0", dims[1], duffing)
    q1 = PerturbativeTransmon("q1", dims[2], transmon)
    # go to rotating frame at duffing.frequency for speed
    # a doubly rotating frame is also possible, but then you can't use floquet
    q0 = RotatingFrameSystem(q0, duffing.frequency)
    q1 = RotatingFrameSystem(q1, duffing.frequency)
    cqs = CompositeQSystem([q0, q1])
    add_hamiltonian!(cqs, q0)
    g = 0.006
    add_hamiltonian!(cqs, t -> g/2 * X_Y([q0, q1], [duffing.frequency, duffing.frequency] * t), [q0,q1])

    # compute iSWAP parameters
    mod_park = 0.0
    mod_amp = 0.323
    ϕ(t, freq) = mod_park + mod_amp * sin(2π * freq * t)
    f(t, freq) = perturbative_transmon_freq(transmon.EC, transmon.EJ1, transmon.EJ2, ϕ(t, freq))
    dummy_freq = 1.0 # use this for computing average frequency
    avg_freq = real(fourier_coefficient(t -> f(t, dummy_freq), dummy_freq, 0))
    detuning = avg_freq - duffing.frequency
    harmonic = 2
    freq = -detuning/harmonic
    fs = FourierSeries(t -> 2π * f(t, freq), freq, collect(-50:50))
    g_eff = abs(g * rotating_frame_series(fs, [harmonic]).terms[harmonic])
    t_gate = 1/(4 * g_eff)

    # create basis states
    basis = TensorProductBasis(dims)
    ψ₀₁ = TensorProductBasisState(basis, (0,1))
    ψ₁₀ = TensorProductBasisState(basis, (1,0))

    # perform time evolution
    add_hamiltonian!(cqs, parametric_drive(q1, t -> ϕ(t, freq)), q1)
    times = collect(range(0, stop=2 * t_gate, length=1000))
    prop = floquet_propagator(unitary_propagator, 1/abs(freq))
    ψs = [u * vec(ψ₁₀) for u in prop(cqs, times)]
    pop_10 = [abs2(ψ[index(ψ₁₀)]) for ψ in ψs]
    pop_01 = [abs2(ψ[index(ψ₀₁)]) for ψ in ψs]
    # population should oscillate between 01 and 10
    expected_10 = cos.(2π * g_eff * times).^ 2
    # there is additional lab frame jaggedness so relax tolerance
    @test isapprox(pop_10, expected_10, rtol=2e-2)
    @test isapprox(pop_01, 1 .- expected_10, rtol=2e-2)
end

########################## ME Solver #################################
# This tests the decay of a driven two level system. Added to the hamiltonian is
# is a drive term of the form (J σ exp(iωt) + h.c.). When ω is chosen to
# coincide with the splitting of the two level system, J expressed in Hz is the
# rate of oscillation. For a system that only has T1 decay, the exponential rate
# of decay is Td = ⁴/₃ T1.

@testset "decaying driven two levels" begin
    dims = 3
    qubit_freq = 5.0
    anharm = -0.2
    q0 = DuffingTransmon("q0", dims, DuffingSpec(qubit_freq, anharm))
    cqs = CompositeQSystem([q0])
    add_hamiltonian!(cqs, q0)
    add_hamiltonian!(cqs, dipole_drive(q0, t -> 0.02*cos(2π*qubit_freq * t)), q0)

    T1 = 50. # in ns
    γ = 1/(2π * T1) # in GHz
    add_lindblad!(cqs, decay(q0, γ), [q0])

    ψ0 = ComplexF64[1; 0; 0]
    ρ0 = ψ0 * ψ0'
    times = collect(range(0, stop=100, length=101))
    # use floquet and me_propagator for speed
    prop = floquet_propagator(me_propagator, 1/qubit_freq)
    ρs = [reshape(u * vec(ρ0), dims, dims) for u in prop(cqs, times)]

    ρ00 = [real(ρ[1, 1]) for ρ in ρs]

    model_t2(x, p) = (exp.(-x ./ p[1]) .* cos.(2π .* p[2] .* x .- p[3]) .+ 1.0) / 2.0

    fit = curve_fit(model_t2, times, ρ00, [T1, .02, 0])
    @test abs(fit.param[1] - 4. * T1 / 3.) < 1.
end

ext: jl


text:
module LibRealSense

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibRealSense was not build properly. Please run Pkg.build(\"LibRealSense\").")
end
include(depsjl_path)

# Module initialization function
function __init__()
    check_deps()
end

include("CEnum.jl")
using .CEnum

include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "rs2_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "rs2_api.jl"))

foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "rs2_") || startswith(string(s), "RS2_")
        @eval export $s
    end
end

const RS2_API_VERSION = RS2_API_MAJOR_VERSION * 10000 + RS2_API_MINOR_VERSION * 100 + RS2_API_PATCH_VERSION
const RS2_API_VERSION_STR = "$(RS2_API_MAJOR_VERSION).$(RS2_API_MINOR_VERSION).$(RS2_API_PATCH_VERSION)"

export RS2_API_VERSION, RS2_API_VERSION_STR

end # module


text:
using TimerOutputs, BenchmarkTools

"""
    assess_source_reg_quality(target, source_reg; title="", display::Bool=true)

Run a few tests on comparing `target` and `source_reg`.

Tests: [`ncc`, `mae`, `rmse`, `mse`]
"""
function assess_source_reg_quality(target, source_reg; title="", display::Bool=true)
    functions = [ncc, mae, rmse, mse]
    names_short = ["ncc", "mae", "rmse", "mse"]
    vals = map(x -> x(target, source_reg), functions)
    names_vals_dict = Dict(names_short[i] => vals[i] for i in 1:4)
    return names_vals_dict
end

"""
    assess_flow_quality(flow, flow_est; title="", display::Bool=true)

Run a few tests on comparing `flow_est` and `flow`.

Tests: [`angle_mae`, `angle_rmse`, `mae`, `mse`]
"""
function assess_flow_quality(flow, flow_est; title="", display::Bool=true)
    functions = [angle_mae, angle_rmse, mae, mse]
    names_short = ["angle-mae", "angle-rmse", "flow_mae", "flow_rmse"]
    # names = ["angle mean absolute error", "angle root mean squared error", "mean absolute error", "root mean squared error"]
    vals = map(x -> x(flow, flow_est), functions)
    names_vals_dict = Dict(names_short[i] => vals[i] for i in 1:4)
    return names_vals_dict
end

# function average_dicts(dicts)
#     avg_dict = Dict()
#     bad_count = 0
#     for dict in dicts
#         for key in keys(dict)
#             if dict[key] == NaN
#                 bad_count += 1
#                 break;
#             end
#             if !haskey(avg_dict, key)
#                 avg_dict[key] = dict[key]
#             else
#                 avg_dict[key] += dict[key]
#             end
#         end
#     end
#     for key in keys(avg_dict)
#         avg_dict[key] = avg_dict[key]/(length(dicts)-bad_count)
#     end
#     return avg_dict
# end

"""
    fun_on_dict_values(dicts, fun)

Run function `fun` on all same keys of all dictionaries `dicts`.

```jldoctest
using BenchmarkTools
fun_on_dict_values((Dict(:a => 1), Dict(:a => 12)), mean)

# output

Dict{Any,Any} with 1 entry:
  :a => 6.5

```
"""
function fun_on_dict_values(dicts, fun)
    out_dict = Dict()

    for key in keys(dicts[1])
        out_dict[key] = fun([dict[key] for dict in dicts])
    end
    return out_dict
end

fun_on_dict_values((Dict(:a => 1), Dict(:a => 12)), mean)

"""
    compare_dicts(old_dict, new_dict)

Compare two dicts by dividing their corresponding values, return the comparation dict.

Note: (old / new)
"""
function compare_dicts(old_dict, new_dict)
    comparation_dict = Dict()
    for key in keys(old_dict)
        if haskey(new_dict, key)
            comparation_dict[key] = old_dict[key]/new_dict[key]
        end
    end
    return comparation_dict
end


text:
struct DashConfig
    external_stylesheets ::Vector{ExternalSrcType}
    external_scripts ::Vector{ExternalSrcType}
    url_base_pathname ::Union{String, Nothing} #TODO This looks unused
    requests_pathname_prefix ::String
    routes_pathname_prefix ::String
    assets_folder ::String
    assets_url_path ::String
    assets_ignore ::String
    serve_locally ::Bool
    suppress_callback_exceptions ::Bool
    prevent_initial_callbacks ::Bool
    eager_loading ::Bool
    meta_tags ::Vector{Dict{String, String}}
    assets_external_path ::Union{String, Nothing}
    include_assets_files ::Bool
    show_undo_redo ::Bool
    compress ::Bool
end


text:
# # LaserTypes

module LaserTypes

export E, B,
    GaussLaser, LaguerreGaussLaser,
    ConstantProfile, GaussProfile, Cos²Profile, QuasiRectangularProfile,
    setup_laser, Fμν, S

using Unitful
using UnitfulAtomic
using UnPack
using LinearAlgebra
using HypergeometricFunctions
using StaticArrays
using CoordinateTransformations
using AutoHashEquals
using ParallelProcessingTools
import PhysicalConstants.CODATA2018: c_0, e, m_e, ε_0, μ_0

const _₁F₁ = HypergeometricFunctions.drummond1F1
const pochhammer = HypergeometricFunctions.pochhammer

abstract type AbstractLaser end
abstract type AbstractTemporalProfile end

include("envelopes.jl")
include("constants.jl")
include("coords.jl")
include("polarization.jl")
include("electricfield.jl")
include("magneticfield.jl")
include("faraday.jl")
include("gauss.jl")
include("laguerre-gauss.jl")
include("setup.jl")
include("derived.jl")
include("utils.jl")
include("show.jl")

end # module


text:
# functions to download the data

const FILESIZES = CSV.read(joinpath(DIR,"metadata","filesizes.csv"), DataFrame)

# _download is a modified MIT expat licensed code from HTTP.jl
# https://github.com/JuliaWeb/HTTP.jl/blob/668e7e68747bb333ebde13af8d16add5b82b3b8a/src/download.jl#L92

function _download(url::AbstractString,
                   file_path::AbstractString;
                   headers = [],
                   progress = nothing, # nothing or Progress or ParallelProgress
                   filesize = typemax(Int),
                   update_period = 1,
                   kw...)

    if isnothing(progress)
        update_period = Inf
    end

    HTTP.open("GET", url, headers; kw...) do stream
        eof(stream) && return  # don't do anything for streams we can't read (yet)

        downloaded_bytes = 0
        prev_time = Dates.now()

        function report_callback()
            prev_time = Dates.now()
            completion_progress = downloaded_bytes
            # we stay stuck at 99% if file is bigger than expected
            update!(progress, min(completion_progress, filesize-1))
        end

        Base.open(file_path, "w") do fh
            while !eof(stream)
                downloaded_bytes += write(fh, readavailable(stream))
                if !isinf(update_period)
                    if Dates.now() - prev_time > Dates.Millisecond(round(1000update_period))
                        report_callback()
                    end
                end
            end
        end
        if !isinf(update_period)
            finish!(progress)
        end
    end
    return
end

function get_file_size(siteid::Integer, compressed::Bool)
    col = compressed ? ".csv.gz" : ".csv"
    return ceil(Int, FILESIZES[siteid, col])
end

function download_site_csv(source, siteid, path_to_data_folder; kw...)
    if source === :zenodo
        return download_site_csv_from_zenodo(siteid, path_to_data_folder; kw...)
    elseif source === :schneider
        return download_site_csv_from_schneider(siteid, path_to_data_folder; kw...)
    else
        error("Source $(repr(source)) not found")
    end
end

function download_sites_data(path_to_data_folder,
                             sitesid = 1:70;
                             source = :zenodo,
                             kw...)
    for siteid in sitesid
        download_site_csv(source, siteid, path_to_data_folder; kw...)
    end
    return
end

function download_pv(path_to_data_folder; kw...)
    download_pv_csv_from_zenodo(path_to_data_folder; kw...)
    return
end

function download_sites_data_parallel(path_to_data_folder,
                                      sitesid = 1:70;
                                      progress = true,
                                      max_threads = 4,
                                      compressed = true,
                                      source = :zenodo,
                                      kw...)

    file_sizes = get_file_size.(sitesid, compressed)

    if progress
        ext = compressed ? ".csv.gz" : ".csv"
        mprog = MultipleProgress(length(sitesid), file_sizes;
                                 kws = [(:desc => "Downloading $j$ext ",) for j in sitesid],
                                 desc = "Downloading files... ",
                                 dt = 0.1)
    end

    i = firstindex(sitesid) - 1
    @sync for p in 1:max_threads
        @async while true
            idx = (i += 1)
            idx > lastindex(sitesid) && break
            download_site_csv(source, sitesid[idx], path_to_data_folder;
                              progress = progress ? mprog[idx] : false,
                              compressed = compressed,
                              file_size = file_sizes[idx],
                              kw...)
        end
    end
    progress && finish!(mprog)
    return
end


text:
@testset "step!" begin
    # get environment and simulate
    env = Dojo.get_environment("pendulum",
        timestep=0.1,
        gravity=0.0);
    reset(env);

    # step (no control)
    z1 = Dojo.get_maximal_state(env.mechanism)
    u1 = zeros(Dojo.input_dimension(env.mechanism))
    z2 = Dojo.step!(env.mechanism, z1, u1)
    @test norm(z2 - z1) < 1.0e-6

    # step (control)
    z1 = Dojo.get_maximal_state(env.mechanism)
    u1 = rand(Dojo.input_dimension(env.mechanism))
    z2 = Dojo.step!(env.mechanism, z1, u1)
    @test norm(z2 - z1) > 1.0e-6
end

@testset "Storage" begin
    # get environment and simulate
    env = Dojo.get_environment("pendulum",
        timestep=0.1);
    reset(env);
    Dojo.initialize_pendulum!(env.mechanism,
        angle=0.25 * π)
    storage = Dojo.simulate!(env.mechanism, 1.0,
        record=true,
        verbose=false)
    # check storage length
    @test Dojo.length(storage) == 10
    # convert to vector of vectors
    z = Dojo.get_maximal_state(storage)
    @test length(z) == 10
    # convert back to storage
    s = Dojo.generate_storage(env.mechanism, z)
    @test typeof(s) <: Dojo.Storage
    @test Dojo.length(s) == 10
end


text:
all_stimuli = [
    (word=:bank, number=0), (word=:card, number=1), (word=:dart, number=2), (word=:face, number=3),
    (word=:game, number=4), (word=:hand, number=5), (word=:jack, number=6), (word=:king, number=7),
    (word=:lamb, number=8), (word=:mask,number= 9), (word=:neck, number=0), (word=:pipe, number=1),
    (word=:quip, number=2), (word=:rope, number=3), (word=:sock, number=4), (word=:tent,number=5),
    (word=:vent, number=6), (word=:wall, number=7), (word=:xray,number= 8), (word=:zinc, number=9)
]


text:
using NextGenSeqUtils
using Distances

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "nextgensequtils" begin
    @testset "align" begin
        include("test_align.jl")
    end

    @testset "demux" begin
        include("test_demux.jl")
    end

    @testset "evo_dist" begin
        include("test_evodist.jl")
    end

    @testset "FAD" begin
        include("test_FAD.jl")
    end

    @testset "hmm" begin
        include("test_hmm.jl")
    end

    @testset "kmers" begin
        include("test_kmers.jl")
    end

    @testset "orient" begin
        include("test_orient.jl")
    end

    @testset "phreds" begin
        include("test_phreds.jl")
    end

    @testset "utils" begin
        include("test_utils.jl")
    end

end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_energy_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfPower::DataFrame, dfPrice::DataFrame, dfCharge::DataFrame)

Function for writing energy revenue from the different generation technologies.
"""
function write_energy_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfPower::DataFrame, dfPrice::DataFrame, dfCharge::DataFrame)
	dfGen = inputs["dfGen"]
	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	# the price is already US$/MWh, and dfPower and dfCharge is already in MW, so no scaling is needed
	dfEnergyRevenue = DataFrame(Resource = inputs["RESOURCES"], AnnualSum = Array{Union{Missing,Float32}}(undef, G))
	# initiation
	i = 1
	dfEnergyRevenue_ = (
		DataFrame(
			[[names(dfPower)]; collect.(eachrow(dfPower))],
			[:column; Symbol.(axes(dfPower, 1))])[4:T+3,2] .*
	DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[1,:][:Zone]+1].*
	inputs["omega"])
	if i in inputs["FLEX"]
		dfEnergyRevenue_ = (DataFrame([[names(dfCharge)]; collect.(eachrow(dfCharge))], [:column; Symbol.(axes(dfCharge, 1))])[4:T+3,2] .*
		DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[1,:][:Zone]+1].*
		inputs["omega"])
	end
	for i in 2:G
		if i in inputs["FLEX"]
			dfEnergyRevenue_1 = (DataFrame([[names(dfCharge)]; collect.(eachrow(dfCharge))], [:column; Symbol.(axes(dfCharge, 1))])[4:T+3,i+1] .*
			DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1].*
			inputs["omega"])
		else
			dfEnergyRevenue_1 = (DataFrame([[names(dfPower)]; collect.(eachrow(dfPower))], [:column; Symbol.(axes(dfPower, 1))])[4:T+3,i+1] .*
			DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1].*
			inputs["omega"])
		end
		dfEnergyRevenue_ = hcat(dfEnergyRevenue_, dfEnergyRevenue_1)
	end
	dfEnergyRevenue = hcat(dfEnergyRevenue, DataFrame(dfEnergyRevenue_', :auto))
	for i in 1:G
		dfEnergyRevenue[!,:AnnualSum][i] = sum(dfEnergyRevenue[i,3:T+2])
	end
	dfEnergyRevenue_annualonly = dfEnergyRevenue[!,1:2]
	CSV.write(joinpath(path, "EnergyRevenue.csv"), dfEnergyRevenue_annualonly)
	return dfEnergyRevenue
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using Insolation, Documenter

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "100"

pages = Any[
        "Home" => "index.md"
        "Zenith Angle Equations" => "ZenithAngleEquations.md"
        "Insolation Examples" => "InsolationExamples.md"
        "Milankovitch Cycles" => "Milankovitch.md"
        "APIs" => "library.md"
]

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    collapselevel = 1,
)

makedocs(
    sitename = "Insolation.jl",
    format = format,
    clean = true,
    strict = true,
    modules = [Insolation],
    pages = pages,
)

deploydocs(
    repo = "github.com/CliMA/Insolation.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
function optcontract(ixs, xs, iy)
    size_dict = get_size_dict(ixs, xs)

    print("finding optimal tree ...")
    tree, cost = optimaltree(ixs, size_dict)
    @show tree, cost
    treecontract(tree, ixs, xs, iy)
end


text:
using ClimateMaps
using PyPlot
using Test


text:
include("discrete_laplacian.jl")

using Test

@testset "MALS linsolv" begin
    n = 10
    L = Lap(n,3)
    b = ones(n,n,n)
    x0 = randn(n,n,n)
    L_tt = tto_decomp(reshape(L,n,n,n,n,n,n))
    b_tt = ttv_decomp(b)
    x0_tt = ttv_decomp(x0)
    x = L\b[:]
    x_tt = mals_linsolv(L_tt,b_tt,x0_tt;tol=0.0)
    y = ttv_to_tensor(x_tt)
    @test isapprox(y[:],x)
end

@testset "MALS eigsolv" begin
    n = 10
    L = Lap(n,3)
    E = eigmin(L)
    x0 = randn(n,n,n)
    L_tt = tto_decomp(reshape(L,n,n,n,n,n,n))
    x0_tt = ttv_decomp(x0)
    E_tt,x_tt,r_hist = mals_eigsolv(L_tt,x0_tt;tol=0.0)
    y = ttv_to_tensor(x_tt)
    @test isapprox(E,E_tt[end])
    E_tt,x_tt,r_hist = mals_eigsolv(L_tt,x0_tt,it_solver=true)
    @test isapprox(E,E_tt[end],atol=1e-8)
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using H3.API
using H3.Lib
using Documenter

makedocs(
    build = joinpath(@__DIR__, "local" in ARGS ? "build_local" : "build"),
    modules = [API, Lib],
    clean = false,
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    sitename = "H3.jl ⬡",
    authors = "WooKyoung Noh",
    pages = Any[
        "Home" => "index.md",
        "API" => "API.md",
        "Lib" => "Lib.md",
    ],
)


text:
using VLConstraintBasedModelGenerationUtilities

# setup path to JSON file -
path_to_json_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/iJO1366.json"


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Unrestricted Open-Shell Hartree-Fock
#
# In the first two tutorials in this module, we wrote programs which implement a closed-shell formulation of Hartree-Fock theory using restricted orbitals, aptly named Restricted Hartree-Fock (RHF).  In this tutorial, we will abandon strictly closed-shell systems and the notion of restricted orbitals, in favor of a more general theory known as Unrestricted Hartree-Fock (UHF) which can accommodate more diverse molecules.  In UHF, the orbitals occupied by spin up ($\alpha$) electrons and those occupied by spin down ($\beta$) electrons no longer have the same spatial component, e.g.,
#
# $$\chi_i({\bf x}) = \begin{cases}\psi^{\alpha}_j({\bf r})\alpha(\omega) \\ \psi^{\beta}_j({\bf r})\beta(\omega)\end{cases},$$
#
# meaning that they will not have the same orbital energy.  This relaxation of orbital constraints allows for more variational flexibility, which leads to UHF always being able to find a lower total energy solution than RHF.
#
# ## I. Theoretical Overview
# In UHF, we seek to solve the coupled equations
#
# \begin{align}
# {\bf F}^{\alpha}{\bf C}^{\alpha} &= {\bf SC}^{\alpha}{\bf\epsilon}^{\alpha} \\
# {\bf F}^{\beta}{\bf C}^{\beta} &= {\bf SC}^{\beta}{\bf\epsilon}^{\beta},
# \end{align}
#
# which are the unrestricted generalizations of the restricted Roothan equations, called the Pople-Nesbet-Berthier equations.  Here, the one-electron Fock matrices are given by
#
# \begin{align}
# F_{\mu\nu}^{\alpha} &= H_{\mu\nu} + (\mu\,\nu\mid\lambda\,\sigma)[D_{\lambda\sigma}^{\alpha} + D_{\lambda\sigma}^{\beta}] - (\mu\,\lambda\,\mid\nu\,\sigma)D_{\lambda\sigma}^{\beta}\\
# F_{\mu\nu}^{\beta} &= H_{\mu\nu} + (\mu\,\nu\mid\,\lambda\,\sigma)[D_{\lambda\sigma}^{\alpha} + D_{\lambda\sigma}^{\beta}] - (\mu\,\lambda\,\mid\nu\,\sigma)D_{\lambda\sigma}^{\alpha},
# \end{align}
#
# where the density matrices $D_{\lambda\sigma}^{\alpha}$ and $D_{\lambda\sigma}^{\beta}$ are given by
#
# \begin{align}
# D_{\lambda\sigma}^{\alpha} &= C_{\sigma i}^{\alpha}C_{\lambda i}^{\alpha}\\
# D_{\lambda\sigma}^{\beta} &= C_{\sigma i}^{\beta}C_{\lambda i}^{\beta}.
# \end{align}
#
# Unlike for RHF, the orbital coefficient matrices ${\bf C}^{\alpha}$ and ${\bf C}^{\beta}$ are of dimension $M\times N^{\alpha}$ and $M\times N^{\beta}$, where $M$ is the number of AO basis functions and $N^{\alpha}$ ($N^{\beta}$) is the number of $\alpha$ ($\beta$) electrons.  The total UHF energy is given by
#
# \begin{align}
# E^{\rm UHF}_{\rm total} &= E^{\rm UHF}_{\rm elec} + E^{\rm BO}_{\rm nuc},\;\;{\rm with}\\
# E^{\rm UHF}_{\rm elec} &= \frac{1}{2}[({\bf D}^{\alpha} + {\bf D}^{\beta}){\bf H} +
# {\bf D}^{\alpha}{\bf F}^{\alpha} + {\bf D}^{\beta}{\bf F}^{\beta}].
# \end{align}
#
# ## II. Implementation
#
# In any SCF program, there will be several common elements which can be abstracted from the program itself into separate modules, classes, or functions to 'clean up' the code that will need to be written explicitly; examples of this concept can be seen throughout the Psi4Julia reference implementations.  For the purposes of this tutorial, we can achieve some degree of code cleanup without sacrificing readabilitiy and clarity by focusing on abstracting only the parts of the code which are both
# - Lengthy subroutines, and
# - Used repeatedly.
#
# In our UHF program, let's use what we've learned in the last tutorial by also implementing DIIS convergence accelleration for our SCF iterations.  With this in mind, two subroutines in particular would benefit from abstraction are
#
# 1. Orthogonalize & diagonalize Fock matrix
# 2. Extrapolate previous trial vectors for new DIIS solution vector
#
# Before we start writing our UHF program, let's try to write functions which can perform the above tasks so that we can use them in our implementation of UHF.  Recall that defining functions in Julia has the following syntax:
# ~~~julia
# function function_name(args; kwargs)
#     # function block
#     return_values
# end
# ~~~
# A thorough discussion of defining functions in Julia can be found [here](https://docs.julialang.org/en/v1/manual/functions/index.html "Go to Julia docs").  First, let's write a function which can diagonalize the Fock matrix and return the orbital coefficient matrix **C** and the density matrix **D**.  From our RHF tutorial, this subroutine is executed with:
# ~~~julia
# F_p =  A * F * A
# e, C_p = eigen(Hermitian(F_p))
# C = A * C_p
# C_occ = C[:, 1:ndocc]
# D = C_occ * C_occ'
# ~~~
# Examining this code block, there are three quantities which must be specified beforehand:
# - Fock matrix, **F**
# - Orthogonalization matrix, ${\bf A} = {\bf S}^{-1/2}$
# - Number of doubly occupied orbitals, `ndocc`
#
# However, since the orthogonalization matrix **A** is a static quantity (only built once, then left alone) we may choose to leave **A** as a *global* quantity, instead of an argument to our function.  In the cell below, using the code snippet given above, write a function `diag_F()` which takes **F** and the number of orbitals `norb` as arguments, and returns **C** and **D**:

# ==> Define function to diagonalize F <==
function diag_F(F, norb, A)
    F_p = A * F * A
    e, C_p = eigen(Hermitian(F_p))
    C = A * C_p
    C_occ = C[:, 1:norb]
    D = C_occ * C_occ'
    C, D
end

# Next, let's write a function to perform DIIS extrapolation and generate a new solution vector.  Recall that the DIIS-accellerated SCF algorithm is:
# #### Algorithm 1: DIIS within a generic SCF Iteration
# 1. Compute **F**, append to list of previous trial vectors
# 2. Compute AO orbital gradient **r**, append to list of previous residual vectors
# 3. Compute RHF energy
# 3. Check convergence criteria
#     - If RMSD of **r** sufficiently small, and
#     - If change in SCF energy sufficiently small, break
# 4. Build **B** matrix from previous AO gradient vectors
# 5. Solve Pulay equation for coefficients $\{c_i\}$
# 6. Compute DIIS solution vector **F_DIIS** from $\{c_i\}$ and previous trial vectors
# 7. Compute new orbital guess with **F_DIIS**
#
# In our function, we will perform steps 4-6 of the above algorithm.  What information will we need to provide our function in order to do so?  To build **B** (step 4 above) in the previous tutorial, we used:
# ~~~julia
# # Build B matrix
# B_dim = length(F_list) + 1
# B = zeros(B_dim, B_dim)
# B[end,   :] .= -1
# B[:  , end] .= -1
# B[end, end]  = 0
# for i in eachindex(F_list), j in eachindex(F_list)
#     B[i, j] = dot(DIIS_RESID[i], DIIS_RESID[j])
# end
# ~~~
# Here, we see that we must have all previous DIIS residual vectors (`DIIS_RESID`), as well as knowledge about how many previous trial vectors there are (for the dimension of **B**).  To solve the Pulay equation (step 5 above):
# ~~~julia
# # Build RHS of Pulay equation
# rhs = zeros(B_dim)
# rhs[end] = -1
#
# # Solve Pulay equation for c_i's with NumPy
# coeff = B \ rhs
# ~~~
# For this step, we only need the dimension of **B** (which we computed in step 4 above) and a Julia routine, so this step doesn't require any additional arguments.  Finally, to build the DIIS Fock matrix (step 6):
# ~~~julia
# # Build DIIS Fock matrix
# F = zeros(size(F_list[0]))
# for x in 1:length(coeff) - 1
#     F += coeff[x] * F_list[x]
# end
# ~~~
# Clearly, for this step, we need to know all the previous trial vectors (`F_list`) and the coefficients we generated in the previous step.  In the cell below, write a funciton `diis_xtrap()` according to Algorithm 1 steps 4-6, using the above code snippets, which takes a list of previous trial vectors `F_list` and residual vectors `DIIS_RESID` as arguments and returns the new DIIS solution vector `F_DIIS`:

# ==> Build DIIS Extrapolation Function <==
function diis_xtrap(F_list, DIIS_RESID)
    # Build B matrix
    B_dim = length(F_list) + 1
    B = zeros(B_dim, B_dim)
    B[end,   :] .= -1
    B[:  , end] .= -1
    B[end, end]  =  0
    for i in eachindex(F_list), j in eachindex(F_list)
       B[i, j] = dot(DIIS_RESID[i],  DIIS_RESID[j])
    end

    # Build RHS of Pulay equation
    rhs = zeros(B_dim)
    rhs[end] = -1

    # Solve Pulay equation for c_i's with Julia
    coeff = B \ rhs

    # Build DIIS Fock matrix
    F = zeros(size(F_list[1]))
    for i in 1:length(coeff) - 1
       F += coeff[i] * F_list[i]
    end
    F
end

# We are now ready to begin writing our UHF program!  Let's begin by importing  Psi4 , NumPy, TensorOperations, LinearAlgebra, and defining our molecule & basic options:

# ==> Import Psi4 & NumPy <==
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
using TensorOperations: @tensor
using LinearAlgebra: Diagonal, Hermitian, eigen, tr, norm, dot
using Printf: @printf

# +
# ==> Set Basic Psi4 Options <==
# Memory specification
psi4.set_memory(Int(5e8))
numpy_memory = 2

# Set output file
psi4.core.set_output_file("output.dat", false)

# Define Physicist's water -- don't forget C1 symmetry!
mol = psi4.geometry("""
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")

# Set computation options
psi4.set_options(Dict("basis"         => "cc-pvdz",
                      "scf_type"      => "pk",
                      "e_convergence" => 1e-8,
                      "guess"         => "core",
                      "reference"     => "uhf"))
# -

# You may notice that in the above `psi4.set_options()` block, there are two additional options -- namely, `'guess': 'core'` and `'reference': 'uhf'`.  These options make sure that when we ultimately check our program against  Psi4, the options  Psi4  uses are identical to our implementation.  Next, let's define the options for our UHF program; we can borrow these options from our RHF implementation with DIIS accelleration that we completed in our last tutorial.

# ==> Set default program options <==
# Maximum SCF iterations
MAXITER = 40
# Energy convergence criterion
E_conv = 1.0e-6
D_conv = 1.0e-3

# Static quantities like the ERI tensor, core Hamiltonian, and orthogonalization matrix have exactly the same form in UHF as in RHF.  Unlike in RHF, however, we will need the number of $\alpha$ and $\beta$ electrons.  Fortunately, both these values are available through querying the Wavefunction object.  In the cell below, generate these static objects and compute each of the following:
# - Number of basis functions, `nbf`
# - Number of alpha electrons, `nalpha`
# - Number of beta electrons, `nbeta`
# - Number of doubly occupied orbitals, `ndocc` (Hint: In UHF, there can be unpaired electrons!)

# +
# ==> Compute static 1e- and 2e- quantities with Psi4 <==
# Class instantiation
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option("basis"))
mints = psi4.core.MintsHelper(wfn.basisset())

# Overlap matrix
S = np.asarray(mints.ao_overlap()) # we only need a copy

# Number of basis Functions, alpha & beta orbitals, and # doubly occupied orbitals
nbf = wfn.nso()
nalpha = wfn.nalpha()
nbeta = wfn.nbeta()
ndocc = min(nalpha, nbeta)

println("Number of basis functions: ", nbf)
println("Number of singly occupied orbitals: ", abs(nalpha-nbeta))
println("Number of doubly occupied orbitals: ", ndocc)

# Memory check for ERI tensor
I_size = nbf^4 * 8.e-9
println("\nSize of the ERI tensor will be $I_size GB.")
memory_footprint = I_size * 1.5
if I_size > numpy_memory
    psi4.core.clean()
    throw(OutOfMemoryError("Estimated memory utilization ($memory_footprint GB) exceeds " *
                           "allotted memory limit of $numpy_memory GB."))
end

# Build ERI Tensor
I = np.asarray(mints.ao_eri()) # we only need a copy

# Build core Hamiltonian
T = np.asarray(mints.ao_kinetic()) # we only need a copy
V = np.asarray(mints.ao_potential()) # we only need a copy
H = T + V;

# Construct AO orthogonalization matrix A
A = mints.ao_overlap()
A.power(-0.5, 1.e-16) # ≈ Julia's A^(-0.5) after psi4view()
A = np.asarray(A);
# -

# Unlike the static quantities above, the CORE guess in UHF is slightly different than in RHF.  Since the $\alpha$ and $\beta$ electrons do not share spatial orbitals, we must construct a guess for *each* of the $\alpha$ and $\beta$ orbitals and densities.  In the cell below, using the function `diag_F()`, construct the CORE guesses and compute the nuclear repulsion energy:
#
# (Hint: The number of $\alpha$ orbitals is the same as the number of $\alpha$ electrons!)

# +
# ==> Build alpha & beta CORE guess <==
Ca, Da = diag_F(H, nalpha, A)
Cb, Db = diag_F(H, nbeta, A)

# Get nuclear repulsion energy
E_nuc = mol.nuclear_repulsion_energy()
# -

# We are almost ready to perform our SCF iterations; beforehand, however, we must initiate variables for the current & previous SCF energies, and the lists to hold previous residual vectors and trial vectors for the DIIS procedure.  Since, in UHF, there are Fock matrices ${\bf F}^{\alpha}$ and ${\bf F}^{\beta}$ for both $\alpha$ and $\beta$ orbitals, we must apply DIIS to each of these matrices separately.  In the cell below, define empty lists to hold previous Fock matrices and residual vectors for both $\alpha$ and $\beta$ orbitals:

# ==> Pre-Iteration Setup <==
# SCF & Previous Energy
SCF_E = 0.0
E_old = 0.0

# We are now ready to write the SCF iterations.  The algorithm for UHF-SCF iteration, with DIIS convergence accelleration, is:
# #### Algorithm 2: DIIS within UHF-SCF Iteration
# 1. Build ${\bf F}^{\alpha}$ and ${\bf F}^{\beta}$, append to trial vector lists
# 2. Compute the DIIS residual for $\alpha$ and $\beta$, append to residual vector lists
# 3. Compute UHF energy
# 4. Convergence check
#     - If average of RMSD of $\alpha$ and $\beta$ residual sufficiently small, and
#     - If change in UHF energy sufficiently small, break
# 5. DIIS extrapolation of ${\bf F}^{\alpha}$ and ${\bf F}^{\beta}$ to form new solution vector
# 6. Compute new ${\alpha}$ and ${\beta}$ orbital & density guesses
#
# In the cell below, write the UHF-SCF iteration according to Algorithm 2:
#
# (Hint: Use your functions `diis_xtrap()` and `diag_F` for Algorithm 2 steps 5 & 6, respectively)

SCF_E = let SCF_E = SCF_E, E_old = E_old, Da = Da, Db = Db, A = A, I = I, H = H, S = S

   # Trial & Residual Vector Lists -- one each for α & β
   F_list_a = []
   F_list_b = []
   R_list_a = []
   R_list_b = []

   # ==> UHF-SCF Iterations <==
   println("==> Starting SCF Iterations <==")

   # Begin Iterations
   for scf_iter in 1:MAXITER
       # Build Fa & Fb matrices
       @tensor Ja[p,q] := I[p,q,r,s] * Da[r,s]
       @tensor Jb[p,q] := I[p,q,r,s] * Db[r,s]
       @tensor Ka[p,q] := I[p,r,q,s] * Da[r,s]
       @tensor Kb[p,q] := I[p,r,q,s] * Db[r,s]
       Fa = H + (Ja + Jb) - Ka
       Fb = H + (Ja + Jb) - Kb

       # Compute DIIS residual for Fa & Fb
       diis_r_a = A * (Fa * Da * S - S * Da * Fa) * A
       diis_r_b = A * (Fb * Db * S - S * Db * Fb) * A

       # Append trial & residual vectors to lists
       push!(F_list_a, Fa)
       push!(F_list_b, Fb)
       push!(R_list_a, diis_r_a)
       push!(R_list_b, diis_r_b)

       # Compute UHF Energy
       SCF_E = 0.5*tr( H*(Da + Db) + Fa*Da + Fb*Db) + E_nuc

       dE = SCF_E - E_old
       dRMS = 0.5(norm(diis_r_a) + norm(diis_r_b))
       @printf("SCF Iteration %3d: Energy = %4.16f dE = %1.5e dRMS = %1.5e \n",
                           scf_iter,        SCF_E, SCF_E - E_old,     dRMS)

       # Convergence Check
       if abs(dE) < E_conv && dRMS < D_conv
           break
       end
       E_old = SCF_E

       # DIIS Extrapolation
       if scf_iter >= 2
           Fa = diis_xtrap(F_list_a, R_list_a)
           Fb = diis_xtrap(F_list_b, R_list_b)
       end

       # Compute new orbital guess
       Ca, Da = diag_F(Fa, nalpha, A)
       Cb, Db = diag_F(Fb, nbeta, A)

       # MAXITER exceeded?
       if scf_iter == MAXITER
           psi4.core.clean()
           throw(MethodError("Maximum number of SCF iterations exceeded."))
       end
   end
   SCF_E
end
# Post iterations
println("\nSCF converged.")
println("Final RHF Energy: $SCF_E [Eh]")
println()

# Congratulations! You've written your very own Unrestricted Hartree-Fock program with DIIS convergence accelleration!  Finally, let's check your final UHF energy against  Psi4:

# Compare to Psi4
SCF_E_psi = psi4.energy("SCF")
SCF_E
psi4.compare_values(SCF_E_psi, SCF_E, 6, "SCF Energy")

# ## References
# 1. A. Szabo and N. S. Ostlund, *Modern Quantum Chemistry*, Introduction to Advanced Electronic Structure Theory. Courier Corporation, 1996.
# 2. I. N. Levine, *Quantum Chemistry*. Prentice-Hall, New Jersey, 5th edition, 2000.
# 3. T. Helgaker, P. Jorgensen, and J. Olsen, *Molecular Electronic Structure Theory*, John Wiley & Sons Inc, 2000.


text:
# Run all tests.
module TestRunner

using FactCheck

include("test_util.jl")
include("test_loss.jl")
include("test_gb.jl")
include("test_gb_dt.jl")
include("test_gb_bl.jl")
include("test_ml.jl")
include("test_system.jl")

exitstatus()

end # module


text:
#***************************
#...LaGrange 3- and 4-point interpolation
#...arrays A and B are the npt data points,  given aa, a value of the
#...A variable, the routine will find the corresponding bb value
#
#...input:  aa
#...output: bb
using Random

function AtoB(aa,A,B,npt)
  for I in 2:(npt+1)

    if A[I-1] >= aa
      if I < 3 || I == npt
        J = I
        if I < 3
          J = 3
        end
        if I == npt
          J = npt
        end

        # J = J - 1
        A0D1 = A[J-2] - A[J-1]

        if A0D1 == 0.0
          A0D1=0.0001
        end

        A0D2=A[J-2]-A[J]

        if A0D2 == 0.0
          A0D2=0.0000
        end

        A1D1=A[J-1]-A[J-2]

        if A1D1 == 0.0
          A1D1=0.0001
        end

        A1D2=A[J-1]-A[J]

        if A1D2 == 0.0
          A1D2=0.0001
        end

        A2D1=A[J]-A[J-2]

        if A2D1 == 0.0
          A2D1=0.0001
        end

        A2D2=A[J]-A[J-1]

        if A2D2 == 0.0
          A2D2=0.0001
        end

        A0=(aa-A[J-1])*(aa-A[J])/(A0D1*A0D2)
        A1=(aa-A[J-2])*(aa-A[J])/(A1D1*A1D2)
        A2=(aa-A[J-2])*(aa-A[J-1])/(A2D1*A2D2)

        bb = A0*B[J-2] + A1*B[J-1] + A2*B[J]

        return bb

      else

        J = I
        # J = J-1   # zero index correction
        A0D1=A[J-2]-A[J-1]

        if A0D1 == 0.0
          A0D1=0.0001
        end
        A0D2=A[J-2]-A[J]
        if A0D2 == 0.0
          A0D2=0.0001
        end
        A0D3 = (A[J-2]-A[J+1])

        if A0D3 == 0.0
          A0D3=0.0001
        end
        A1D1=A[J-1]-A[J-2]
        if A1D1 == 0.0
          A1D1=0.0001
        end
        A1D2=A[J-1]-A[J]
        if A1D2 == 0.0
          A1D2=0.0001
        end
        A1D3 = A[J-1]-A[J+1]

        if A1D3 == 0.0
          A1D3=0.0001
        end

        A2D1=A[J]-A[J-2]
        if A2D1 == 0.0
          A2D1=0.0001
        end
        A2D2=A[J]-A[J-1]
        if A2D2 == 0.0
          A2D2=0.0001
        end
        A2D3 = A[J]-A[J+1]
        if A2D3 == 0.0
          A2D3=0.0001
        end

        A3D1 = A[J+1]-A[J-2]
        if A3D1 == 0.0
          A3D1=0.0001
        end
        A3D2 = A[J+1]-A[J-1]
        if A3D2 == 0.0
          A3D2=0.0001
        end
        A3D3 = A[J+1]-A[J]
        if A3D3 == 0.0
          A3D3=0.0001
        end

        A0=(aa-A[J-1])*(aa-A[J])*(aa-A[J+1])
        A0=A0/(A0D1*A0D2*A0D3)
        A1=(aa-A[J-2])*(aa-A[J])*(aa-A[J+1])
        A1=A1/(A1D1*A1D2*A1D3)
        A2=(aa-A[J-2])*(aa-A[J-1])*(aa-A[J+1])
        A2=A2/(A2D1*A2D2*A2D3)
        A3=(aa-A[J-2])*(aa-A[J-1])*(aa-A[J])
        A3=A3/(A3D1*A3D2*A3D3)

        bb = A0*B[J-2] + A1*B[J-1] + A2*B[J] + A3*B[J+1]

        return bb
      end
      break
    end
  end

  return bb

end

function qoft!(M, I, T,T_ref, result)

    # Get temperature grid
    TT = get_TT(M, I) #TT_[(M,I)]
    TQ = get_TQ(M, I)

    Tmin = minimum(TT); Tmax = maximum(TT)

    # Error if out of temperature range
    if T<Tmin || T>Tmax
      error("TIPS2017: T ($T) must be between $Tmin K and $Tmax K.")
    end

    # Interpolate partition sum for specified isotopologue
    # Lagrange interpolation from HAPI for now -- will eventually replace with
    # Interpolations.jl function
    # LI = LinearInterpolation(TT, TQ, extrapolation_bc=NaN)
    Qt = AtoB(T,TT,TQ,length(TT))
    Qt2 = AtoB(T_ref,TT,TQ,length(TT))

    # Save the ratio result
    result[1] = Qt2/Qt
    return nothing
end


text:
# Autogenerated wrapper script for libudfread_jll for armv7l-linux-gnueabihf
export libudfread

JLLWrappers.@generate_wrapper_header("libudfread")
JLLWrappers.@declare_library_product(libudfread, "libudfread.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libudfread,
        "lib/libudfread.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LinearAlgebra
using TensorOperations

#
# vhf tools
#
function get_coloumb(g, D)
    @tensor J[μ, ν] := g[μ, ν, ρ, σ] * D[ρ, σ]
    return J
end

function get_exchange(g, D)
    @tensor K[μ, ν] := g[μ, σ, ρ, ν] * D[ρ, σ]
    return K
end

function get_vhf(g, D)
    return get_coloumb(g, D) - 0.5 * get_exchange(g, D)
end

#
# RDM tools
#
function make_rdm1(mo_occ::Vector{Float64}, mo_coeff::Array{Float64})
    occ_coeff = mo_coeff[:, mo_occ.>0]
    return occ_coeff * diagm(0 => mo_occ[mo_occ.>0]) * occ_coeff'
end

#
# Fock tools
#
function get_fock(h, g, D)
    # This has to be Hermitian or eigenvalue solver with get messed up
    return Hermitian(h + get_vhf(g, D))
end

function solve_fock(h, g, D, s)
    res = eigen(get_fock(h, g, D), s)
    mo_energy = res.values
    mo_coeff = res.vectors
    return mo_energy, mo_coeff
end

function calc_energy_elec(h, g, dm)
    @tensor one_body = scalar(h[μ, ν] * dm[ν, μ])
    vhf = get_vhf(g, dm)
    @tensor two_body = 0.5 * scalar(vhf[μ, ν] * dm[ν, μ])
    return one_body + two_body, one_body, two_body
end


text:
module KeyChains

using Pkg.TOML
using DemeNet: Deme, DemeSpec, Signer, ID, DemeID, Profile, keydir, Certificate
using PeaceVote: AbstractVote, AbstractProposal, AbstractChain, load, attest, voters

import PeaceVote: braid!

struct KeyChain
    deme::Deme ### This is necessary to make braid! function obvious
    account::AbstractString
    member::Signer
    signers::Vector{Signer}
end

function KeyChain(deme::Deme,account::AbstractString)
    uuid = deme.spec.uuid
    notary = deme.notary

    ### Each app may have it's own subdirectory to preven conflicts.
    ### The member key could be the same!
    member = Signer(deme,account * "/member")
    voterdir = keydir(uuid) * account * "/voters/"
    mkpath(voterdir)

    t = []
    voters = Signer[]
    for file in readdir(voterdir)
        voter = Signer(uuid,notary,account * "/voters/$file")
        push!(t,mtime("$voterdir/$file"))
        push!(voters,voter)
    end

    if length(voters) > 0
        voters=voters[sortperm(t)]
    end

    return KeyChain(deme,account,member,voters)
end

### I could actually define a constructor for New so it would construct the type.

KeyChain(deme::Deme) = KeyChain(deme,"")

voter(kc::KeyChain) = kc.signers[end]

function voter(kc::KeyChain,vset::Set)
    for v in kc.signers
        if v.id in vset
            return v
        end
    end
end

Certificate(stuff::Union{AbstractProposal,ID},keychain::KeyChain) = Certificate(stuff,keychain.member)

### In place of this one I need to add Certificate methods for KeyChain
### The chain part could be omitted by storing the pid with the keys.
function Certificate(option::AbstractVote,chain::AbstractChain,keychain::KeyChain)
    loadedledger = load(chain)
    messages = attest(loadedledger,keychain.deme.notary)
    vset = voters(messages[1:option.pid])
    v = voter(keychain,vset)
    return Certificate(option,v)
end

function braid!(chain::AbstractChain,kc::KeyChain)
    config = chain.config

    if length(kc.signers)==0
        oldvoter = kc.member
    else
        oldvoter = kc.signers[end]
    end

    newvoter = Signer(kc.deme,kc.account * "/voters/$(string(oldvoter.id))")

    braid!(chain,newvoter,oldvoter)
    # if fails, delete the newvoter
    push!(kc.signers,newvoter)
end

end


text:
using Gadfly

set_default_plot_size(6inch, 6inch)

plot(x=rand(100), y=rand(100), shape=rand(["foo","bar","pooh"], 100), Geom.point)


text:
export RBF_SPD_regularization
function RBF_SPD_regularization(m::Array{Float64,1},mref,nRBF)

R = 0.0;
dR = zeros(size(m));
II = ones(Int64,nRBF*36);
JJ = ones(Int64,nRBF*36);
VV = zeros(Float64,nRBF*36);
I3 = Matrix(1.0I,3,3);
o6 = ones(Int64,6);
itrilA = idx_trilA();

ineqj = [2;3;5];

P = zeros(9,6);
P[itrilA,1:6] = Matrix(1.0I,6,6);
P[4,2] = 1.0; P[7,3] = 1.0; P[8,5] = 1.0;

for k=1:nRBF
	(a1,a2,a3,a4,a5,a6) = getA(k,m);
	A     = [a1 a2 a3 ; a2 a4 a5 ; a3 a5 a6];
	chol = 0;
	try
		chol = cholesky(A);
	catch
		R = Inf;
		break;
	end
	R      -= 2.0*log(det(chol.L));
	invA    = chol\I3;
	# invA = inv(A);
	# R -= log(det(A));
	iAiA 	= (kron(invA,invA)*P)[itrilA,:];
	offset  = 10*(k-1);
	AIDXs = (2+offset):(7+offset);
	invA = -invA[itrilA];
	invA[ineqj].*=2.0;
	iAiA[ineqj,:].*=2.0;
	dR[AIDXs] = invA;
	II[((k-1)*36+1) : (k*36)] = [o6*AIDXs[1];o6*AIDXs[2];o6*AIDXs[3];o6*AIDXs[4];o6*AIDXs[5];o6*AIDXs[6]];
	JJ[((k-1)*36+1) : (k*36)] = [AIDXs;AIDXs;AIDXs;AIDXs;AIDXs;AIDXs];
	VV[((k-1)*36+1) : (k*36)] = iAiA[:];
end
d2R = sparse(II,JJ,VV,length(m),length(m));
return R,dR,d2R;
end


text:
function train!(w, s, images, x, y, batchsizes, opts, o)
    values = []
    gloss = lossgradient(
        w, s, images, x, y, batchsizes; o=o, values=values)
    update!(w, gloss, opts)
    return values
end

function bulkloss(w, srnn, o, data, vocab)
    nsamples = length(data)
    ids = [1:nsamples...]
    nbatches = div(nsamples, o[:batchsize])

    total, nwords = 0.0, 0
    newo = Dict(:finetune => get(o, :finetune, false),
                :atype => get(o, :atype, AutoGrad.getval(typeof(w["wsoft"]))))
    for i = 1:nbatches
        lower = (i-1)*o[:batchsize]+1
        upper = min(lower+o[:batchsize]-1, nsamples)
        samples = data[ids[lower:upper]]
        images, x, y, batchsizes = make_batch(o, samples, vocab)
        values = []
        loss(w, srnn, images, x, y, batchsizes; o=newo, values=values)
        total  += values[1]
        nwords += values[2]
        images, captions, masks = 0, 0, 0; gc(); flush(STDOUT)
    end
    return total/nwords
end

# Training started (nsamples=30000, nbatches=120, loss=Inf, score=0). [2017-03-12T00:57:09.813]

# (epoch/iter): 1/120, loss: 4.64231/3.83807 [2017-03-12T00:58:30.048]
# BLEU = 60.5/33.7/18.9/11.5 (BP=0.941512, ratio=0.943157, hyp_len=8263, ref_len=8761) [2017-03-12T00:58:44.928]
# Model saved to /mnt/kufs/scratch/ikesen16/data/image-captioning/trained/nic-flickr8k-vgg16-v6-nofinetune-a-iter-120.jld

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
## load the libraries
# option load before IIF
# using Cairo, Fontfconfig, Gadfly

# the multimodal iSAM library
using IncrementalInference
using RoMEPlotting
# build some factor graph
fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0], Prior(Normal(0,1)))
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0, :x1], LinearRelative(Normal(10.0,1)))
addVariable!(fg, :x2, ContinuousScalar)
mmo = Mixture(LinearRelative, [Rayleigh(3); Uniform(30,55)], Categorical([0.4; 0.6]))
addFactor!(fg, [:x1, :x2], mmo)

# show the factor graph
drawGraph(fg, show=true)
# show the tree
tree = buildTreeReset!(fg, drawpdf=true, show=true)

# solve the factor graph and show solving progress on tree in src/JunctionTree.jl
fg.solverParams.showtree = true
fg.solverParams.drawtree = true
tree = solveTree!(fg)

## building a new tree -- as per IIF.prepBatchTree(...)

resetFactorGraphNewTree!(fg)

# Look at variable ordering used to build the Bayes net/tree
p = getEliminationOrder(fg, ordering=:qr)

fge = deepcopy(fg)

# Building Bayes net.
buildBayesNet!(fge, p)

# prep and build tree
tree = emptyBayesTree()
buildTree!(tree, fge, p)

# Find potential functions for each clique
cliq = tree.cliques[1] # start at the root
buildCliquePotentials(fg, tree, cliq);

drawTree(tree, show=true)

# println("Bayes Net")
# sleep(0.1)
#fid = open("bn.dot","w+")
#write(fid,to_dot(fge.bn))
#close(fid)

## can also show the Clique Association matrix by first importing Cairo, Fontconfig, Gadfly

cliq = tree.cliques[1]
cliq = getClique(tree, :x0) # where is :x0 a frontal variable
spyCliqMat(cliq)

tree = drawTree(tree, show=true, imgs=true)


text:
#define the critical parts of the interface
#note- we may need another layer if we make a version without an Intercept
abstract type AbstractVolumeParts{
    TM<:AbstractMatrix, TV<:AbstractVector, T<:Real} <: AbstractVector{T} end

#convenience assignment operator
function updatexfromy!(dest::T, source::T) where T<:AbstractVolumeParts
  dest .= source
  return nothing
end

#this is the predict function
function (Θ::AbstractVolumeParts)(X::NamedTuple, ::AbstractVector, ::Val{:level})
  @unpack ws, RLws = X
  #first map each coefficient to its appropriate
  A = view(Θ.A, Θ.expand.Aidx...)
  LA = view(Θ.A, Θ.expand.LAidx...)
  V₀ = view((exp).(Θ.V₀), Θ.expand.V₀idx)

  absμₖ = abs.(A .* ws' .- LA .* RLws')
  absμ = vec(sum(absμₖ,dims=1)) .+ V₀

  return absμ
end

#this is the predict function
function (Θ::AbstractVolumeParts)(X::NamedTuple, ::AbstractVector, ::Val{:levellog})
  @unpack ws, RLws = X
  #first map each coefficient to its appropriate
  A = view(exp.(Θ.A), Θ.expand.Aidx...)
  LA = view(exp.(Θ.A), Θ.expand.LAidx...)
  V₀ = view((exp).(Θ.V₀), Θ.expand.V₀idx)

  absμₖ = abs.(A .* ws' .- LA .* RLws')
  absμ = vec(sum(absμₖ,dims=1)) .+ V₀

  return absμ
end

genabsμₖ(Aᵢ, LAᵢ, wsᵢ, RLwsᵢ) = abs(Aᵢ * wsᵢ - LAᵢ * RLwsᵢ)
signabsμₖ(Aᵢ, LAᵢ, wsᵢ, RLwsᵢ) = sign(Aᵢ * wsᵢ - LAᵢ * RLwsᵢ)
cudagenabsμₖ(Aᵢ, LAᵢ, wsᵢ, RLwsᵢ) = CUDA.abs(Aᵢ * wsᵢ - LAᵢ * RLwsᵢ) #one of these should go away
cudaabsμₖ(Aᵢ, LAᵢ, wsᵢ, RLwsᵢ) = CUDA.abs(Aᵢ * wsᵢ - LAᵢ * RLwsᵢ)

#makes volumeparts callable for a level estimate
function (Θ::AbstractVolumeParts)(X::NamedTuple, ::CuVector{T}, ::Val{:level}) where T
  @unpack ws, RLws = X
  #re-creating the matrices gets around
  #https://github.com/JuliaGPU/CUDA.jl/issues/684
  #and https://github.com/FluxML/Zygote.jl/issues/600

  #REMINDER- A, LA, and V₀ are in logs here
  A = CuMatrix{T}(Θ.A[Θ.expand.Aidx...])
  LA = CuMatrix{T}(Θ.A[Θ.expand.LAidx...])
  V₀ = CUDA.exp.(Flux.CuArrays.CuVector{T}((Θ.V₀)[Θ.expand.V₀idx]))

  #absμₖ = CUDA.abs.(A .* ws' .- LA .* RLws')
  absμₖ = cudaabsμₖ.(A, LA, ws', RLws')
  absμ = CUDA.vec(CUDA.sum(absμₖ,dims=1))

  return absμ .+ V₀
end

#makes volumeparts callable for a level estimate
function (Θ::AbstractVolumeParts)(X::NamedTuple, ::CuVector{T}, ::Val{:levellog}) where T
  @unpack ws, RLws = X
  #re-creating the matrices gets around
  #https://github.com/JuliaGPU/CUDA.jl/issues/684
  #and https://github.com/FluxML/Zygote.jl/issues/600

  #REMINDER- A, LA, and V₀ are in logs here
  A = CUDA.exp.(CuMatrix{T}(Θ.A[Θ.expand.Aidx...]))
  LA = CUDA.exp.(CuMatrix{T}(Θ.A[Θ.expand.LAidx...]))
  V₀ = CUDA.exp.(CuVector{T}((Θ.V₀)[Θ.expand.V₀idx]))

  #absμₖ = CUDA.abs.(A .* ws' .- LA .* RLws')
  absμₖ = cudaabsμₖ.(A, LA, ws', RLws')
  absμ = CUDA.vec(CUDA.sum(absμₖ,dims=1)) .+ V₀

  return absμ
end

formresults(panel::AbstractDataFrame, ms::MeasureSpec, Θ::AbstractVolumeParts,
  ::Val{:levellog}) = (
  formresults(panel, ms, Θ, exp.(Θ.A)' .* grossleverages(ms)'))

formresults(panel::AbstractDataFrame, ms::MeasureSpec, Θ::AbstractVolumeParts,
  ::Val{:level}) = (
  formresults(panel, ms, Θ, (Θ.A |> Matrix!)' .* grossleverages(ms)'))

function formresults(panel::AbstractDataFrame, ms::MeasureSpec, Θ::AbstractVolumeParts,
    ::Val{:iter})

  Amat = (Θ.A |> Matrix!)' .* grossleverages(ms)'
  Gmat = vcat(zeros(eltype(Amat), size(Amat,2))', (Θ.G|>Matrix!)')
  return formresults(panel, ms, Θ, Amat, Gmat)
end

###The below provides 1d indexing of the parameters
Base.length(Θ::AbstractVolumeParts) = length(Θ.Π)
Base.size(Θ::AbstractVolumeParts) = (Base.length(Θ),)
Base.IndexStyle(::Type{<:AbstractVolumeParts}) = IndexLinear()

Base.getindex(Θ::AbstractVolumeParts, inds) = Θ.Π[inds]
Base.setindex!(Θ::AbstractVolumeParts, x::Real, i::Int) = (Θ.Π[i] = x)

∇rms(Θ::AbstractVolumeParts, ∇Θ) = sqrt(mean((vec∇(Θ, ∇Θ)) .^ 2))

####Expanders

levelexpandlsq(A,V₀, ts) = (
  Aidx = to_indices(A, (:, ts)),
  LAidx = to_indices(A, (:, ts .- 1)),
  V₀idx = ts .- 1)

levelexpandmle(A, V₀, ts) = (
  Aidx = to_indices(A, (:, ts)),
  LAidx = to_indices(A, (:, ts .- 1)),
  V₀idx = ts.-1,
  τidx = ts.-1)

volumepartsexpanders = Dict{Symbol, Function}(
  :levelmle => levelexpandmle,
  :levellsq => levelexpandlsq,
)

cudasignμₖ(Aᵢ, LAᵢ, wsᵢ, RLwsᵢ) = CUDA.sign(Aᵢ * wsᵢ - LAᵢ * RLwsᵢ)
cudanoop(x) = x

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Test
using Gumbo

include("basics.jl")
include("comparison.jl")
include("parsing.jl")
include("traversal.jl")
include("io.jl")


text:
"""
    AbstractMod{T}

Abstract type for all Modueles of type T

"""
abstract type AbstractMod{T} end

"""
    Mod{T} <: AbstractMod{T}

Structure to store one specific modules of type{T}

"""

# problem input as array in unique gives an error if dims = 2 not added... ne
struct Mod{T} <: AbstractMod{T}
    m::T
end

# some base base function for the new type
"""
    isless(mod1::Mod, mod2::Mod)

Evaluates the order of two modules based on the input characters,
needed to allow sorting.
How sorting is done isn't important for the model;
it only assures a specific rank between different modules.
For more information on how the different inputs are sorted exactly
[see](https://docs.julialang.org/en/v1/base/base/)
"""
isless(mod1::Mod, mod2::Mod) = isless(mod1.m,mod2.m)

"""
    isequal(mod1::Mod, mod2::Mod)

Evaluates if two modules are equal
for more infromation [see](https://docs.julialang.org/en/v1/base/base/)
"""
isequal(mod1::Mod, mod2::Mod) = isequal(mod1.m,mod2.m)

"""
    length(mod1::Mod)

length of module, is 1
"""

Base.length(mod1::Mod) = 1

"""
    GroupMod{N <: Mod{T} where T}

Structure to group multiple modules.
The input values are first filtered to prevent duplicated modules.
Afterwards, the modules are sorted to give consistent results even if modules are
oad in a different order.
"""
struct GroupMod{N <: Mod{T} where T} <: AbstractMod{N}
    m::Array{N}
    GroupMod(m) =  m |> Set |> collect |> sort |> (y -> new{eltype(m)}(y))
end

"""
    group_mod(input::Array{T} where T)

A Function to facilitate the input of multiple modules. It returns a "Group_mod"
structure. The input is an array containing the data that needs to be transformed
to an ``Mod`` and grouped afterwards. See [`Group_mod`](@ref)

"""
groupmod(input::Array{T} where T) = GroupMod([Mod(newmod) for newmod in input])

"""
    iterate(Group::GroupMod,state =1 )

Extensions of ``Base.iterate`` for modules of type ``Group_Mod``
for more infromation [see](https://docs.julialang.org/en/v1/manual/interfaces/)
"""
function Base.iterate(Group::GroupMod,state =1 )
    if  state <= length(Group.m)
        mod = Group.m[state]
        state += 1
        return (mod,state)
    else
        return
    end
end

"""
    length(m::GroupMod)

Returns the number of `Mod` in `m`.
"""
Base.length(m::GroupMod) = length(m.m)

"""
    eltype(m::GroupMod)

Returns the type of the iterator.
"""
Base.eltype(::GroupMod{T}) where {T} = T

"""

"""
Base.getindex(m::GroupMod,i) = m.m[i]

function Base.show(io::IO,m::GroupMod)
    s = [mi.m for mi in m.m]
    s = replace(string(s), "[" => "{") |> x->  replace(x, "]" => "}")
    print(io,"{GroupMod}$s")
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


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

using Markdown
using InteractiveUtils

# ╔═╡ cd892ad4-fdba-11ea-2b87-41c0d68c6a27
using Pkg, DrWatson

# ╔═╡ cd89683c-fdba-11ea-0ff4-3d500f086774
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ f2c5b4bc-fdb9-11ea-1f2f-13802209834b
md"## Clip-05-41s.jl"

# ╔═╡ cd89eb2e-fdba-11ea-3689-f933b28ab670
begin
	df1 = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
	df1 = filter(row -> !(row[:neocortex_perc] == "NA"), df1);

	df = DataFrame()
	df[!, :NC] = parse.(Float64, df1[:, :neocortex_perc])
	df[!, :M] = log.(df1[:, :mass])
	df[!, :K] = df1[:, :kcal_per_g]
	scale!(df, [:K, :NC, :M])
end;

# ╔═╡ 2f8bd97c-fdbe-11ea-0855-55a321d0e010
# Define the Stan language model

stan5_7_A = "
data {
  int N;
  vector[N] K;
  vector[N] M;
  vector[N] NC;
}
parameters {
  real a;
  real bN;
  real bM;
  real aNC;
  real bMNC;
  real<lower=0> sigma;
  real<lower=0> sigma_NC;
}
model {
  // M -> K <- NC
  vector[N] mu = a + bN * NC + bM * M;
  a ~ normal( 0 , 0.2 );
  bN ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  K ~ normal( mu , sigma );
  // M -> NC
  vector[N] mu_NC = aNC + bMNC * M;
  aNC ~ normal( 0 , 0.2 );
  bMNC ~ normal( 0 , 0.5 );
  sigma_NC ~ exponential( 1 );
  NC ~ normal( mu_NC , sigma_NC );
}
";

# ╔═╡ 664db714-fdbe-11ea-20ba-b18074f27490
begin
	m5_7_As = SampleModel("m5.7_A", stan5_7_A);
	m5_7_A_data = Dict(
	  "N" => size(df, 1),
	  "K" => df[:, :K_s],
	  "M" => df[:, :M_s],
	  "NC" => df[:, :NC_s]
	);
	rc5_7_As = stan_sample(m5_7_As, data=m5_7_A_data);
	post5_7_As_df = read_samples(m5_7_As, :dataframe);
end;

# ╔═╡ cd977756-fdba-11ea-0b92-79d3b2439454
md"### Snippet 5.22"

# ╔═╡ cda28e98-fdba-11ea-2e23-95ab5e6d382a
a_seq = range(-2, stop=2, length=100)

# ╔═╡ cda9c6f4-fdba-11ea-359b-a5464143c8a6
m_sim, d_sim = simulate(post5_7_As_df, [:aNC, :bMNC, :sigma_NC], a_seq, [:bM, :sigma]);

# ╔═╡ cdab170c-fdba-11ea-1c9b-17e26c5d30f5
md"### Snippet 5.24"

# ╔═╡ cdb7171e-fdba-11ea-350c-f3cc664af8d9
begin
	plot(xlab="Manipulated M", ylab="Counterfactual K",
		title="Total counterfactual effect of M on K")
	plot!(a_seq, mean(d_sim, dims=1)[1, :], leg=false)
	hpdi_array = zeros(length(a_seq), 2)
	for i in 1:length(a_seq)
		hpdi_array[i, :] =  hpdi(d_sim[i, :])
	end
	plot!(a_seq, mean(d_sim, dims=1)[1, :]; ribbon=(hpdi_array[:, 1], -hpdi_array[:, 2]))
end

# ╔═╡ cdb7b26e-fdba-11ea-3abb-2714f92bb8c8
md"## End of clip-05-41s.jl"

# ╔═╡ Cell order:
# ╠═f2c5b4bc-fdb9-11ea-1f2f-13802209834b
# ╠═cd892ad4-fdba-11ea-2b87-41c0d68c6a27
# ╠═cd89683c-fdba-11ea-0ff4-3d500f086774
# ╠═cd89eb2e-fdba-11ea-3689-f933b28ab670
# ╠═2f8bd97c-fdbe-11ea-0855-55a321d0e010
# ╠═664db714-fdbe-11ea-20ba-b18074f27490
# ╟─cd977756-fdba-11ea-0b92-79d3b2439454
# ╠═cda28e98-fdba-11ea-2e23-95ab5e6d382a
# ╠═cda9c6f4-fdba-11ea-359b-a5464143c8a6
# ╟─cdab170c-fdba-11ea-1c9b-17e26c5d30f5
# ╠═cdb7171e-fdba-11ea-350c-f3cc664af8d9
# ╟─cdb7b26e-fdba-11ea-3abb-2714f92bb8c8


text:
module GardenOptim

using Logging

using DocStringExtensions
using Tables

export loadplants, loadgarden, loadclassification, loadcosts
export update!, fillgardenrandomly!, gardenevolution!, outputgarden

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(TYPEDSIGNATURES)
    $(DOCSTRING)
    """

include("classification.jl")
include("loaddata.jl")
include("mcmc.jl")

"Save the garden to a CSV file."
function outputgarden(garden::Matrix{Int}, plants::Vector{Symbol})
    output = vcat([""], String.(plants))[garden .+ 1]
    CSV.write("out/output.csv", Tables.table(output), writeheader=false)
end

end # module

max_stars_repo_name: dlozeve/GardenOptim

max_issues_repo_name: dlozeve/GardenOptim

max_forks_repo_name: dlozeve/GardenOptim


text:
function ChainRulesCore.rrule(::Type{<:ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11
                }},u,args...) where {T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}

                function ODESolutionAdjoint(ȳ)
                  (NoTangent(),ȳ,ntuple(_->NoTangent(), length(args))...)
                end

                ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}(u,args...),ODESolutionAdjoint
end

ZygoteRules.@adjoint function ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11
                }(u,args...) where {T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}

                function ODESolutionAdjoint(ȳ)
                  (ȳ,ntuple(_->nothing, length(args))...)
                end

                ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}(u,args...),ODESolutionAdjoint
end

function ChainRulesCore.rrule(::Type{<:ODESolution{uType,tType,isinplace,P,NP,F,G,K,ND
                }},u,args...) where {uType,tType,isinplace,P,NP,F,G,K,ND}

                function SDESolutionAdjoint(ȳ)
                  (NoTangent(),ȳ,ntuple(_->NoTangent(), length(args))...)
                end

                SDESolution{uType,tType,isinplace,P,NP,F,G,K,ND}(u,args...),SDESolutionAdjoint
end

ZygoteRules.@adjoint function SDEProblem{uType,tType,isinplace,P,NP,F,G,K,ND}(u,args...) where
                              {uType,tType,isinplace,P,NP,F,G,K,ND}

                function SDESolutionAdjoint(ȳ)
                  (ȳ,ntuple(_->nothing, length(args))...)
                end

                SDESolution{uType,tType,isinplace,P,NP,F,G,K,ND}(u,args...),SDESolutionAdjoint
end

ZygoteRules.@adjoint function NonlinearSolution{T,N,uType,R,P,A,O,uType2}(u,
                                      args...) where {T,N,uType,R,P,A,O,uType2}
    function NonlinearSolutionAdjoint(ȳ)
      (ȳ,ntuple(_->nothing, length(args))...)
    end
    NonlinearSolution{T,N,uType,R,P,A,O,uType2}(u,args...),NonlinearSolutionAdjoint
end

ZygoteRules.@adjoint function ZygoteRules.literal_getproperty(sol::AbstractTimeseriesSolution, ::Val{:u})
  function solu_adjoint(Δ)
        zerou = zero(sol.prob.u0)
        _Δ = @. ifelse(Δ == nothing,(zerou,),Δ)
        (DiffEqBase.build_solution(sol.prob,sol.alg,sol.t,_Δ),)
  end
  sol.u,solu_adjoint
end

ZygoteRules.@adjoint function ZygoteRules.literal_getproperty(sol::AbstractNoTimeSolution, ::Val{:u})
  function solu_adjoint(Δ)
        zerou = zero(sol.prob.u0)
        _Δ = @. ifelse(Δ == nothing,zerou,Δ)
        (DiffEqBase.build_solution(sol.prob,sol.alg,_Δ,sol.resid),)
  end
  sol.u,solu_adjoint
end

function ChainRulesCore.rrule(::DiffEqBase.EnsembleSolution,sim,time,converged)
  out = EnsembleSolution(sim,time,converged)
  function EnsembleSolution_adjoint(p̄::AbstractArray{T,N}) where {T,N}
    arrarr = [[p̄[ntuple(x->Colon(),Val(N-2))...,j,i] for j in 1:size(p̄)[end-1]] for i in 1:size(p̄)[end]]
    (NoTangent(),EnsembleSolution(arrarr, 0.0, true),NoTangent(),NoTangent())
  end
  function EnsembleSolution_adjoint(p̄::EnsembleSolution)
    (NoTangent(),p̄,NoTangent(),NoTangent())
  end
  out,EnsembleSolution_adjoint
end

ZygoteRules.@adjoint function DiffEqBase.EnsembleSolution(sim,time,converged)
  out = EnsembleSolution(sim,time,converged)
  function EnsembleSolution_adjoint(p̄::AbstractArray{T,N}) where {T,N}
    arrarr = [[p̄[ntuple(x->Colon(),Val(N-2))...,j,i] for j in 1:size(p̄)[end-1]] for i in 1:size(p̄)[end]]
    (EnsembleSolution(arrarr, 0.0, true),nothing,nothing)
  end
  function EnsembleSolution_adjoint(p̄::AbstractArray{<:AbstractArray,1})where {T}
    (EnsembleSolution(p̄, 0.0, true),nothing,nothing)
  end
  function EnsembleSolution_adjoint(p̄::EnsembleSolution)
    (p̄,nothing,nothing)
  end
  out,EnsembleSolution_adjoint
end

ZygoteRules.@adjoint ZygoteRules.literal_getproperty(sim::EnsembleSolution, ::Val{:u}) = sim.u, p̄ -> (EnsembleSolution(p̄, 0.0, true),)

#=
ChainRulesCore.frule(f::ODEFunction,u,p,t)
  if f.jvp === nothing
    ChainRulesCore.frule(f.f,u,p,t)
  else
    function ode_jvp(f,du,dp,dt)
      f.jvp_u(du,u,p,t) + f.jvp_p(dp,u,p,t) + f.tgrad(u,p,t)*dt
    end
    f.f(u,p,t),ode_jvp
end
=#

function ChainRulesCore.rrule(f::ODEFunction,u,p,t)
  if f.vjp === nothing
    ChainRulesCore.rrule(f.f,u,p,t)
  else
    f.vjp(u,p,t)
  end
end

ZygoteRules.@adjoint function (f::ODEFunction)(u,p,t)
  if f.vjp === nothing
    ZygoteRules._pullback(f.f,u,p,t)
  else
    f.vjp(u,p,t)
  end
end

#=
ZygoteRules.@adjoint! function (f::ODEFunction)(du,u,p,t)
  if f.vjp === nothing
    ZygoteRules._pullback(f.f,du,u,p,t)
  else
    f.vjp(du,u,p,t)
  end
end
=#

ChainRulesCore.rrule(::typeof(numargs),f) = (numargs(f),df->(NoTangent(),NoTangent()))


text:
function get_pages(resp)

    headers = Dict(resp.headers)

    before, after = parse.(Int64, (get(headers, "Cb-Before", ""), get(headers, "Cb-After", "")))
    (before, (before+after)/2, after)
end


text:
using DashBootstrapComponents, DashHtmlComponents

items = [
    dbc_dropdownmenuitem("Item 1"),
    dbc_dropdownmenuitem("Item 2"),
    dbc_dropdownmenuitem("Item 3"),
];

dropdowns = html_div(
    [
        dbc_dropdownmenu(items, label="Primary", color="primary", className="m-1"),
        dbc_dropdownmenu(
            items,
            label="Secondary",
            color="secondary",
            className="m-1",
        ),
        dbc_dropdownmenu(items, label="Success", color="success", className="m-1"),
        dbc_dropdownmenu(items, label="Warning", color="warning", className="m-1"),
        dbc_dropdownmenu(items, label="Danger", color="danger", className="m-1"),
        dbc_dropdownmenu(items, label="Info", color="info", className="m-1"),
        dbc_dropdownmenu(items, label="Link", color="link", className="m-1"),
    ],
    style=Dict("display" => "flex", "flexWrap" => "wrap"),
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# error handling

export KernelError, InternalCompilerError

struct KernelError <: Exception
    job::CompilerJob
    message::String
    help::Union{Nothing,String}
    bt::StackTraces.StackTrace

    KernelError(@nospecialize(job::CompilerJob), message::String, help=nothing;
                bt=StackTraces.StackTrace()) =
        new(job, message, help, bt)
end

function Base.showerror(io::IO, err::KernelError)
    println(io, "GPU compilation of ", err.job.source, " failed")
    println(io, "KernelError: $(err.message)")
    println(io)
    println(io, something(err.help, "Try inspecting the generated code with any of the @device_code_... macros."))
    Base.show_backtrace(io, err.bt)
end

struct InternalCompilerError <: Exception
    job::CompilerJob
    message::String
    meta::Dict
    InternalCompilerError(job, message; kwargs...) = new(job, message, kwargs)
end

function Base.showerror(io::IO, err::InternalCompilerError)
    println(io, """GPUCompiler.jl encountered an unexpected internal error.
                   Please file an issue attaching the following information, including the backtrace,
                   as well as a reproducible example (if possible).""")

    println(io, "\nInternalCompilerError: $(err.message)")

    println(io, "\nCompiler invocation: ", err.job)

    if !isempty(err.meta)
        println(io, "\nAdditional information:")
        for (key,val) in err.meta
            println(io, " - $key = $(repr(val))")
        end
    end

    let Pkg = Base.require(Base.PkgId(Base.UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f"), "Pkg"))
        println(io, "\nInstalled packages:")
        for (uuid, pkg) in Pkg.dependencies()
            println(io, " - $(pkg.name) = $(repr(pkg.version))")
        end
    end

    println(io)

    let InteractiveUtils = Base.require(Base.PkgId(Base.UUID("b77e0a4c-d291-57a0-90e8-8db25a27a240"), "InteractiveUtils"))
        InteractiveUtils.versioninfo(io)
    end
end

macro compiler_assert(ex, job, kwargs...)
    msg = "$ex, at $(__source__.file):$(__source__.line)"
    return :($(esc(ex)) ? $(nothing)
                        : throw(InternalCompilerError($(esc(job)), $msg;
                                                      $(map(esc, kwargs)...)))
            )
end


text:
import Base: mod, real, imag, reim, conj, promote_rule

export GaussMod, AbstractMod

# mod support for Gaussian integers until officially adopted into Base
mod(z::Complex{<:Integer}, n::Integer) = Complex(mod(real(z), n), mod(imag(z), n))

"""
`GaussMod{N,T}` is an alias of `Mod{N,Complex{T}}`.
It is for computing Gaussian Modulus.
"""
const GaussMod{N,T} = Mod{N,Complex{T}}
GaussMod{N}(x::T) where {N,T<:Integer} = Mod{N,Complex{T}}(x)
GaussMod{N}(x::T) where {N,T<:Complex} = Mod{N,T}(x)
Mod{N}(x::Complex{Rational{T}}) where {N,T} = Mod{N}(real(x)) + Mod{N}(imag(x)) * im
Mod{N}(re::Integer, im::Integer) where {N} = Mod{N}(Complex(re, im))

reim(x::AbstractMod) = (real(x), imag(x))
real(x::Mod{N}) where {N} = Mod{N}(real(x.val))
imag(x::Mod{N}) where {N} = Mod{N}(imag(x.val))
conj(x::Mod{N}) where {N} = Mod{N}(conj(x.val))

# ARITHMETIC
function (+)(x::GaussMod{N,T}, y::GaussMod{N,T}) where {N,T}
    xx = widen(x.val)
    yy = widen(y.val)
    zz = mod(xx + yy, N)
    return GaussMod{N,T}(zz)
end

(-)(x::GaussMod{N}) where {N} = Mod{N}(-x.val)

function (*)(x::GaussMod{N,T}, y::GaussMod{N,T}) where {N,T}
    xx = widen(x.val)
    yy = widen(y.val)
    zz = mod(xx * yy, N)
    return GaussMod{N,T}(zz)
end

function inv(x::GaussMod{N}) where {N}
    try
        a, b = reim(x)
        bot = inv(a * a + b * b)
        aa = a * bot
        bb = -b * bot
        return aa + bb * im
    catch
        error("$x is not invertible")
    end
end

is_invertible(x::GaussMod) = is_invertible(real(x * x'))
rand(::Type{GaussMod{N}}, args::Integer...) where {N} = rand(GaussMod{N,Int}, args...)