text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function genbo3(foldername, filename, u, xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, rho, mu, eta, kappa, natomtype)

# n,m = size(u);
# if (u[1]==0)
#     m = 0;
# else
#     n,m = size(u);
#     if (n==0) || (n == 1)
#         m = 0;
#     elseif (n == natomtype+2)
#         if u[n,1] == 0
#             m = 0;
#         end
#     else
#         error("Number of functions for the two-body bond order potential must be equal to 1 plus the number of atom types");
#     end
# end

if length(u)>0
    if (u[1]==0)
        m = 0;
        n = 0;
    else
        n,m = size(u);
    end
else
    n = 0;
    m = 0;
end

gen = 1;
ifile = 0;
stropu = "";
strcpu = "";
strgpu = "";
potnum = 0;
fp = "TripletcPair";
fn = "TripletcDensity";

for i = 1:m
    ui = u[:,i];
    if ui[n] != 0
        potnum = potnum + 1;
        stropu0, strcpu0, strgpu0 = gentriplet(foldername, filename * string(potnum), ui[1:(n-2)], xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, mu, eta, kappa, gen, ifile);
        stropu0 = modifystr3(stropu0, natomtype);
        strcpu0 = modifystr3(strcpu0, natomtype);
        strgpu0 = modifystr3(strgpu0, natomtype);
        stropu1, strcpu1, strgpu1 = genpair(foldername, fp * string(potnum), ui[n-1], xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
        stropu2, strcpu2, strgpu2 = gendensity(foldername, fn * string(potnum), ui[end], rho, mu, eta, kappa, gen, ifile);
        stropu = stropu * stropu0 * "\n" * stropu1 * "\n" * stropu2 * "\n";
        strcpu = strcpu * strcpu0 * "\n" * strcpu1 * "\n" * strcpu2 * "\n";
        strgpu = strgpu * strgpu0 * "\n" * strgpu1 * "\n" * strgpu2 * "\n";
    end
end

if potnum==0
    ifile = 0;
    gen = 0;
    stropu, strcpu, strgpu = gentriplet(foldername, filename, u, xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, mu, eta, kappa, gen, ifile);
    tmqopu, tmqcpu, tmqgpu = genpair(foldername, fp, u, xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
    tmpopu, tmpcpu, tmpgpu = gendensity(foldername, fn, u, mu, mu, eta, kappa, gen, ifile);
    stropu = stropu  * "\n" * tmqopu * "\n" * tmpopu * "\n";
    strcpu = strcpu  * "\n" * tmqcpu * "\n" * tmpcpu * "\n";
    strgpu = strgpu  * "\n" * tmqgpu * "\n" * tmpgpu * "\n";
    cppfiles(foldername, filename, stropu, strcpu, strgpu, 1);
else
    sp0, sp1, sp2, sp3 = getpotstr(3);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(filename, potnum, sp0, sp1, sp2, sp3);
    stropu = stropu * tmpopu * "\n";
    strcpu = strcpu * tmpcpu * "\n";
    strgpu = strgpu * tmpgpu * "\n";

    # here
    tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_xij, T *u_xik,");
    tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_xij, u_xik,");
    tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*, double*");
    tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*, float*");
    strgpu = strgpu * "\n" * tmpgpu;

    sp0, sp1, sp2, sp3 = getpotstr(2);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(fp, potnum, sp0, sp1, sp2, sp3);
    stropu = stropu * tmpopu * "\n";
    strcpu = strcpu * tmpcpu * "\n";
    strgpu = strgpu * tmpgpu * "\n";

    # here
    tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_xij,");
    tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_xij,");
    tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*");
    tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*");
    strgpu = strgpu * "\n" * tmpgpu;

    sp0, sp1, sp2, sp3 = getpotstr(0);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(fn, potnum, sp0, sp1, sp2, sp3);
    stropu = stropu * tmpopu;
    strcpu = strcpu * tmpcpu;
    strgpu = strgpu * tmpgpu;

    # here
    tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_rho,");
    tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_rho,");
    tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*");
    tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*");
    strgpu = strgpu * "\n" * tmpgpu;

    cppfiles(foldername, filename, stropu, strcpu, strgpu, 0);
end

end

function modifystr3(oldstr, natomtype)

newstr = oldstr;
for i = 1:natomtype
    s = "u[" * string(i-1) * " * i*" * string(natomtype) * "]";
    if contains(newstr, s)
        news = "if (tk[i] == " * string(i) * ") \n\t\t\t" * "u[i]";
        newstr = replace(newstr, s, news);
    end
end

return newstr

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module TestSnap
using Test
import GraphIdx.Io.Snap:
  parse_snap_idx, parse_snap_header, parse_snap_edges
import GraphIdx

const TEST_SNAP = """# comments
# ... to be ignored
#
# Nodes: 4 Edges: 3
# FromNodeId    ToNodeId
1       2
1       4
3       4
"""

@testset "SNAP reader                    " begin
    idx = parse_snap_idx(IOBuffer(TEST_SNAP))
    n, m, head, tail = parse_snap_edges(IOBuffer(TEST_SNAP), Int)
    @test parse_snap_header(IOBuffer(TEST_SNAP)) == (4, 3)
    @test GraphIdx.num_nodes(idx) == 4
    @test GraphIdx.num_edges(idx) == 3
    @test n == 4
    @test m == 3
    @test head == [1, 1, 3]
    @test tail == [2, 4, 4]
    @test idx == GraphIdx.BiAdjacentIndex([0, 2, 3, 4, 6] .+ 1,
                                          [1, 3, 0, 3, 0, 2] .+ 1)
end

end

ext: jl

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
droplast(x::Tuple) = (x[1:(end - 1)]...,)

function convolution(
    _src::Memory{T},
    _weights::Memory{T},
    bias::Memory,
    dims::Dims{N};
    kind = Inference(),
    algo = Lib.dnnl_convolution_auto,
    attributes = noattributes(),
    src_md = memorydesc_ptr(_src),
    weights_md = memorydesc_ptr(_weights),
    callback = (x...) -> nothing,
    opdesc = Ref{Lib.dnnl_convolution_desc_t}(),
) where {T,N}
    dst_dims = tuple_replace(output_size(size(_src), dims), size(_weights, 4), Val(3))
    dst_md = memorydesc(T, dst_dims, dnnl_format_any())

    @apicall dnnl_convolution_forward_desc_init(
        opdesc,
        kind,
        algo,
        src_md,
        weights_md,
        bias,
        dst_md,
        dims.strides,
        dims.padding,
        dims.padding,
    )

    return temp_primitive(opdesc, attributes, global_engine(), noforward()) do p, pd
        callback(p, pd)

        # Convert if needed.
        src = reorder_if_any(src_md, pd, _src, @query(src))
        weights = reorder_if_any(weights_md, pd, _weights, @query(weights))

        # Allocate destination
        dst_md = query_md(pd, @query(dst))
        dst = similar(src, T, dst_dims, dst_md)

        # Run kernel
        execute!(p, @dnnl_args src weights bias dst)
        if kind === Inference()
            return dst
        elseif kind === Training()
            return (; dst, forward = pd)
        end
    end
end

function convolution_backward_data(
    diff_src_dims::NTuple,
    _weights::Memory{T},
    _diff_dst::Memory{T},
    dims::Dims;
    forward,
    algo = Lib.dnnl_convolution_auto,
    opdesc = Ref{Lib.dnnl_convolution_desc_t}(),
) where {T}
    diff_src_md = memorydesc(T, diff_src_dims, dnnl_format_any())

    @apicall dnnl_convolution_backward_data_desc_init(
        opdesc,
        algo,
        diff_src_md,
        toany(_weights),
        toany(_diff_dst),
        dims.strides,
        dims.padding,
        dims.padding,
    )

    return temp_primitive(opdesc, noattributes(), global_engine(), forward) do p, pd
        weights = maybe_reorder(pd, _weights, @query(weights))
        diff_dst = maybe_reorder(pd, _diff_dst, @query(diff_dst))

        diff_src_md = query_md(pd, @query(diff_src))
        diff_src = similar(diff_dst, T, diff_src_dims, diff_src_md)
        execute!(p, @dnnl_args diff_src weights diff_dst)
        # Optimization - return the reformatted `diff_dst` for use in the
        # backward_weights kernel.
        return (; diff_src, diff_dst)
    end
end

function convolution_backward_weights(
    diff_weights_dims::NTuple{N},
    _src::Memory{T},
    _diff_dst::Memory{T},
    dims::Dims;
    forward,
    algo = Lib.dnnl_convolution_auto,
    opdesc = Ref{Lib.dnnl_convolution_desc_t}(),
) where {T,N}
    diff_bias_dims = (diff_weights_dims[N],)
    diff_weights_md = memorydesc(T, diff_weights_dims, dnnl_format_any())
    diff_bias_md = memorydesc(T, diff_bias_dims)

    @apicall dnnl_convolution_backward_weights_desc_init(
        opdesc,
        algo,
        toany(_src),
        diff_weights_md,
        diff_bias_md,
        toany(_diff_dst),
        dims.strides,
        dims.padding,
        dims.padding,
    )

    return temp_primitive(opdesc, noattributes(), global_engine(), forward) do p, pd
        # Maybe convert arguments
        src = maybe_reorder(pd, _src, @query(src))
        diff_dst = maybe_reorder(pd, _diff_dst, @query(diff_dst))

        # Materialize destinations
        diff_weights_md = query_md(pd, @query(diff_weights))
        diff_weights = similar(diff_dst, T, diff_weights_dims, diff_weights_md)
        diff_bias = similar(diff_dst, T, diff_bias_dims, diff_bias_md)

        execute!(p, @dnnl_args diff_weights diff_bias src diff_dst)
        return (diff_weights, diff_bias)
    end
end

## Convolutional layer
mutable struct Conv{W<:Memory,B<:Memory,F,N}
    weights::W
    bias::B
    activation::F
    dims::Dims{N}
    attributes::Attributes
    optimized_weights::Bool
    ### Op Descriptors
    opdesc::Base.RefValue{Lib.dnnl_convolution_desc_t}
end

function Conv(
    _weights::AbstractArray{T,N},
    _bias::AbstractVector{T},
    activation = identity;
    stride = 1,
    padding = 0,
    dilation = 0,
) where {T,N}
    kernel = (size(_weights, 1), size(_weights, 2))
    stride = expand(Val(2), stride)
    dilation = expand(Val(2), dilation)
    padding = expand(Val(2), padding)

    dims = Dims{2}(kernel, stride, dilation, padding)
    weights = Memory(_weights)
    bias = Memory(_bias)
    opdesc = Ref{Lib.dnnl_convolution_desc_t}()

    # Create attributes for fusing the postop.
    attributes = Attributes()
    postops = PostOps()
    eltwise!(postops, activation)
    append!(attributes, postops)
    return Conv(weights, bias, activation, dims, attributes, false, opdesc)
end

_slice(x::NTuple{N}) where {N} = ntuple(i -> x[(2 * i) - 1], Val(div(N, 2)))
_subone(x::NTuple{N}) where {N} = ntuple(i -> x[i] - 1, Val(N))

# Note: OneDNN's convolution is Flux's CrossCor
function Conv(m::Flux.CrossCor)
    # Flux stores its padding as a 2N tuple.
    # Since we're just using symmetrical padding for now, we grab slice every other
    # padding value.
    return Conv(
        m.weight,
        m.bias,
        m.σ;
        m.stride,
        padding = _slice(m.pad),
        dilation = _subone(m.dilation),
    )
end

Flux.@functor Conv (weights, bias)

function (conv::Conv)(_src, fuse_activation = true; kw...)
    src = Memory(_src)
    attributes = fuse_activation ? conv.attributes : noattributes()

    if conv.optimized_weights == false
        # Prepare a callback to capture the optimized memory format.
        weights_md = Ref{MemoryDesc}()
        function weight_opt_callback(_, pd)
            weights_md[] = query_md(pd, @query(weights))
            return nothing
        end

        dst = convolution(
            src,
            conv.weights,
            conv.bias,
            conv.dims;
            attributes = attributes,
            src_md = toany(src),
            weights_md = toany(conv.weights),
            callback = weight_opt_callback,
            conv.opdesc,
            kw...,
        )

        conv.weights = maybe_reorder(weights_md[], conv.weights)
        conv.optimized_weights = true
        return dst
    else
        return convolution(
            src,
            conv.weights,
            conv.bias,
            conv.dims;
            attributes,
            src_md = toany(src),
            conv.opdesc,
            kw...,
        )
    end
end

function ChainRulesCore.rrule(conv::Conv, _src::AbstractArray, _fuse_activation = false)
    # The result of `canfuse` is known at compile time, so Julia can optimize out the branch.
    src = Memory(_src)
    return canfuse(conv.activation) ? rrule_fused(conv, src) : rrule_unfused(conv, src)
end

function rrule_fused(conv::T, src) where {T<:Conv}
    src_size = size(src)
    nt = conv(src, true; kind = Training())
    @unpack dst, forward = nt
    function conv_fused_pullback(_diff_dst)
        # Maybe convert argument
        diff_dst = Memory(_diff_dst)
        # Reverse activation function.
        diff_dst_pre = eltwise_backward(conv.activation, diff_dst, dst)

        # Backprop convolution kernel
        diff_src, diff_dst_reordered = convolution_backward_data(
            src_size, conv.weights, diff_dst_pre, conv.dims; forward, conv.opdesc
        )
        (diff_weights, diff_bias) = convolution_backward_weights(
            size(conv.weights), src, diff_dst_reordered, conv.dims; forward, conv.opdesc
        )

        return (
            ChainRulesCore.Tangent{T}(; weights = diff_weights, bias = diff_bias),
            diff_src,
            ChainRulesCore.NoTangent(),
        )
    end
    return dst, conv_fused_pullback
end

function rrule_unfused(conv::T, _src) where {T<:Conv}
    src = Memory(_src)
    nt = conv(src, false; kind = Training())
    dst_pre, forward = nt.dst, nt.forward
    dst = eltwise(conv.activation, dst_pre)
    src_size = size(src)
    function conv_pullback(_diff_dst)
        diff_dst = Memory(_diff_dst)
        diff_dst_pre = eltwise_backward(conv.activation, diff_dst, dst_pre)
        (diff_src, diff_dst_reordered) = convolution_backward_data(
            src_size, conv.weights, diff_dst_pre, conv.dims; forward, conv.opdesc
        )

        (diff_weights, diff_bias) = convolution_backward_weights(
            size(conv.weights), src, diff_dst_reordered, conv.dims; forward, conv.opdesc
        )

        return (
            ChainRulesCore.Tangent{T}(; weights = diff_weights, bias = diff_bias),
            diff_src,
            ChainRulesCore.NoTangent(),
        )
    end
    return dst, conv_pullback
end


text:
using Test

@timed_testset "layers/" begin
    include("layers/conv2d.jl")
    include("layers/flatten.jl")
    include("layers/linear.jl")
    include("layers/masked_relu.jl")
    include("layers/pool.jl")
    include("layers/relu.jl")
end


text:
#= #an oop mentality was to tackle the object!
function Account()

end
=#

module TAccount
#Enum accountOp = Dr,Cr
#Enum accountType = "Debit","Credit"
@enum Dr 1 2 # 1 for debit , 2 for Credit
@enum cr 1 2
#@enum accountOp(@enum Dr 1 2, @enum Cr 1 2)
#@enum accountOp(@enum Dr 1, @enum Cr 2)

"""
  1: dr drAccount cr crAccount
  increase drAccount , increase crAccount
  infer: there must be an incoming cash in-flow

  2: dr crAccount cr drAccount
  decrease drAccount , decrease crAccount
  infer: cash outflow

"""

"""
var: TypeVar
  name: Symbol T
  lb: Union{}
  ub: Integer <: Real
body: Enum{T<:Integer} <: Any
"""

dump(Enum)

function increment(account, amount)
  account = account + amount

  return account
end

function decrement(account, amount)
  account = account - amount
end
# maybe those for the  t-account
function DebitAccountAdd(date, drAccount, amount)
  drAccount.
  end

  function CreditAccount(date, account, amount)

  end
  """

  """

  function transaction(Dr, Cr, DebitAccount, CreditAccount, amount)

    """
    """

  end
  """

#transfers between accounts of the same type(Debit, or credit)

dr drAccount1 cr drAccount2

dr crAccount1 cr crAccount2

"""
function transactionRead(Dr, Cr, DebitAccount, CreditAccount, amount)
  # 1: dr drAccount cr crAccount
if Dr == 1 & cr == 1 #increment account
  DebitAccount.add(amount)
  CreditAccount.add(amount)
elseif Dr == 2 & cr == 1

end

function transactionRead(Dr, Cr, CreditAccount, CreditAccount2, amount)

end
function transactionRead(Dr, Cr, DebitAccount, DebitAccount2, amount)

end

  function transaction(DebitAccount, DebitAccount2, amount)
    """
    deduct from first account, add to the second account
    """
    DebitAccount.add

  end
  """
  transaction Module

  each Transation Ideally must get written in 3 locations (at least ):
  forming the
  Accounting Triad
  1. A Journal (daily journal of transactions) -of the self (entity)
  2. A Debit Dr. Book
  3. A Credit Cr. Book

  Books can increase + or decrease in value -

  but each increment in a book has a

  """

  function transaction(accountType, date, DrAccount, CrAccount)

    global Id::Int # Auto-Increments

    function run()
      DebitAmount(date, amount)
      CreditAmount(date, amount)

    end

    function DebitAccount(date, amount)
      #TODO: add data row of CreditAccount here
      #date Float64 amount
    end

    function CreditAccount(date, amount)
      #TODO: add data row of DebitAccount here
      # Date Float64 amount
    end

    function BalanceAccount()
      #TODO: check if today is last day in month
      bal_cd # @end of month Only - balance carried douwnd

    end
    function bal_cd()

    end

  end

  #transaction()
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

#############################################################################

# lattice visualizations are drawn with Luxor because it is easier to do 2D drawings with Luxor than with Makie.

function draw(tilebasis::AbstractBasis,tilesize,i,j,color,name)
    vertices = Repeat.tilevertices(tilebasis)
    tile = tilesize*[Luxor.Point(vertices[:,i]...) for i in 1:size(vertices)[2]]
    pt = tilesize*tilebasis[i,j]
    offset = Luxor.Point(pt[1],-pt[2]) #flip y so indices show up correctly

    Luxor.translate(offset)

    Luxor.sethue(color)

    Luxor.poly(tile, :fill, close=true)
    Luxor.sethue("lightgrey")
    Luxor.poly(tile, :stroke, close=true)
    Luxor.sethue("black")

    #scale and offset text so coordinates are readable
    Luxor.fontsize(tilesize/3)
    Luxor.text("$i, $j",Luxor.Point(-tilesize/3,tilesize/2.5))
    if name !== nothing
        Luxor.text(name,Luxor.Point(-tilesize/4,-tilesize/4))
    end

    Luxor.translate(-offset)
end

"""Draws a list of hexagonal cells, represented by their lattice coordinates"""
function drawcells(tilebasis::AbstractBasis, tilesize,cells; color::Union{AbstractArray,Nothing} = nothing, name::Union{AbstractArray{String},Nothing} = nothing, format=:png, resolution=(500,500))

    numcells = size(cells)[2]
    Luxor.Drawing(resolution[1], resolution[2], format)
    Luxor.origin()
    Luxor.background(Colors.RGBA(0, 1, 1, 0.0))
    if color === nothing
        color = Colors.distinguishable_colors(numcells,lchoices = 30:250) #type unstable but not performance critical code
    end

    for i in 1:numcells
        cell = cells[:,i]
        cellname = name === nothing ? nothing : name[i]
        draw(tilebasis,tilesize,cell[1],cell[2],color[i],cellname)
    end

    Luxor.finish()
    if (format == :svg)
        res = Luxor.svgstring()
        return res
    end
    if (format == :png)
        Luxor.preview()
    end
end

""" draw the LatticeCluster offset to (0,0) """
draw(clstr::LatticeCluster,tilesize = 50.0, cells = hcat([0,0])) = drawcells(clstr.elementbasis,tilesize,clustercoordinates(clstr,0,0))

""" draw the ClusterWithProperties at coordinates specified by cells """
function draw(clstr::Repeat.ClusterWithProperties, cells = hcat([0,0]), scale = 50.0) #this is how you have to make a 2x1 matrix. One would expect [0;0] to work but it doesn't.
    dims = size(cells)
    @assert dims[1] == 2 "dims[1] must be 2 instead was $(dims[1])"
    clstrsize = clustersize(clstr)
    points = Matrix(undef,dims[1],dims[2]*clstrsize)
    for i in 1:dims[2]
        points[:,(i-1)*clstrsize+1:i*clstrsize] = clustercoordinates(clstr,cells[1,i],cells[2,i])
    end
    props = repeat(Repeat.properties(clstr),dims[2])
    drawcells(elementbasis(cluster(clstr)),scale,points,color = props[:,:Color],name = props[:,:Name])
end


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")
include("CESMReader.jl")

using .CESMReader

function mreplace(x)
    return replace(x, missing=>NaN)
end

function integrate(x, dydx)
    y = copy(dydx) * 0.0

    for i = 2:length(x)
        y[i] = y[i-1] + (dydx[i-1] + dydx[i]) * (x[i] - x[i-1]) / 2.0
    end

    return y
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds

    global mask, xc, yc, area, ni, nj
    mask = ds["mask"][:] |> nomissing
    xc   = ds["xc"][:]   |> nomissing
    yc   = ds["yc"][:]   |> nomissing
    area = ds["area"][:] |> nomissing
    ni, nj = size(area)

end

idx       = mask .== 1.0
_area     = area[idx]
_area_sum = sum(_area)

calAreaAvg = ( field , ) -> sum(_area .* (field[idx])) / _area_sum
calAreaWeightedSum = ( field , ) -> sum(_area .* (field[idx]))
function calMonthlyAreaAvg(field)
    return convert(Array{Float64}, [calAreaAvg(view(field, :, :, m)) for m=1:12])
end

let

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    _TEMP, _SALT = getData(fh, ["TEMP", "SALT"], (parsed["beg-year"], parsed["end-year"]), (:, :))

    global Nt = size(_TEMP)[3]

    if mod(Nt, 12) != 0
        throw(ErrorException("Time should be a multiple of 12."))
    end

    global nyears = Int64(Nt / 12)
    println(format("We got {:d} years of data.", nyears))

    global TEMP = zeros(Float64, Nt)
    global SALT = zeros(Float64, Nt)

    for t=1:Nt
        TEMP[t] = view(_TEMP, :, :, t) |> calAreaWeightedSum
        SALT[t] = view(_SALT, :, :, t) |> calAreaWeightedSum
    end

end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in [
        ("TEMP",  TEMP,  ("time",), Dict()),
        ("SALT",  SALT,  ("time",), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
@setup begin
  this_should_not_leak = true
end

@test "this file is loaded" begin
  @assert 1 == 1
end

@teardown begin
  this_also_should_not_leak = true
end


text:
"""
	$(SIGNATURES)

Plots a set of regressions as a series of bar graph. Each bar is a regressor so that different scales are not a visual problem.
Standard errors are shown as error bars.

# Arguments

- coeffNameV:  Regressor names. Common across regressions. Symbols or Strings.
- regrLabelV:  Descriptive labels for the regressions. Become x labels under each bar. Symbols or Strings.
- coeffM: Coefficients by [regressor, regression].
- seM: Standard errors in the same format.
- `onePlot`: if `true`, all regressors are plotted in one grouped bar graph. Only works when scales are similar across regressors.
"""
function plot_regressions(coeffNameV, regrLabelV,
    coeffM :: AbstractMatrix{F1},
    seM :: AbstractMatrix{F1};
    onePlot :: Bool = false) where F1 <: AbstractFloat

    n, nRegr = size(coeffM);
    @assert length(coeffNameV) == n
    @assert length(regrLabelV) == nRegr

    if onePlot
        p = plot_regressions_together(coeffNameV, regrLabelV,
            coeffM, seM);
    else
        pV = Vector{Any}(undef, n);
        for ir = 1 : n
            pV[ir] = plot_regression(string.(regrLabelV), coeffNameV[ir],
                coeffM[ir,:], seM[ir,:]);
        end
        p = plot(pV...);
    end
    return p
end

function plot_regression(xStrV, yStr, coeffV :: AbstractVector{F1},
    seV :: AbstractVector{F1}) where F1 <: AbstractFloat

    @assert length(xStrV) == length(coeffV) == length(seV)  "Size mismatch"

    p = bar(xStrV, coeffV, yerr = seV,
        markerstrokewidth = 2, markerstrokecolor = :black,
        linecolor = :black,  legend = false,  ylabel = yStr,
        lw = 0);
    return p
end

# Plot regressions in one plot. Each group of bars is a regressor.
function plot_regressions_together(coeffNameV, regrLabelV, coeffM, seM)
    p = grouped_bar_xy(coeffM',  string.(coeffNameV);
        legendV = regrLabelV, legendPos = :best,
        yErrorM = seM');
    return p
end

# -----------------


text:
using FastTransforms
using Compat.Test
import FastTransforms.pochhammer

const sqrtpi = 1.772453850905516027298

@testset "Pochhammer" begin
    @test pochhammer(2,3) == 24
    @test pochhammer(0.5,3) == 0.5*1.5*2.5
    @test pochhammer(0.5,0.5) == 1/sqrtpi
    @test pochhammer(0,1) == 0
    @test pochhammer(-1,2) == 0
    @test pochhammer(-5,3) == -60
    @test pochhammer(-1,-0.5) == 0
    @test 1.0/pochhammer(-0.5,-0.5) == 0
    @test pochhammer(-1+0im,-1) == -0.5
end


text:
# Global method definition needs to be at top level in .7

"""
Takes a PowerSystems dictionary and return an array of PowerSystems struct for
Bus, Generator, Branch and load
"""
function ps_dict2ps_struct(data::Dict{String,Any})
    generators = Array{G where {G<:Generator},1}()
    storages = Array{S where {S<:Storage},1}()
    buses = Array{Bus,1}()
    branches = Array{B where {B<:Branch},1}()
    loads = Array{E where {E<:ElectricLoad},1}()
    shunts = Array{FixedAdmittance,1}()
    loadZones = Array{D where {D<:Device},1}()
    services = Array{S where {S<:Service},1}()

    # TODO: should we raise an exception in the following?

    if haskey(data, "bus")
        buses = PowerSystems.bus_dict_parse(data["bus"])
    else
        @warn "key 'bus' not found in PowerSystems dictionary, this will result in an empty Bus array"
    end
    if haskey(data, "gen")
        (generators, storage) = PowerSystems.gen_dict_parser(data["gen"])
    else
        @warn "key 'gen' not found in PowerSystems dictionary, this will result in an empty Generators and Storage array"
    end
    if haskey(data, "branch")
        branches = PowerSystems.branch_dict_parser(data["branch"], branches)
    else
        @warn "key 'branch' not found in PowerSystems dictionary, this will result in an empty Branches array"
    end
    if haskey(data, "load")
        loads = PowerSystems.load_dict_parser(data["load"])
    else
        @warn "key 'load' not found in PowerSystems dictionary, this will result in an empty Loads array"
    end
    if haskey(data, "loadzone")
        loadZones = PowerSystems.loadzone_dict_parser(data["loadzone"])
    else
        @warn "key 'loadzone' not found in PowerSystems dictionary, this will result in an empty LoadZones array"
    end
    if haskey(data, "shunt")
        shunts = PowerSystems.shunt_dict_parser(data["shunt"])
    else
        @warn "key 'shunt' not found in PowerSystems dictionary, this will result in an empty Shunts array"
    end
    if haskey(data, "dcline")
        branches = PowerSystems.dclines_dict_parser(data["dcline"], branches)
    else
        @warn "key 'dcline' not found in PowerSystems dictionary, this will result in an empty DCLines array"
    end
    if haskey(data, "services")
        services = PowerSystems.services_dict_parser(data["services"],generators)
    else
        @warn "key 'services' not found in PowerSystems dictionary, this will result in an empty services array"
    end

    return sort!(buses, by = x -> x.number), generators, storage,  sort!(branches, by = x -> x.connectionpoints.from.number), loads, loadZones, shunts, services

end

function _retrieve(dict::T, key_of_interest::String, output = Dict(), path = []) where T<:AbstractDict
    iter_result = Base.iterate(dict)
    last_element = length(path)
    while iter_result !== nothing
        ((key,value), state) = iter_result
        if key == key_of_interest
            output[value] = !haskey(output,value) ? path[1:end] : push!(output[value],path[1:end])
        end
        if value isa AbstractDict
            push!(path,key)
            _retrieve(value, key_of_interest, output, path)
            path = path[1:last_element]
        end
        iter_result = Base.iterate(dict, state)
    end
    return output
end

function _retrieve(dict::T, type_of_interest, output = Dict(), path = []) where T<:AbstractDict
    iter_result = Base.iterate(dict)
    last_element = length(path)
    while iter_result !== nothing
        ((key,value), state) = iter_result
        if typeof(value) <: type_of_interest
            output[key] = !haskey(output,value) ? path[1:end] : push!(output[value],path[1:end])
        end
        if value isa AbstractDict
            push!(path,key)
            _retrieve(value, type_of_interest, output, path)
            path = path[1:last_element]
        end
        iter_result = Base.iterate(dict, state)
    end
    return output
end

function _access(nesteddict::T,keylist) where T<:AbstractDict
    if !haskey(nesteddict,keylist[1])
        @error "$(keylist[1]) not found in dict"
    end
    if length(keylist) > 1
        nesteddict = _access(nesteddict[keylist[1]],keylist[2:end])
    else
        nesteddict = nesteddict[keylist[1]]
    end
end

function _get_device(name, collection, devices = [])
    if isa(collection,Array)
        fn = fieldnames(typeof(collection[1]))
        if :name in fn
            [push!(devices,d) for d in collection if d.name == name]
        end
    else
        fn = fieldnames(typeof(collection))
        for f in fn
            _get_device(name,getfield(collection,f),devices)
        end
    end
    return devices
end

"""
Arg:
    Dataframes which includes a timerseries columns of either:
        Year, Month, Day, Period
      or
        DateTime
      or
        nothing (creates a today referenced DateTime Column)
Returns:
    Dataframe with a DateTime columns
"""
function read_datetime(df; kwargs...)
    if [c for c in [:Year,:Month,:Day,:Period] if c in names(df)] == [:Year,:Month,:Day,:Period]
        if Dates.Hour(DataFrames.maximum(df[:Period])) <= Dates.Hour(25)
            df[:DateTime] = collect(Dates.DateTime(df[1,:Year],df[1,:Month],df[1,:Day],(df[1,:Period]-1)) :Dates.Hour(1) :
            Dates.DateTime(df[end,:Year],df[end,:Month],df[end,:Day],(df[end,:Period]-1)))
        elseif (Dates.Minute(5) * DataFrames.maximum(df[:Period]) >= Dates.Minute(1440))& (Dates.Minute(5) * DataFrames.maximum(df[:Period]) <= Dates.Minute(1500))
            df[:DateTime] = collect(Dates.DateTime(df[1,:Year],df[1,:Month],df[1,:Day],floor(df[1,:Period]/12),Int(df[1,:Period])-1) :Dates.Minute(5) :
                            Dates.DateTime(df[end,:Year],df[end,:Month],df[end,:Day],floor(df[end,:Period]/12)-1,5*(Int(df[end,:Period])-(floor(df[end,:Period]/12)-1)*12) -5))
        else
            @error "I don't know what the period length is, reformat timeseries"
        end
        DataFrames.deletecols!(df, [:Year,:Month,:Day,:Period])

    elseif :DateTime in names(df)
        df[:DateTime] = Dates.DateTime(df[:DateTime])
    else
        if :startdatetime in keys(kwargs)
            startdatetime = kwargs[:startdatetime]
        else
            @warn "No reference date given, assuming today"
            startdatetime = Dates.today()
        end
        df[:DateTime] = collect(Dates.DateTime(startdatetime):Dates.Hour(1):Dates.DateTime(startdatetime)+Dates.Hour(size(df)[1]-1))
    end
    return df
end

"""
Arg:
    Device dictionary - Generators
    Dataframe contains device Realtime/Forecast TimeSeries
Returns:
    Device dictionary with timeseries added
"""
function add_time_series(Device_dict::Dict{String,Any}, df::DataFrames.DataFrame)
    for (device_key,device) in Device_dict
        if device_key in map(string, names(df))
            ts_raw = df[Symbol(device_key)]
            if maximum(ts_raw) <= 1.0
                @info "assumed time series is a scaling factor for $device_key"
                Device_dict[device_key]["scalingfactor"] = TimeSeries.TimeArray(df[:DateTime],ts_raw)
            else
                @info "assumed time series is MW for $device_key"
                Device_dict[device_key]["scalingfactor"] = TimeSeries.TimeArray(df[:DateTime],ts_raw/device["tech"]["installedcapacity"])
            end
        end
    end
    return Device_dict
end

function add_time_series(Device_dict::Dict{String,Any}, ts_raw::TimeSeries.TimeArray)
    """
    Arg:
        Device dictionary - Generators
        Dict contains device Realtime/Forecast TimeSeries.TimeArray
    Returns:
        Device dictionary with timeseries added
    """

    name = get(Device_dict, "name", "")
    if name == ""
        throw(DataFormatError("input dict to add_time_series in wrong format"))
    end

    if maximum(values(ts_raw)) > 1.0
        @warn "Time series for $name has values > 1.0, expected values in range {0.0,1.0}"
    end
    Device_dict["scalingfactor"] = ts_raw

    return Device_dict
end

"""
Arg:
    Load dictionary
    LoadZones dictionary
    Dataframe contains device Realtime/Forecast TimeSeries
Returns:
    Device dictionary with timeseries added
"""
function add_time_series_load(data::Dict{String,Any}, df::DataFrames.DataFrame)
    load_dict = data["load"]

    load_names = [string(l["name"]) for (k,l) in load_dict]
    ts_names = [string(n) for n in names(df) if n != :DateTime]

    write_sf_by_lz = false
    lzkey = [k for k in ["loadzone","load_zone"] if haskey(data,k)][1]
    if lzkey in keys(data)
        load_zone_dict = data[lzkey]
        z_names = [string(z["name"]) for (k,z) in load_zone_dict]
        if length([n for n in z_names if n in ts_names]) > 0
            write_sf_by_lz = true
        end
    end

    assigned_loads = []
    if write_sf_by_lz
        @info "assigning load scaling factors by load_zone"
        # TODO: make this faster/better
        for (l_key,l) in load_dict
            for (lz_key,lz) in load_zone_dict
                if l["bus"] in lz["buses"]
                    ts_raw = df[lz_key]/lz["maxactivepower"]
                    load_dict[l_key]["scalingfactor"] = TimeSeries.TimeArray(df[:DateTime],ts_raw)
                    push!(assigned_loads,l_key)
                end
            end

        end
    else
        @info "assigning load scaling factors by bus"
        for (l_key,l) in load_dict
            load_dict[l_key]["scalingfactor"] = TimeSeries.TimeArray(df[:DateTime],df[Symbol(l["name"])])
            push!(assigned_loads,l["name"])
        end

    end

    for l in [l for l in load_names if !(l in assigned_loads)]
        @warn "No load scaling factor assigned for $l" maxlog=PS_MAX_LOG
    end

    return load_dict
end

## - Parse Dict to Struct
function bus_dict_parse(dict::Dict{Int,Any})
    Buses = [Bus(b["number"],b["name"], b["bustype"],b["angle"],b["voltage"],b["voltagelimits"],b["basevoltage"]) for (k_b,b) in dict ]
    return Buses
end

## - Parse Dict to Array
function gen_dict_parser(dict::Dict{String,Any})
    Generators = Array{G where {G<:Generator},1}()
    Storage_gen = Array{S where {S<:Storage},1}()
    for (gen_type_key,gen_type_dict) in dict
        if gen_type_key =="Thermal"
            for (thermal_key,thermal_dict) in gen_type_dict
                push!(Generators,ThermalDispatch(string(thermal_dict["name"]),
                                                            Bool(thermal_dict["available"]),
                                                            thermal_dict["bus"],
                                                            TechThermal(thermal_dict["tech"]["activepower"],
                                                                        thermal_dict["tech"]["activepowerlimits"],
                                                                        thermal_dict["tech"]["reactivepower"],
                                                                        thermal_dict["tech"]["reactivepowerlimits"],
                                                                        thermal_dict["tech"]["ramplimits"],
                                                                        thermal_dict["tech"]["timelimits"]),
                                                            EconThermal(thermal_dict["econ"]["capacity"],
                                                                        thermal_dict["econ"]["variablecost"],
                                                                        thermal_dict["econ"]["fixedcost"],
                                                                        thermal_dict["econ"]["startupcost"],
                                                                        thermal_dict["econ"]["shutdncost"],
                                                                        thermal_dict["econ"]["annualcapacityfactor"])
                            ))
            end
        elseif gen_type_key =="Hydro"
            for (hydro_key,hydro_dict) in gen_type_dict
                push!(Generators,HydroCurtailment(string(hydro_dict["name"]),
                                                            hydro_dict["available"],
                                                            hydro_dict["bus"],
                                                            TechHydro(  hydro_dict["tech"]["installedcapacity"],
                                                                        hydro_dict["tech"]["activepower"],
                                                                        hydro_dict["tech"]["activepowerlimits"],
                                                                        hydro_dict["tech"]["reactivepower"],
                                                                        hydro_dict["tech"]["reactivepowerlimits"],
                                                                        hydro_dict["tech"]["ramplimits"],
                                                                        hydro_dict["tech"]["timelimits"]),
                                                            hydro_dict["econ"]["curtailcost"],
                                                            hydro_dict["scalingfactor"]
                            ))
            end
        elseif gen_type_key =="Renewable"
            for (ren_key,ren_dict) in  gen_type_dict
                if ren_key == "PV"
                    for (pv_key,pv_dict) in ren_dict
                        push!(Generators,RenewableCurtailment(string(pv_dict["name"]),
                                                                    Bool( pv_dict["available"]),
                                                                    pv_dict["bus"],
                                                                    pv_dict["tech"]["installedcapacity"],
                                                                    EconRenewable(pv_dict["econ"]["curtailcost"],
                                                                                pv_dict["econ"]["interruptioncost"]),
                                                                    pv_dict["scalingfactor"]
                                    ))
                    end
                elseif ren_key == "RTPV"
                    for (rtpv_key,rtpv_dict) in ren_dict
                        push!(Generators,RenewableFix(string(rtpv_dict["name"]),
                                                                    Bool(rtpv_dict["available"]),
                                                                    rtpv_dict["bus"],
                                                                    rtpv_dict["tech"]["installedcapacity"],
                                                                    rtpv_dict["scalingfactor"]
                                    ))
                    end
                elseif ren_key == "WIND"
                    for (wind_key,wind_dict) in ren_dict
                        push!(Generators,RenewableCurtailment(string(wind_dict["name"]),
                                                                    Bool(wind_dict["available"]),
                                                                    wind_dict["bus"],
                                                                    wind_dict["tech"]["installedcapacity"],
                                                                    EconRenewable(wind_dict["econ"]["curtailcost"],
                                                                                wind_dict["econ"]["interruptioncost"]),
                                                                    wind_dict["scalingfactor"]
                                    ))
                    end
                end
            end
        elseif gen_type_key =="Storage"
            for (storage_key,storage_dict) in  gen_type_dict
                push!(Storage_gen,GenericBattery(string(storage_dict["name"]),
                                                            Bool(storage_dict["available"]),
                                                            storage_dict["bus"],
                                                            storage_dict["energy"],
                                                            storage_dict["capacity"],
                                                            storage_dict["activepower"],
                                                            storage_dict["inputactivepowerlimits"],
                                                            storage_dict["outputactivepowerlimits"],
                                                            storage_dict["efficiency"],
                                                            storage_dict["reactivepower"],
                                                            storage_dict["reactivepowerlimits"]
                            ))
            end
        end
    end
    return (Generators, Storage_gen)
end

# - Parse Dict to Array

function branch_dict_parser(dict::Dict{String,Any},Branches::Array{B,1}) where {B<:Branch}
    for (branch_key,branch_dict) in dict
        if branch_key == "Transformers"
            for (trans_key,trans_dict) in branch_dict
                if trans_dict["tap"] ==1.0
                    push!(Branches,Transformer2W(string(trans_dict["name"]),
                                                Bool(trans_dict["available"]),
                                                trans_dict["connectionpoints"],
                                                trans_dict["r"],
                                                trans_dict["x"],
                                                trans_dict["primaryshunt"],
                                                trans_dict["rate"]
                                                ))
                elseif trans_dict["tap"] !=1.0
                    alpha = "α" in keys(trans_dict) ? trans_dict["α"] : 0.0
                    if alpha !=0.0 #TODO : 3W Transformer
                        push!(Branches,PhaseShiftingTransformer(string(trans_dict["name"]),
                                                    Bool(trans_dict["available"]),
                                                    trans_dict["connectionpoints"],
                                                    trans_dict["r"],
                                                    trans_dict["x"],
                                                    trans_dict["primaryshunt"],
                                                    trans_dict["tap"],
                                                    trans_dict["α"],
                                                    trans_dict["rate"]
                                                    ))
                    else
                        push!(Branches,TapTransformer(string(trans_dict["name"]),
                                                    Bool(trans_dict["available"]),
                                                    trans_dict["connectionpoints"],
                                                    trans_dict["r"],
                                                    trans_dict["x"],
                                                    trans_dict["primaryshunt"],
                                                    trans_dict["tap"],
                                                    trans_dict["rate"]
                                                    ))
                    end
                end
            end
        else branch_key == "Lines"
            for (line_key,line_dict) in branch_dict
                push!(Branches,Line(string(line_dict["name"]),
                                    Bool(line_dict["available"]),
                                    line_dict["connectionpoints"],
                                    line_dict["r"],
                                    line_dict["x"],
                                    line_dict["b"],
                                    float(line_dict["rate"]),
                                    line_dict["anglelimits"]
                                    ))
            end
        end
    end
    return Branches
end

function load_dict_parser(dict::Dict{String,Any})
    Loads =Array{L where {L<:ElectricLoad},1}()
    for (load_key,load_dict) in dict
        push!(Loads,PowerLoad(string(load_dict["name"]),
                Bool(load_dict["available"]),
                load_dict["bus"],
                load_dict["maxactivepower"],
                load_dict["maxreactivepower"],
                load_dict["scalingfactor"]
                ))
    end
    return Loads
end

function loadzone_dict_parser(dict::Dict{Int64,Any})
    LoadZs =Array{D where {D<:Device},1}()
    for (lz_key,lz_dict) in dict
        push!(LoadZs,LoadZones(lz_dict["number"],
                                string(lz_dict["name"]),
                                lz_dict["buses"],
                                lz_dict["maxactivepower"],
                                lz_dict["maxreactivepower"]
                                ))
    end
    return LoadZs
end

function shunt_dict_parser(dict::Dict{String,Any})
    Shunts = Array{FixedAdmittance,1}()
    for (s_key,s_dict) in dict
        push!(Shunts,FixedAdmittance(string(s_dict["name"]),
                            Bool(s_dict["available"]),
                            s_dict["bus"],
                            s_dict["Y"]
                            )
            )
    end
    return Shunts
end

function dclines_dict_parser(dict::Dict{String,Any},Branches::Array{Branch,1})
    for (dct_key,dct_dict) in dict
        if dct_key == "HVDCLine"
            for (dcl_key,dcl_dict) in dct_dict
                push!(Branches,HVDCLine(string(dcl_dict["name"]),
                                    Bool(dcl_dict["available"]),
                                    dcl_dict["connectionpoints"],
                                    dcl_dict["activepowerlimits_from"],
                                    dcl_dict["activepowerlimits_to"],
                                    dcl_dict["reactivepowerlimits_from"],
                                    dcl_dict["reactivepowerlimits_to"],
                                    dcl_dict["loss"]
                                    ))
            end
        elseif dct_key == "VSCDCLine"
            for (dcl_key,dcl_dict) in dct_dict
                push!(Branches,VSCDCLine(string(dcl_dict["name"]),
                                    Bool(dcl_dict["available"]),
                                    dcl_dict["connectionpoints"],
                                    dcl_dict["rectifier_taplimits"],
                                    dcl_dict["rectifier_xrc"],
                                    dcl_dict["rectifier_firingangle"],
                                    dcl_dict["inverter_taplimits"],
                                    dcl_dict["inverter_xrc"],
                                    dcl_dict["inverter_firingangle"]
                                    ))
            end
        end
    end
    return Branches
end

function services_dict_parser(dict::Dict{String,Any},generators::Array{Generator,1})
    Services = Array{D where {D <: Service},1}()

    for (k,d) in dict
        contributingdevices = Array{D where {D<:PowerSystems.Device},1}()
        [PowerSystems._get_device(dev,generators) for dev in d["contributingdevices"]] |> (x->[push!(contributingdevices,d[1]) for d in x if length(d)==1])
        push!(Services,ProportionalReserve(d["name"],
                            contributingdevices,
                            Float64(d["timeframe"]),
                            TimeSeries.TimeArray(Dates.today(),ones(1))  #TODO : fix requirement
                            ))
    end
    return Services
end


text:
"""
    power_method(x0,Hop,PARAM)
POWER_METHOD: Power iteration method to computes max eigenvalue of H'H
              where H and H' are given by Hop. This function is
              needed to evaluate the step parameter of FISTA.
# Arguments
- `x0`:initial seed with dimensions such that H'*H*x0 does not abort
- `Hop`:liner operator that encapsulates H and H' such that
        Hop(x,PARAM, 1) = H x
        Hop(y,PARAM,-1) = H'y
- `PARAM`:parameters to run Hop and it's adjoint

"""
function power_method(x0,Hop,PARAM)
   #x = x0;
   x=copy(x0);
   value=0;
    for k = 1:10;
       aux = Hop(x,PARAM,1);
       y = Hop(aux,PARAM,-1);
       n = norm(y);
       x = y/n;
       value = n;
    end
    return value
end


text:
#################################
## Generar cromosoma
##########################

"""
    GeneraCromosoma( tamaño, umbral )
Cromosoma con componentes uniformes
"""
function GeneraCromosoma( tamaño, umbral )
    cromosoma = rand(tamaño)
    return map( x -> x < umbral ? 0.0 : x, cromosoma)
end


text:
module ACPlots
export
    wiggle,
    wiggle!
using ..ACWT, ..ACTransforms
using LinearAlgebra, Plots, Wavelets, Statistics

## Visualizations
"""
    wiggle(wav; taxis=1:size(wav,1), zaxis=1:size(wav,2), sc=1, EdgeColor=:black, FaceColor=:black, Orient=:across, Overlap=true, ZDir=:normal)

Plot a set of shaded wiggles

# Arguments
- `wav::Array`: matrix of waveform columns.
- `taxis::Array=1:size(wav,1)`: time axis vector
- `zaxis::Array=1:size(wav,2)`: space axis vector
- `sc::Float=1`: scale factor/magnification.
- `EdgeColor::Symbol=:black`: Sets edge of wiggles color.
- `FaceColor::Symbol=:black`: Sets shading color of wiggles.
- `Overlap::bool=true`: How signals are scaled.
        true  - Signals overlap (default);
        false - Signals are scaled so they do not overlap.
- `Orient::Symbol=:across`: Controls orientation of wiggles.
        :across - from left to right
        :down   - from top to down
- `ZDir::Symbol=:normal`: Direction of space axis.
        :normal  - First signal at bottom (default)
        :reverse - First signal at top.

Translated by Nicholas Hausch -- MATLAB file provided by Naoki Saito
The previous MATLAB version contributors are:
 Anthony K. Booer (SLB) and Bradley Marchand (NSWC-PC)
Revised by Naoki Saito, Feb. 05, 2018
"""
function wiggle(wav::AbstractArray{Float64};
                taxis=1:size(wav,1), zaxis=1:size(wav,2),
                sc::Integer=1, EdgeColor::Symbol=:black, FaceColor::Symbol=:black,
                Overlap::Bool=true, Orient::Symbol=:across, ZDir::Symbol=:normal)

    # Set axes
    (n,m) = size(wav)

    # Sanity check
    if length(taxis) != n
        error("Inconsistent taxis dimension!")
    end
    if length(zaxis) != m
        error("Inconsistent zaxis dimension!")
    end

    # For calculation purposes
    maxrow = zeros(m); minrow = zeros(m)
    for k = 1:m
      maxrow[k] = maximum(wav[:,k]); minrow[k] = minimum(wav[:,k])
    end

    # Scale the data for plotting
    wamp = deepcopy(wav)
    dt = mean(diff(taxis))
    dz = mean(diff(zaxis))
    if Overlap
      wamp *= 2 * dz * (sc/maximum(maxrow-minrow))
    else
      wmax = maximum(maxrow); wmin = minimum(minrow);
      if wmax<=0
        wmax = 0
      end
      if wmin>=0
        wmin = 0
      end
        wamp = sc*wav/(wmax-wmin)
    end

    # Set initial plot
    t0 = minimum(taxis)
    t1 = maximum(taxis)
    z0 = minimum(zaxis)
    z1 = maximum(zaxis)
    if Orient == :down
     plot(xlims=(z0-dz,z1+dz), ylims=(t0,t1), yflip=true, legend=:none)
    else
     plot(xlims=(t0,t1), ylims=(z0-dz,z1+dz), legend=:none)
    end
    if ZDir == :reverse
        wamp = flipdim(wamp,2)
    end

    # Plot each wavelet
    for k = 1:m
      sig = wamp[:,k]
      t = deepcopy(taxis)
      w_sign = sign.(sig)
      for j=1:n-1
        if (w_sign[j]!=w_sign[j+1] && w_sign[j]!=0 && w_sign[j+1]!=0)
          sig = [sig; 0]
          t = [t; t[j]-sig[j]*(t[j+1]-t[j])/(sig[j+1]-sig[j])]
        end
      end
      IX = sortperm(t)
      t = t[IX]
      sig = sig[IX]
      len = length(t)
      len1 = collect(len:-1:1)
      indperm = [1:len;len1]
      inputx = t[indperm]
      inputy = zaxis[k] .+ [sig;min.(sig[len1],0)]
        # In the plot! functions below, theoretically speaking, either
        # fillrange = zaxis[k] or fillrange=[zaxis[k], zaxis[k]+dz] should be used.
        # However, those do not generate the desired plots as of O3/19/2018.
        # Somehow, the relative value of 0, i.e., fillrange=0, works well,
        # which is used temporarily.
      if Orient == :down
        plot!(inputy, inputx, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor, orientation=:v)
      else
        plot!(inputx, inputy, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor)
      end
    end
    plot!() # flushing the display.
end

"""
    wiggle!(wav; taxis=1:size(wav,1), zaxis=1:size(wav,2), sc=1, EdgeColor=:black, FaceColor=:black, Orient=:across, Overlap=true, ZDir=:normal)

Plot a set of shaded wiggles on the current displayed graphics

# Arguments
- `wav::Array`: matrix of waveform columns.
- `taxis::Array=1:size(wav,1)`: time axis vector
- `zaxis::Array=1:size(wav,2)`: space axis vector
- `sc::Float=1`: scale factor/magnification.
- `EdgeColor::Symbol=:black`: Sets edge of wiggles color.
- `FaceColor::Symbol=:black`: Sets shading color of wiggles.
- `Overlap::bool=true`: How signals are scaled.
        true  - Signals overlap (default);
        false - Signals are scaled so they do not overlap.
- `Orient::Symbol=:across`: Controls orientation of wiggles.
        :across - from left to right
        :down   - from top to down
- `ZDir::Symbol=:normal`: Direction of space axis.
        :normal  - First signal at bottom (default)
        :reverse - First signal at top.

Translated by Nicholas Hausch -- MATLAB file provided by Naoki Saito
The previous MATLAB version contributors are:
 Anthony K. Booer (SLB) and Bradley Marchand (NSWC-PC)
Revised by Naoki Saito, Feb. 05, 2018
"""
function wiggle!(wav; taxis=1:size(wav,1), zaxis=1:size(wav,2), sc=1, EdgeColor=:black, FaceColor=:black, Overlap=true, Orient=:across, ZDir=:normal)

    # Set axes
    (n,m) = size(wav)

    # Sanity check
    if length(taxis) != n
        error("Inconsistent taxis dimension!")
    end
    if length(zaxis) != m
        error("Inconsistent zaxis dimension!")
    end

    # For calculation purposes
    maxrow = zeros(m); minrow = zeros(m)
    for k = 1:m
      maxrow[k] = maximum(wav[:,k]); minrow[k] = minimum(wav[:,k])
    end

    # Scale the data for plotting
    wamp = deepcopy(wav)
    dt = mean(diff(taxis))
    dz = mean(diff(zaxis))
        if Overlap
      wamp *= 2 * dz * (sc/maximum(maxrow-minrow))
    else
      wmax = maximum(maxrow); wmin = minimum(minrow);
      if wmax<=0
        wmax = 0
      end
      if wmin>=0
        wmin = 0
      end
        wamp = sc*wav/(wmax-wmin)
    end

    # Set initial plot
    t0 = minimum(taxis)
    t1 = maximum(taxis)
    z0 = minimum(zaxis)
    z1 = maximum(zaxis)
    if Orient == :down
     plot!(xlims=(z0-dz,z1+dz), ylims=(t0,t1), yflip=true, legend=:none)
    else
     plot!(xlims=(t0,t1), ylims=(z0-dz,z1+dz), legend=:none)
    end
    if ZDir == :reverse
        wamp = flipdim(wamp,2)
    end

    # Plot each wavelet
    for k = 1:m
      sig = wamp[:,k]
      t = deepcopy(taxis)
      w_sign = sign.(sig)
      for j=1:n-1
        if (w_sign[j]!=w_sign[j+1] && w_sign[j]!=0 && w_sign[j+1]!=0)
          sig = [sig; 0]
          t = [t; t[j]-sig[j]*(t[j+1]-t[j])/(sig[j+1]-sig[j])]
        end
      end
      IX = sortperm(t)
      t = t[IX]
      sig = sig[IX]
      len = length(t)
      len1 = collect(len:-1:1)
      indperm = [1:len;len1]
      inputx = t[indperm]
      inputy = zaxis[k] .+ [sig;min.(sig[len1],0)]
        # In the plot! functions below, theoretically speaking, either
        # fillrange = zaxis[k] or fillrange=[zaxis[k], zaxis[k]+dz] should be used.
        # However, those do not generate the desired plots as of O3/19/2018.
        # Somehow, the relative value of 0, i.e., fillrange=0, works well,
        # which is used temporarily.
      if Orient == :down
        plot!(inputy, inputx, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor, orientation=:v)
      else
        plot!(inputx, inputy, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor)
      end
    end
    plot!() # flushing the display.
end

end # End module

max_stars_repo_name: ShozenD/AutocorrelationShell

max_issues_repo_name: ShozenD/AutocorrelationShell

max_forks_repo_name: ShozenD/AutocorrelationShell


text:
include("distribution.jl")
include("const.jl")

export entropy

function entropy_ML(data::Vector{Int64}, base::Number)
  p = fe1p(data)
  return -sum([ p[x] > 0 ? (p[x] * log(base, p[x])) : 0 for x=1:size(p)[1]])
end

# implemented from [1] (see below)
function entropy_MLBC(data::Vector{Int64}, base::Number)
  p = fe1p(data)
  n = float(size(data)[1])
  S = float(size(p)[1])
  H = -sum([ p[x] > ϵ ? (p[x] * log(base, p[x])) : 0 for x=1:size(p)[1]])
  return H + (S-1) / (2.0 * n)
end

# implemented from [1] (see below)
function entropy_HT(data::Vector{Int64}, base::Number)
  p = fe1p(data)
  n = size(data)[1]
  return -sum([ p[x] > ϵ ? ((p[x] * log(base, p[x])) / (1.0 - ((1.0 - p[x])^n))) : 0 for x=1:size(p)[1]])
end

# implemented from [1] (see below)
function entropy_CS(data::Vector{Int64}, base::Number)
  m = maximum(v)
  n  = size(v)[1]
  c  = counts(v, 1:m)
  c = c ./ n
  # just to get rid of the numerical inaccuracies and make sure its a probability distribution
  s = sum(c)
  p = c ./ s
  C = 1.0 - float(sum(filter(x == 1, c))) / float(n)
  p = p .* C
  return -sum([ p[x] > ϵ ? ((p[x] * log(base, p[x])) / (1.0 - ((1.0 - p[x])^l))) : 0 for x=1:size(p)[1]])
end

function entropy(data::Vector{Int64}; base=2, mode="ML")
  modes  = ["ML", "Maximum Likelihood",
            "MLBC", "Maximum Likelihood with Bias Correction",
            "Horovitz-Thompson", "HT",
            "ChaoShen", "Chao-Shen", "CS"]
  umodes = map(x->uppercase(x), modes)
  known_mode = uppercase(mode) in umodes
  pmodes = map(x->", $x",modes)
  pmodes = foldl(*, pmodes[1][3:end], pmodes[2:end])
  @assert known_mode "Mode may be any of the following: [$pmodes]."

  if uppercase(mode) in umodes[1:2]
    return entropy_ML(data, base)
  elseif uppercase(mode) in umodes[3:4]
    return entropy_MLBC(data, base)
  elseif uppercase(mode) in umodes[5:6]
    return entropy_HT(data, base)
  elseif uppercase(mode) in umodes[7:9]
    return entropy_CS(data, base)
  end
  return nothing
end

# [1] A. Chao and T.-J. Shen. Nonparametric estimation of shannon’s index of diversity when there are unseen species in sample. Environmental and Ecological Statistics, 10(4):429–443, 2003.


text:
using JuMP
using Clp
using LinearAlgebra
using Plots

model = Model(Clp.Optimizer)

T = 100
d = 50
f = 1000
δ = 1e-3
m = 0
p = zeros(T)
q = zeros(T)

@variables(model, begin
        m >= 0
        p[1:T] >= 0
        q[1:T] >= 0
    end
)

@objective(model, Min, m)

@constraints(
	model,
	begin
        m .>= (p + q)
		sum(p) - sum(q) == 0
		sum((T - t) * (p[t] - q[t]) for t in 1:T) == d
        sum(p) + sum(q) <= f
        [t=1:T-1], -δ <= ((p[t] - q[t]) - (p[t + 1] - q[t + 1])) <= δ
    end
)

optimize!(model)

a = value.(p - q)
v = zeros(T + 1)
x = zeros(T + 1)
for t in 1:T
    v[t + 1] = v[t] + a[t]
    x[t + 1] = x[t] + v[t]
end

plot(1:T, a, label="a")
savefig("a.png")
plot(1:T+1, v, label="v")
savefig("v.png")
plot(1:T+1, x, label="x")
savefig("x.png")

max_stars_repo_name: tansongchen/learn-optimization

max_issues_repo_name: tansongchen/learn-optimization

max_forks_repo_name: tansongchen/learn-optimization


text:
#TODO really needed ? It seems that we can return the given space if fix_constant == false
abstract type ConstantApproach end
struct FixConstant      <: ConstantApproach end
struct DoNotFixConstant <: ConstantApproach end

struct FESpaceWithConstantFixed{CA<:ConstantApproach,S<:SingleFieldFESpace} <: SingleFieldFESpace
  space::S
  dof_to_fix::Int
  @doc """
      FESpaceWithConstantFixed(space::SingleFieldFESpace, fix_constant::Bool,
      dof_to_fix::Int=num_free_dofs(space))
  """
  function FESpaceWithConstantFixed(
    space::SingleFieldFESpace,
    fix_constant::Bool,
    dof_to_fix::Int=num_free_dofs(space))

    S = typeof(space)
    if (fix_constant && num_dirichlet_dofs(space)==0)
      new{FixConstant,S}(space,dof_to_fix)
    else
      new{DoNotFixConstant,S}(space,dof_to_fix)
    end
  end
end

ConstantApproach(::FESpaceWithConstantFixed{CA}) where CA = CA()

# Genuine functions
function num_free_dofs(f::FESpaceWithConstantFixed{FixConstant})
  num_free_dofs(f.space)-1
end

function num_free_dofs(f::FESpaceWithConstantFixed{DoNotFixConstant})
  num_free_dofs(f.space)
end

function get_cell_dof_ids(f::FESpaceWithConstantFixed{FixConstant})
  cell_dofs = get_cell_dof_ids(f.space)
  CellDofIdsWithDofFixed(cell_dofs,f.dof_to_fix)
end

function get_cell_dof_ids(f::FESpaceWithConstantFixed{DoNotFixConstant})
  get_cell_dof_ids(f.space)
end

num_dirichlet_dofs(f::FESpaceWithConstantFixed{FixConstant}) = 1

num_dirichlet_dofs(f::FESpaceWithConstantFixed{DoNotFixConstant}) = 0

num_dirichlet_tags(f::FESpaceWithConstantFixed{FixConstant}) = 1

num_dirichlet_tags(f::FESpaceWithConstantFixed{DoNotFixConstant}) = 0

get_dirichlet_dof_tag(f::FESpaceWithConstantFixed{FixConstant}) = Int8[1,]

get_dirichlet_dof_tag(f::FESpaceWithConstantFixed{DoNotFixConstant}) = Int8[]

function scatter_free_and_dirichlet_values(f::FESpaceWithConstantFixed{FixConstant},fv,dv)
  @assert length(dv) == 1
  _dv = similar(dv,eltype(dv),0)
  _fv = VectorWithEntryInserted(fv,f.dof_to_fix,dv[1])
  scatter_free_and_dirichlet_values(f.space,_fv,_dv)
end

function scatter_free_and_dirichlet_values(f::FESpaceWithConstantFixed{DoNotFixConstant},fv,dv)
  @assert length(dv) == 0
 scatter_free_and_dirichlet_values(f.space,fv,dv)
end

function gather_free_and_dirichlet_values(f::FESpaceWithConstantFixed{FixConstant},cv)
  _fv, _dv = gather_free_and_dirichlet_values(f.space,cv)
  @assert length(_dv) == 0
  fv = VectorWithEntryRemoved(_fv,f.dof_to_fix)
  dv = _fv[f.dof_to_fix:f.dof_to_fix]
  (fv, dv)
end

function gather_free_and_dirichlet_values(f::FESpaceWithConstantFixed{DoNotFixConstant},cv)
  gather_free_and_dirichlet_values(f.space,cv)
end

function gather_free_and_dirichlet_values!(fv,dv,f::FESpaceWithConstantFixed{FixConstant},cv)
  _fv, _dv = gather_free_and_dirichlet_values(f.space,cv)
  @assert length(_dv) == 0
  fv    .= VectorWithEntryRemoved(_fv,f.dof_to_fix)
  dv[1]  = _fv[f.dof_to_fix]
  (fv, dv)
end

function gather_free_and_dirichlet_values!(fv,dv,f::FESpaceWithConstantFixed{DoNotFixConstant},cv)
  gather_free_and_dirichlet_values(f.space,cv)
end

function TrialFESpace(f::FESpaceWithConstantFixed{CA}) where CA
  U = TrialFESpace(f.space)
  FESpaceWithConstantFixed(U,CA==FixConstant,f.dof_to_fix)
end

# Delegated functions

function get_cell_shapefuns(f::FESpaceWithConstantFixed)
  get_cell_shapefuns(f.space)
end

function get_cell_shapefuns_trial(f::FESpaceWithConstantFixed)
  get_cell_shapefuns_trial(f.space)
end

function get_cell_dof_basis(f::FESpaceWithConstantFixed)
  get_cell_dof_basis(f.space)
end

get_vector_type(f::FESpaceWithConstantFixed) = get_vector_type(f.space)

get_dof_value_type(f::FESpaceWithConstantFixed) = get_dof_value_type(f.space)

CellField(t::FESpaceWithConstantFixed,cell_vals) = CellField(t.space,cell_vals)

ConstraintStyle(::Type{<:FESpaceWithConstantFixed{CA,S}}) where {CA,S} = ConstraintStyle(S)

get_triangulation(f::FESpaceWithConstantFixed) = get_triangulation(f.space)

# Helpers

struct CellDofIdsWithDofFixed{A<:AbstractArray} <: AbstractVector{Vector{Int32}}
  cell_dofs::A
  dof_to_fix::Int
end

Base.size(a::CellDofIdsWithDofFixed) = (length(a.cell_dofs),)

Base.IndexStyle(::Type{<:CellDofIdsWithDofFixed}) = IndexLinear()

function Base.getindex(a::CellDofIdsWithDofFixed,i::Integer)
  cache = array_cache(a)
  getindex!(cache,a,i)
end

function array_cache(a::CellDofIdsWithDofFixed)
  @assert eltype(a.cell_dofs) == Vector{Int32}
  b = testitem(a.cell_dofs)
  c = CachedArray(b)
  cache = array_cache(a.cell_dofs)
  (c, cache)
end

@inline function getindex!(d,a::CellDofIdsWithDofFixed,i::Integer)
  c, cache = d
  b = getindex!(cache,a.cell_dofs,i)
  setsize!(c,size(b))
  r = c.array
  for j in 1:length(b)
    bj = b[j]
    if bj == a.dof_to_fix
      rj = -1
    elseif bj < a.dof_to_fix
      rj = bj
    else
      rj = bj-1
    end
    r[j] = Int32(rj)
  end
  r
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9b5db9f2-3974-11ec-28f2-65388bfe0fbf
md"""
# The secrets (`config/secrets.jl`) file

Confidential configuration data (like API keys, usernames, passwords, etc) should be added to the `config/secrets.jl` file. This file is by default added to `.gitignore` when creating a Genie app, so it won't be added to source control -- avoid that it is accidentally exposed.

## Scope

All the definitions (variables, constants, functions, modules, etc) added to the `secrets.jl` file are loaded into your app's module. So if your app is called `MyGenieApp`, the definitions will be available under the `MyGenieApp` module.

---
**HEADS UP**

Given the your app's name is variable, you can also access your app's module through the `UserApp` constant. So all the definitions added to `secrets.jl` can also be accessed through the `UserApp` module (`UserApp === MyGenieApp`).

---
"""

# ╔═╡ Cell order:
# ╟─9b5db9f2-3974-11ec-28f2-65388bfe0fbf


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: cost_and_usage_report_service
using AWS.Compat
using AWS.UUIDs

"""
    delete_report_definition()
    delete_report_definition(params::Dict{String,<:Any})

Deletes the specified report.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ReportName"`: The name of the report that you want to delete. The name must be unique,
  is case sensitive, and can't include spaces.
"""
function delete_report_definition(; aws_config::AbstractAWSConfig=global_aws_config())
    return cost_and_usage_report_service(
        "DeleteReportDefinition"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function delete_report_definition(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return cost_and_usage_report_service(
        "DeleteReportDefinition",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_report_definitions()
    describe_report_definitions(params::Dict{String,<:Any})

Lists the AWS Cost and Usage reports available to this account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`:
- `"NextToken"`:
"""
function describe_report_definitions(; aws_config::AbstractAWSConfig=global_aws_config())
    return cost_and_usage_report_service(
        "DescribeReportDefinitions"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_report_definitions(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return cost_and_usage_report_service(
        "DescribeReportDefinitions",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    modify_report_definition(report_definition, report_name)
    modify_report_definition(report_definition, report_name, params::Dict{String,<:Any})

Allows you to programatically update your report preferences.

# Arguments
- `report_definition`:
- `report_name`:

"""
function modify_report_definition(
    ReportDefinition, ReportName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return cost_and_usage_report_service(
        "ModifyReportDefinition",
        Dict{String,Any}(
            "ReportDefinition" => ReportDefinition, "ReportName" => ReportName
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function modify_report_definition(
    ReportDefinition,
    ReportName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return cost_and_usage_report_service(
        "ModifyReportDefinition",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ReportDefinition" => ReportDefinition, "ReportName" => ReportName
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_report_definition(report_definition)
    put_report_definition(report_definition, params::Dict{String,<:Any})

Creates a new report using the description that you provide.

# Arguments
- `report_definition`: Represents the output of the PutReportDefinition operation. The
  content consists of the detailed metadata and data file information.

"""
function put_report_definition(
    ReportDefinition; aws_config::AbstractAWSConfig=global_aws_config()
)
    return cost_and_usage_report_service(
        "PutReportDefinition",
        Dict{String,Any}("ReportDefinition" => ReportDefinition);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_report_definition(
    ReportDefinition,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return cost_and_usage_report_service(
        "PutReportDefinition",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("ReportDefinition" => ReportDefinition), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
"""
    get_dummy(df, cVar, cType. trtRef)

Convert categorical variable to dummy indicators using specified contrast
type. This covers all cases except for treatment contrasts with a specified
reference level.

# Arguments

- df = DataFrame of variables
- cVar = symbol for the categorical variable in df to be converted
- cType = character string indicating the type of contrast to use for `cVar`
- trtRef = nothing

# Value

DataFrame of dummy variables for the specified categorical variable

"""
function get_dummy(df::DataFrames.DataFrame, cVar::Symbol, cType::String,
                   trtRef::Nothing)
    # Obtain the levels to use for the dummy indicators, depending on
    # contrast type
    thisVar = string.(df[:,cVar])
    if cType=="treat"
        levs = unique(thisVar)[2:end]
    elseif (cType=="sum")
        levs = unique(thisVar)[(2:end).-1]
    elseif cType=="noint" || (cType=="sumnoint")
        levs = unique(thisVar)
    else
        error(string("Did not recognize contrast type for ", cVar))
    end

    # Iterate through levels to make dummy indicators
    namedict = Dict(zip(levs, 1:length(levs)))
    dummies = zeros(length(thisVar), length(namedict))
    for i=1:length(thisVar)
        if haskey(namedict, thisVar[i])
            dummies[i, namedict[thisVar[i]]] = 1
        end
    end

    # Some additional modifications for sum contrasts
    if (cType=="sum")
        dummies[thisVar.==(unique(thisVar)[end]),:] .= -1
    end
    if cType=="sumnoint"
        dummies = dummies - (1/length(unique(thisVar)))
    end

    # Convert results to a DataFrame and rename columns
    newDf = convert(DataFrame, dummies)
    rename!(newDf, [Symbol("$(cVar)_$k") for k in levs])
    return newDf
end

"""

    get_dummy(df, cVar, cType. trtRef)

Convert categorical variables to for treatment contrasts with a specified
reference level.

# Arguments

- df = DataFrame of variables
- cVar = symbol for the categorical variable in df to be converted
- cType = character string indicating the type of contrast to use for `cVar`
- trtRef = character string specifying the level in cVar to use as the
reference

# Value

DataFrame of dummy variables for the specified categorical variable

"""
function get_dummy(df::DataFrames.DataFrame, cVar::Symbol, cType::String,
                   trtRef::String)

    # Obtain the levels to use for the dummy indicators.
    thisVar = string.(df[:,cVar])
    if cType=="treat"
        levs = unique(thisVar)[unique(thisVar) .!= trtRef]
    else
        error("Can only specify trtRef for treatment contrasts.")
    end

    # Iterate through levels to make dummy indicators
    namedict = Dict(zip(levs, 1:length(levs)))
    dummies = zeros(length(thisVar), length(namedict))
    for i=1:length(thisVar)
        if haskey(namedict, thisVar[i])
            dummies[i, namedict[thisVar[i]]] = 1
        end
    end

    # Convert results to a DataFrame and rename columns
    newDf = convert(DataFrame, dummies)
    rename!(newDf, [Symbol("$(cVar)_$k") for k in levs])
    return newDf
end

"""
    contr(df, cVars, cTypes, trtRefs)

Converts categorical variables in a DataFrame to specified contrast types.
All other variables are left as-is.

# Arguments

- df = DataFrame of variables
- cVars = 1d array of symbols corresponding to categorical variable names
  in df to be converted
- cTypes = 1d array of character strings of the same length as `cVars`,
  indicating the types of contrasts to use. Defaults to treatment contrasts
  ("treat") for all variables in `cVars`. Other options include "sum" for sum
  contrasts, "noint" for treatment contrasts with no intercept, and
  "sumnoint" for sum contrasts with no intercept. For "treat" `cTypes`, you
  can also specify the level to use as the reference treatment using `trtRefs`.
- trtRefs = optional 1d array of character strings of the same length as
  `cVars`, specifying the level to use as the references for treatment
  contrasts. Defaults to nothing for all variables in `cVars`.

# Value

DataFrame with same variables as the original DataFrame, but categorical
variables converted to dummy contrasts.

# Some notes

If `cVars` consists of only an empty Symbol, i.e. `cVars=[Symbol()]`, this
will signal to the function that no contrasts should be created. The
original DataFrame will be returned.

"""
function contr(df::DataFrames.DataFrame, cVars::AbstractArray{Symbol,1},
               cTypes::AbstractArray{String,1}=repeat(["treat"],
                                                      inner=length(cVars)),
               trtRefs::AbstractArray= repeat([nothing], inner=length(cVars)))

    # If cVars only contains an empty Symbol, stop and return df
    if cVars == [Symbol()]
        return df
    end

    # Ensure that all Symbols in cVars are variables in df
    for cVar in cVars
        if !in(cVar, propertynames(df))
            error(string(cVar, " is not a variable in the DataFrame"))
        end
    end

    # Initialize new DataFrame
    newDf = DataFrame()

    # Iterate through variables in df
    for var in propertynames(df)
        if !in(var, cVars)
            # Add non-categorical variables to the new DataFrame
            newDf[!,var] = df[:,var]
        else
            # Convert categorical variables to specified dummy contrasts
            dummyDf = get_dummy(df, var, cTypes[var.==cVars][1],
                                trtRefs[var.==cVars][1])
            for dummy in propertynames(dummyDf)
                newDf[!,dummy] = dummyDf[:,dummy]
            end
        end
    end

    return newDf
end


text: dummy_func(x,y) = 2x+y


text:
"""
    FactorType

`fixed` `manipulated` - A crossed fixed-effects factor

`random` `measured` - A crossed random-effects factor

`nested` - A random factor fully nested within another factor

`subject` `block` - A random factor subjected to multiple levels of another factor
"""
@enum FactorType fixed manipulated random measured nested subject block

Broadcast.broadcastable(a::FactorType) = (a,) # workaround for current behavior

isnested(x) = x == nested
isfixed(x) = x == fixed || x == manipulated
israndom(x) = x == random || x == measured
issubject(x) = x == subject || x == block


text:
function my_f(x,y)
    2x+y
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TOML

module Internals
    # The parser is defined in Base
    using Base.TOML: Parser, parse, tryparse, ParserError, isvalid_barekey_char, reinit!
    # Put the error instances in this module
    for errtype in instances(Base.TOML.ErrorType)
        @eval using Base.TOML: $(Symbol(errtype))
    end
    # We put the printing functionality in a separate module since It
    # defines a function `print` and we don't want that to collide with normal
    # usage of `(Base.)print` in other files
    module Printer
        include("print.jl")
    end
end

# https://github.com/JuliaLang/julia/issues/36605
readstring(f::AbstractString) = isfile(f) ? read(f, String) : error(repr(f), ": No such file")

"""
    Parser()

Constructor for a TOML `Parser`.  Note that in most cases one does not need to
explicitly create a `Parser` but instead one directly use use
[`TOML.parsefile`](@ref) or [`TOML.parse`](@ref).  Using an explicit parser
will however reuse some internal data structures which can be beneficial for
performance if a larger number of small files are parsed.
"""
const Parser = Internals.Parser

"""
    parsefile(f::AbstractString)
    parsefile(p::Parser, f::AbstractString)

Parse file `f` and return the resulting table (dictionary). Throw a
[`ParserError`](@ref) upon failure.

See also: [`TOML.tryparsefile`](@ref)
"""
parsefile(f::AbstractString) =
    Internals.parse(Parser(readstring(f); filepath=abspath(f)))
parsefile(p::Parser, f::AbstractString) =
    Internals.parse(Internals.reinit!(p, readstring(f); filepath=abspath(f)))

"""
    tryparsefile(f::AbstractString)
    tryparsefile(p::Parser, f::AbstractString)

Parse file `f` and return the resulting table (dictionary). Return a
[`ParserError`](@ref) upon failure.

See also: [`TOML.parsefile`](@ref)
"""
tryparsefile(f::AbstractString) =
    Internals.tryparse(Parser(readstring(f); filepath=abspath(f)))
tryparsefile(p::Parser, f::AbstractString) =
    Internals.tryparse(Internals.reinit!(p, readstring(f); filepath=abspath(f)))

"""
    parse(x::Union{AbstractString, IO})
    parse(p::Parser, x::Union{AbstractString, IO})

Parse the string  or stream `x`, and return the resulting table (dictionary).
Throw a [`ParserError`](@ref) upon failure.

See also: [`TOML.tryparse`](@ref)
"""
parse(str::AbstractString) =
    Internals.parse(Parser(String(str)))
parse(p::Parser, str::AbstractString) =
    Internals.parse(Internals.reinit!(p, String(str)))
parse(io::IO) = parse(read(io, String))
parse(p::Parser, io::IO) = parse(p, read(io, String))

"""
    tryparse(x::Union{AbstractString, IO})
    tryparse(p::Parser, x::Union{AbstractString, IO})

Parse the string or stream `x`, and return the resulting table (dictionary).
Return a [`ParserError`](@ref) upon failure.

See also: [`TOML.parse`](@ref)
"""
tryparse(str::AbstractString) =
    Internals.tryparse(Parser(String(str)))
tryparse(p::Parser, str::AbstractString) =
    Internals.tryparse(Internals.reinit!(p, String(str)))
tryparse(io::IO) = tryparse(read(io, String))
tryparse(p::Parser, io::IO) = tryparse(p, read(io, String))

"""
    ParserError

Type that is returned from [`tryparse`](@ref) and [`tryparsefile`](@ref)
when parsing fails. It contains (among others) the following fields:

- `pos`, the position in the string when the error happened
- `table`, the result that so far was successfully parsed
- `type`, an error type, different for different types of errors
"""
const ParserError = Internals.ParserError

"""
    print([to_toml::Function], io::IO [=stdout], data::AbstractDict; sorted=false, by=identity)

Write `data` as TOML syntax to the stream `io`. If the keyword argument `sorted` is set to `true`,
sort tables according to the function given by the keyword argument `by`.

The following data types are supported: `AbstractDict`, `Integer`, `AbstractFloat`, `Bool`,
`Dates.DateTime`, `Dates.Time`, `Dates.Date`. Note that the integers and floats
need to be convertible to `Float64` and `Int64` respectively. For other data types,
pass the function `to_toml` that takes the data types and returns a value of a
supported type.
"""
const print = Internals.Printer.print

end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
export i_mean_sum, i_var_mean_sum, i_normal_logpdf, i_cor_cov
export VarianceInfo

"""
    i_mean_sum(out!, sum!, x)

get the `mean` and `sum` of `x`.
"""
@i function i_mean_sum(out!, sum!, x)
    for i=1:length(x)
        sum! += x[i]
    end
    out! += sum!/(@const length(x))
end

struct VarianceInfo{T}
    variance::T
    variance_accumulated::T
    mean::T
    sum::T
end

function VarianceInfo(::Type{T}) where T
    VarianceInfo(zero(T), zero(T), zero(T), zero(T))
end

"""
    i_var_mean_sum(varinfo, sqv)
    i_var_mean_sum(var!, varsum!, mean!, sum!, v)

Compute the variance, the accumulated variance, mean and sum.
`varinfo` is the `VarianceInfo` object to store outputs.
"""
@i function i_var_mean_sum(varinfo::VarianceInfo{T}, v::AbstractVector{T}) where T
    i_var_mean_sum(varinfo.variance, varinfo.variance_accumulated, varinfo.mean, varinfo.sum, v)
end

@i function i_var_mean_sum(var!, varsum!, mean!, sum!, v::AbstractVector{T}) where T
    i_mean_sum(mean!, sum!, v)
     for i=1:length(v)
        @routine @invcheckoff begin
            x ← zero(T)
            x += v[i] - mean!
        end
        varsum! += x ^ 2
        ~@routine
     end
    var! += varsum! / (@const length(v)-1)
 end

"""
    i_normal_logpdf(out, x, μ, σ)

get the pdf of `Normal(μ, σ)` at point `x`.
"""
@i function i_normal_logpdf(out, x::T, μ, σ) where T
    @zeros T anc1 anc2 anc3

    @routine begin
        anc1 += x
        anc1 -= μ
        anc2 += anc1 / σ  # (x- μ)/σ
        anc3 += anc2^2 # (x-μ)^2/σ^2
    end

    out -= anc3 * 0.5 # -(x-μ)^2/2σ^2
    out -= log(σ) # -(x-μ)^2/2σ^2 - log(σ)
    out -= log(2π)/2 # -(x-μ)^2/2σ^2 - log(σ) - log(2π)/2

    ~@routine
end

"""
     i_cor_cov(rho!,cov!,a,b)

get Pearson correlation and covariance of two vectors `a` and `b`

"""

@i function i_cor_cov(rho!::T, cov!::T, a::AbstractVector{T}, b::AbstractVector{T}) where T
    @safe @assert length(a) == length(b)
    @routine  @invcheckoff begin
        @zeros T std1 std2
        info1 ← _zero(VarianceInfo{T})
        i_var_mean_sum(info1, a)
        std1 += sqrt(info1.variance)
        info2 ← _zero(VarianceInfo{T})
        i_var_mean_sum(info2, b)
        std2 += sqrt(info2.variance)
        @zeros T anc5 anc6 anc7
        @inbounds for i=1:length(b)
            @routine begin
                @zeros T anc3 anc4
                anc3 += a[i] - info1.mean
                anc4 += b[i] - info2.mean
            end
            anc5 += anc3 * anc4
            ~@routine
        end
        anc6 += std1 * std2
        anc7 += anc6 * (@const length(b)-1)
    end
    cov! += anc5 / (@const length(b)-1)
    rho! += anc5 / anc7
    ~@routine
end


text:
import Observables
import Makie
import FileIO

# Map from S² into its upper hemisphere
s2tos2map(b::S²) = begin
    p = Geographic(b)
    r = sqrt((1 - sin(p.θ)) / 2)
    Geographic(p.r, r * cos(p.ϕ), r * sin(p.ϕ))
end

scene = Makie.Scene()
circle = U1(rand() * 2pi - pi)
s3rotation = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
config = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
segments = rand(5:10)
color = FileIO.load("../data/basemap_color.png")
transparency = rand(1:2) == 1 ? true : false
s2tos3map = rand() > 0.5 ? σmap : τmap
frame = Frame(scene,
              circle,
              s2tos3map,
              s2tos2map,
              color,
              s3rotation = s3rotation,
              config = config,
              segments = segments,
              transparency = transparency)

## update circle

circle2 = U1(rand() * 2pi - pi)
update(frame, circle2)

@test isapprox(frame.circle, circle2)

## update maps f: S² → S³ and g: S² → S²

s2tos3map = rand() > 0.5 ? σmap : τmap
s2tos2map′(x) = Geographic(0.5 * x.r, 0.5 * x.ϕ, 0.5 * x.θ)
update(frame, s2tos3map, s2tos2map′)
p = Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2)

@test isapprox(frame.s2tos3map(p), s2tos3map(p))
@test isapprox(frame.s2tos2map(p), s2tos2map′(p))

## update S³ rotation

s3rotation2 = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
update(frame, s3rotation2)

@test isapprox(frame.s3rotation, s3rotation2)

## update configuration

config2 = config * Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
value1 = getsurface(frame.observable, frame.segments)
update(frame, config2)
value2 = getsurface(frame.observable, frame.segments)

@test isapprox(frame.config, config2)
@test isapprox(value1, value2) == false

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta