text:
#primary structure for portfolio math
#some auxillery structures to improve missing value efficiency
struct Portfolio
  U::Universe
  wₜ::Vector{Float64}
  R::Vector{Float64}
  pμ::Base.RefValue{Float64}

  #the below vector contains all valid weights, and does not include invalid weights
  #thus the shape of wvec is a ragged array, with each array of length equal to the
  #number of valid values
  wvec::Vector{Vector{Float64}} #store this so we don't have to reallocate
  wvecsize::Vector{Int} #holds the lengths so we don't ahve to call length
end

#update the portfolio object
#innermost loop and likely a performance bottleneck, so highly optimized
function Portfolio!(P, wₜ)::Portfolio

  P.wₜ .= wₜ
  #P.wmat .= P.U.missingreturns .* wₜ
  #P.wmat ./= sum(P.wmat, dims=1)

  for c ∈ 1:P.U.N
    for r ∈ 1:P.wvecsize[c]
      P.wvec[c][r] = wₜ[P.U.missingindices[c][r]]
    end

    adjfactor = sum(P.wvec[c])
    P.wvec[c] ./= adjfactor

  end

  #build the returns (Note: Threads.@threads provides signficant performance gains)
  P.R .= 0.0
  for c ∈ 1:P.U.N
    for r ∈ 1:P.wvecsize[c]
      P.R[c] += P.wvec[c][r] * P.U.returns[P.U.missingindices[c][r],c]
    end
  end

  #update the pointer
  P.pμ[] = mean(P.R)
  #print(P.R)
  return P
end

#construct a brand new portfolio object
function Portfolio(U, wₜ)::Portfolio
  wvecsize::Vector{Int} = (v->length(v)).(U.missingindices)
  wvec::Vector{Vector{Float64}} = (n->Vector{Float64}(undef, n)).(wvecsize)

  P::Portfolio = Portfolio( #pre-allocate if needed, reduces code reuse
    U,
    Vector{Float64}(undef, U.M),
    Vector{Float64}(undef, U.N),
    Ref(-1.0),
    #Matrix{Float64}(undef, U.M, U.N))
    wvec,
    wvecsize)

  return Portfolio!(P, wₜ)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import DelayEmbeddings: Dataset

"""
    cembed(data::Dataset)

Returns an embedding consisting of a zero-lagged, unmodified version of `data`.
"""
function cembed(data::Dataset)
    if size(data, 1) > size(data, 2)
        #info("Treating each row as a point")
    	dim = size(data, 2)
        which_pos = [i for i = 1:dim]
        which_lags = [0 for i in 1:dim]
        cembed([data[:, i] for i = 1:dim], which_pos, which_lags)
    else
    	#info("Treating each column of data as a point")
    	dim = size(data, 1)
        which_pos = [i for i = 1:dim]
        which_lags = [0 for i in 1:dim]
        cembed([data[i, :] for i = 1:dim], which_pos, which_lags)
    end
end

"""
    cembed(d::Dataset,
        which_pos::Vector{Int},
        which_lags::Vector{Int})

Returns a state space embedding of the columns of `d`.

## Arguments
* `d::Dataset`: The columns of `d` contains the data series to use
    for the embedding.
* `which_pos::Vector{Int}``. The length of in_which_pos gives the dimension
    of the embedding. The value of the ith element of `in_which_pos`
    indicates which column of `d` goes in the i-th column of the embedding.
* `which_lags::Vector{Int}` sets the lag in each column of the reconstrution.
    Must be the same length as `dimension(d)`
    - **Example**: if `in_which_pos = [2, 2, 1]`, then
        `at_what_lags = [1, 0, -1]` means that the lag in column 1 is 1, the
        lag in the second column is 0 and the lag in the third column is -1.
"""
function cembed(data::Dataset,
        which_pos,
        which_lags)

    if size(data, 1) > size(data, 2)
        #info("Treating each row as a point")
    	dim = size(data, 2)
    	cembed([data[:, i] for i = 1:dim], which_pos, which_lags)
    else
    	#info("Treating each column of data as a point")
    	dim = size(data, 1)
        cembed([data[i, :] for i = 1:dim], which_pos, which_lags)
    end
end

function Embedding(data::Dataset,
        which_pos,
        which_lags)
    if size(data, 1) > size(data, 2)
        #info("Treating each row as a point")
        dim = size(data, 2)
        cembed([data[:, i] for i = 1:dim], which_pos, what_lags)
    else
        #info("Treating each column of data as a point")
        dim = size(data, 1)
        cembed([data[i, :] for i = 1:dim], which_pos, what_lags)
    end
end

export cembed

ext: jl


text:
const ERFA_DPI = 3.141592653589793
const ERFA_D2PI = 6.283185307179586
const ERFA_DR2D = 57.29577951308232
const ERFA_DD2R = 0.017453292519943295
const ERFA_DR2AS = 206264.80624709636
const ERFA_DAS2R = 4.84813681109536e-6
const ERFA_DS2R = 7.27220521664304e-5
const ERFA_TURNAS = 1.296e6
const ERFA_DMAS2R = ERFA_DAS2R / 1000.0
const ERFA_DTY = 365.242198781
const ERFA_DAYSEC = 86400.0
const ERFA_DJY = 365.25
const ERFA_DJC = 36525.0
const ERFA_DJM = 365250.0
const ERFA_DJ00 = 2.451545e6
const ERFA_DJM0 = 2.4000005e6
const ERFA_DJM00 = 51544.5
const ERFA_DJM77 = 43144.0
const ERFA_TTMTAI = 32.184
const ERFA_DAU = 1.4959787e11
const ERFA_CMPS = 2.99792458e8
const ERFA_AULT = 499.004782
const ERFA_DC = ERFA_DAYSEC / ERFA_AULT
const ERFA_ELG = 6.969290134e-10
const ERFA_ELB = 1.550519768e-8
const ERFA_TDB0 = -6.55e-5
const ERFA_SRS = 1.97412574336e-8

# Skipping MacroDefinition: ERFA_DINT ( A ) ( ( A ) < 0.0 ? ceil ( A ) : floor ( A ) )
# Skipping MacroDefinition: ERFA_DNINT ( A ) ( ( A ) < 0.0 ? ceil ( ( A ) - 0.5 ) : floor ( ( A ) + 0.5 ) )
# Skipping MacroDefinition: ERFA_DSIGN ( A , B ) ( ( B ) < 0.0 ? - fabs ( A ) : fabs ( A ) )
# Skipping MacroDefinition: ERFA_GMAX ( A , B ) ( ( ( A ) > ( B ) ) ? ( A ) : ( B ) )
# Skipping MacroDefinition: ERFA_GMIN ( A , B ) ( ( ( A ) < ( B ) ) ? ( A ) : ( B ) )

const ERFA_WGS84 = 1
const ERFA_GRS80 = 2
const ERFA_WGS72 = 3

immutable Array_3_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
end

zero(::Type{Array_3_Cdouble}) = Array_3_Cdouble(fill(zero(Cdouble),3)...)

immutable Array_3_Array_3_Cdouble
    d1::Array_3_Cdouble
    d2::Array_3_Cdouble
    d3::Array_3_Cdouble
end

zero(::Type{Array_3_Array_3_Cdouble}) = Array_3_Array_3_Cdouble(fill(zero(Array_3_Cdouble),3)...)

type eraASTROM
    pmt::Cdouble
    eb::Array_3_Cdouble
    eh::Array_3_Cdouble
    em::Cdouble
    v::Array_3_Cdouble
    bm1::Cdouble
    bpn::Array_3_Array_3_Cdouble
    along::Cdouble
    phi::Cdouble
    xpl::Cdouble
    ypl::Cdouble
    sphi::Cdouble
    cphi::Cdouble
    diurab::Cdouble
    eral::Cdouble
    refa::Cdouble
    refb::Cdouble
end

immutable Array_2_Array_3_Cdouble
    d1::Array_3_Cdouble
    d2::Array_3_Cdouble
end

zero(::Type{Array_2_Array_3_Cdouble}) = Array_2_Array_3_Cdouble(fill(zero(Array_3_Cdouble),2)...)

immutable eraLDBODY
    bm::Cdouble
    dl::Cdouble
    pv::Array_2_Array_3_Cdouble
end

immutable Array_4_Cint
    d1::Cint
    d2::Cint
    d3::Cint
    d4::Cint
end

zero(::Type{Array_4_Cint}) = Array_4_Cint(fill(zero(Cint),4)...)

immutable Array_2_Cdouble
    d1::Cdouble
    d2::Cdouble
end

zero(::Type{Array_2_Cdouble}) = Array_2_Cdouble(fill(zero(Cdouble),2)...)


text:
"""
    Rotational{T} <: Joint{T}

    constraint limiting rotational degrees of freedom

    axis: rotation axis in parent offset frame
    axis_mask1: rotation axis mask in pbody's frame
    axis_mask2: rotation axis mask in pbody's frame
    axis_mask3: rotation axis mask in pbody's frame
    axis_offset: rotation axis offset from pbody's frame
    spring :coefficient for joint spring
    damper: coefficient for joint damper
    spring_offset: nominal joint configuration
    joint_limits: lower and upper limits on the joint configuration
    spring_type: can be :linear or :sinusoidal (currently not implemented), if linear then we need joint_limits to avoid the 360° singularity.
    input: external impulse torque
"""
mutable struct Rotational{T,Nλ,Nb,N,Nb½,N̄λ} <: Joint{T,Nλ,Nb,N,Nb½}
    axis::SVector{3,T} # rotation axis in parent offset frame
    axis_mask1::Adjoint{T,SVector{3,T}} # in pbody's frame
    axis_mask2::Adjoint{T,SVector{3,T}} # in pbody's frame
    axis_mask3::Adjoint{T,SVector{3,T}} # in pbody's frame
    axis_offset::Quaternion{T} # in pbody's frame
    spring::T
    damper::T
    spring_offset::SVector{N̄λ,T}
    joint_limits::Vector{SVector{Nb½,T}} # lower and upper limits on the joint minimal coordinate angles
    spring_type::Symbol # the rotational springs can be :linear or :sinusoidal (currently not implemented), if linear then we need joint_limits to avoid the 360° singularity.
    input::SVector{3,T}
end

function Rotational{T,Nλ}(pbody::Node, cbody::Node;
        axis::AbstractVector=szeros(T,3),
        axis_offset::Quaternion=one(Quaternion{T}),
        spring=zero(T),
        damper=zero(T),
        spring_offset=szeros(T,3-Nλ),
        joint_limits=[szeros(T,0), szeros(T,0)],
        spring_type::Symbol=:linear) where {T,Nλ}
    V1, V2, V3 = orthogonal_rows(axis)
    input = zeros(T,3)
    Nb½ = length(joint_limits[1])
    Nb = 2Nb½
    N̄λ = 3 - Nλ
    N = Nλ + 2Nb
    Rotational{T,Nλ,Nb,N,Nb½,N̄λ}(axis, V1, V2, V3, axis_offset, spring, damper, spring_offset, joint_limits, spring_type, input), pbody.id, cbody.id
end


text:
"""
    frac(startdate::Date, enddate::Date, dc::DayCount)

Compute the fractional number of years between `startdate` and `enddate`, according to the
[`DayCount` object](@ref daycount_types) `dc`.

- If `startdate == enddate`, then the result is zero.
- If `startdate > enddate`, then the result is `-frac(enddate, startdate, dc)`.
"""
function frac
end

abstract type DayCount end

"""
    Actual365Fixed()

**Actual/365 (Fixed)** day count convention.

The year fraction is computed as:
```math
\\frac{\\text{# of days}}{365}
```

# Reference
 - 2006 ISDA definitions, §4.16 (d)
"""
struct Actual365Fixed <: DayCount end
const Actual365F = Actual365Fixed
function frac(startdate::Date, enddate::Date, ::Actual365Fixed)
    return Dates.value(enddate-startdate)/365
end

"""
    Actual360()

**Actual/360** day count convention.

The year fraction is computed as:
```math
\\frac{\\text{# of days}}{360}
```

# Reference
 - 2006 ISDA definitions, §4.16 (e)

"""
struct Actual360 <: DayCount end
function frac(startdate::Date, enddate::Date, ::Actual360)
    return Dates.value(enddate-startdate)/360
end

"""
    ActualActualISDA()

**Actual/Actual (ISDA)** day count convention.

The year fraction is computed as:
```math
\\frac{\\text{# of days in standard year}}{365} +
\\frac{\\text{# of days in leap year}}{366}
```

For the purposes of above, the start date is included and the end date is excluded.

# Reference
 - 2006 ISDA definitions, §4.16 (b)
"""
struct ActualActualISDA <: DayCount end
function frac(startdate::Date, enddate::Date, dc::ActualActualISDA)
    if startdate > enddate
        return -frac(enddate, startdate, dc)
    end
    startyear = year(startdate)
    endyear = year(enddate)
    if startyear == endyear
        return Dates.value(enddate - startdate) / daysinyear(startdate)
    else
        d1 = daysinyear(startyear)
        n1 = d1 - dayofyear(startdate) + 1
        d2 = daysinyear(endyear)
        n2 = dayofyear(enddate) - 1
        return n1/d1 + n2/d2 + (endyear - startyear - 1)
    end
end

"""
    ActualActualICMA(schedule::StepRange{Date,Month})

**Actual/Actual (ICMA)**, **Actual/Actual (ISMA)** or **ISMA-99** day count convention.

Note that this is dependent on the coupon or payment schedule of the underlying
security. This is provided via the `schedule` argument, currently only date ranges with
steps of months are supported.

The year fraction is computed as:
```math
\\sum_{\\text{schedule period}} \\frac{\\text{# of days in period}}{\\text{length of period} \\times \\text{periods per year}}
```

This ensures that:
* all days in a period are of equal length, and
* all periods are of equal length.

# Reference
 - 2006 ISDA definitions, §4.16 (c).
 - [ICMA Rule Book](https://www.isda.org/a/NIJEE/ICMA-Rule-Book-Rule-251-reproduced-by-permission-of-ICMA.pdf), Rule 251.1 (iii).
 - [EMU and Market Conventions: Recent Developments, ISDA - BS:9951.1](https://www.isda.org/a/AIJEE/1998-ISDA-memo-EMU-and-Market-Conventions-Recent-Developments.pdf), §4. The Actual/Actual Day Count Convention.
"""
struct ActualActualICMA{T<:AbstractVector{Date}} <: DayCount
    frequency::Int
    schedule::T
end
ActualActualICMA(schedule::StepRange{Date,Month}) = ActualActualICMA(Int(Month(12)/step(schedule)),schedule)
const ISMA99 = ActualActualICMA
### const ActualActualISMA = ActualActualICMA

function frac(startdate::Date, enddate::Date, dc::ActualActualICMA)
    if startdate > enddate
        return -frac(enddate,startdate,dc)
    end
    startdate < first(dc.schedule) || enddate > last(dc.schedule) && error("The start and end dates must be contained within the scheduled dates.")

    schedule = dc.schedule

    periodfrac = 0.

    # Fraction of the last period
    istart = findlast(x -> x < enddate, schedule)
    periodfrac += schedule[istart+1] == enddate ? 1. : Dates.value(enddate-schedule[istart])/Dates.value(schedule[istart+1]-schedule[istart])

    # Fraction of the first period
    inext = findfirst(x -> x >= startdate, schedule)
    periodfrac += schedule[inext] == startdate ? 0. : Dates.value(schedule[inext]-startdate)/Dates.value(schedule[inext]-schedule[inext-1])

    # Number of full period between the first and last periods
    periodfrac += istart-inext

    return periodfrac/dc.frequency
end

# helper function
thirty360(dy,dm,dd) = (360*dy + 30*dm + dd)/360

"""
    Thirty360()

**30/360** or **Bond Basis** day count convention.

The year fraction is computed as:
```math
\\frac{360 \\times (y_2 - y_1) + 30 \\times (m_2 - m_1) + (d_2 - d_1)}{360}
```
where
- ``y_1`` and ``y_2`` are the years of the start and end date, respectively.
- ``m_1`` and ``m_2`` are the months of the start and end date, respectively.
- ``d_1`` is the day of the month at the start date, unless it is 31, in which case it is
  30.
- ``d_2`` is the day of the month at the end date, unless it is 31 and ``d_1 ≥ 30``, in
  which case it is 30.

# Reference
 - 2006 ISDA definitions, §4.16 (f)
"""
struct Thirty360 <: DayCount end
const BondBasis = Thirty360
function frac(startdate::Date, enddate::Date, dc::Thirty360)
    if startdate > enddate
        return -frac(enddate, startdate, dc)
    end

    dy = year(enddate)-year(startdate)
    dm = month(enddate)-month(startdate)

    d1 = day(startdate)
    d2 = day(enddate)
    if d1 >= 30
        d1 = 30
        if d2 >= 30
            d2 = 30
        end
    end
    return thirty360(dy,dm,d2-d1)
end

"""
    ThirtyE360()

**30E/360** or **Eurobond Basis** day count convention.

The year fraction is computed as:
```math
\\frac{360 \\times (y_2 - y_1) + 30 \\times (m_2 - m_1) + (d_2 - d_1)}{360}
```
where
- ``y_1`` and ``y_2`` are the years of the start and end date, respectively.
- ``m_1`` and ``m_2`` are the months of the start and end date, respectively.
- ``d_1`` is the day of the month at the start date, unless it is 31st day of the month,
  in which case it is 30.
- ``d_2`` is the day of the month at the end date,  unless it is 31st day of the month,
  in which case it is 30.

# Reference
 - 2006 ISDA definitions, §4.16 (g)
"""
struct ThirtyE360 <: DayCount end
const EurobondBasis = ThirtyE360

function frac(startdate::Date, enddate::Date, dc::ThirtyE360)
    if startdate > enddate
        return -frac(enddate, startdate, dc)
    end
    dy = year(enddate)-year(startdate)
    dm = month(enddate)-month(startdate)
    d1 = min(day(startdate),30)
    d2 = min(day(enddate),30)
    return thirty360(dy,dm,d2-d1)
end

"""
    ThirtyE360ISDA(maturity::Date)

**30E/360 (ISDA)** day count convention. Note that this is dependant of the maturity date
of the underlying security, provided as the `maturity` argument.

The year fraction is computed as:
```math
\\frac{360 \\times (y_2 - y_1) + 30 \\times (m_2 - m_1) + (d_2 - d_1)}{360}
```
where
- ``y_1`` and ``y_2`` are the years of the start and end date, respectively.
- ``m_1`` and ``m_2`` are the months of the start and end date, respectively.
- ``d_1`` is the day of the month at the start date, unless it is:
  * the last day of February, or
  * the 31st day of the month,
  in which case it is 30.
- ``d_2`` is the day of the month at the end date, unless it is:
  * the last day of February and not the maturity date, or
  * the 31st day of the month,
  in which case it is 30.

# Reference
 - 2006 ISDA definitions, §4.16 (h)
"""
struct ThirtyE360ISDA <: DayCount
    maturity::Date
end
function frac(startdate::Date, enddate::Date, dc::ThirtyE360ISDA)
    if startdate > enddate
        return -frac(enddate, startdate, dc)
    elseif startdate == enddate
        # in case startdate == enddate == maturitydate == end of Feb.
        return 0.0
    end
    y1 = year(startdate)
    y2 = year(enddate)
    m1 = month(startdate)
    m2 = month(enddate)
    d1 = startdate == lastdayofmonth(startdate) ? 30 : day(startdate)
    d2 = enddate == lastdayofmonth(enddate) && !(m2 == 2 && enddate == dc.maturity) ? 30 : day(enddate)
    return thirty360(y2-y1,m2-m1,d2-d1)
end


text:
const DEFAULT_MAX_MEMORY = 50 * 1024 * 1024
const DEFAULT_MAX_LATENCY = 2 * 60 * 1000
const DEFAULT_TIMEOUT = typemax(Int64)
const DEFAULT_MAX_WRITE_THREADS = 3

function _col_delete(col_family, col_qualifier)
    thriftbuild(ColumnUpdate, Dict(:colFamily => bytes(col_family),
        :colQualifier => bytes(col_qualifier),
        :deleteCell => true))
end

function _col_update(col_family, col_qualifier, col_visibility, value, timestamp::Integer=time_ms())
    thriftbuild(ColumnUpdate, Dict(:colFamily => bytes(col_family),
        :colQualifier => bytes(col_qualifier),
        :colVisibility => bytes(col_visibility),
        :timestamp => timestamp,
        :value => bytes(value)))
end

const ColumnUpdates         = Vector{ColumnUpdate}
const AbstractUpdate        = Union{ColumnUpdates, ConditionalUpdates}
const AbstractWriterOptions = Union{WriterOptions, ConditionalWriterOptions}

##
# BatchUpdates collect mutations to be applied with a batch operation
type BatchUpdates{T <: AbstractUpdate}
    mutations::Dict{Vector{UInt8}, T}
end
function (::Type{BatchUpdates{T}}){T}()
    BatchUpdates{T}(Dict{Vector{UInt8}, T}())
end

batch() = BatchUpdates{ColumnUpdates}()
conditional_batch() = BatchUpdates{ConditionalUpdates}()

function mutate{T<:ColumnUpdates}(upd::BatchUpdates{T}, row, action::ColumnUpdate)
    rowbytes = bytes(row)
    if rowbytes in keys(upd.mutations)
        mut = upd.mutations[rowbytes]
    else
        mut = upd.mutations[rowbytes] = ColumnUpdate[]
    end
    push!(mut, action)
    mut
end

function mutate{T<:ConditionalUpdates}(upd::BatchUpdates{T}, row, action::ColumnUpdate)
    rowbytes = bytes(row)
    if rowbytes in keys(upd.mutations)
        cu = upd.mutations[rowbytes]
    else
        cu = upd.mutations[rowbytes] = ConditionalUpdates()
    end
    mutate(cu, action)
end

function mutate(cu::ConditionalUpdates, condition::ColumnUpdate)
    isfilled(cu, :updates) || set_field!(cu, :updates, ColumnUpdate[])
    push!(cu.updates, condition)
    cu
end

function mutate{T<:ConditionalUpdates}(upd::BatchUpdates{T}, row, condition::Accumulo.proxy.Condition)
    rowbytes = bytes(row)
    if rowbytes in keys(upd.mutations)
        cu = upd.mutations[rowbytes]
    else
        cu = upd.mutations[rowbytes] = ConditionalUpdates()
    end
    mutate(cu, condition)
end

function mutate(cu::ConditionalUpdates, condition::Accumulo.proxy.Condition)
    isfilled(cu, :conditions) || set_field!(cu, :conditions, Accumulo.proxy.Condition[])
    push!(cu.conditions, condition)
    cu
end

# returns a Condition type
function where{T<:ConditionalUpdates}(upd::BatchUpdates{T}, row, col_family, col_qualifier, col_visibility=UInt8[]; value=nothing, timestamp::Integer=-1, iterators::Vector{IteratorSetting}=IteratorSetting[])
    col = thriftbuild(Column, Dict(:colFamily => bytes(col_family),
            :colQualifier => bytes(col_qualifier),
            :colVisibility => bytes(col_visibility)))

    cond = Accumulo.proxy.Condition()
    set_field!(cond, :column, col)

    (value === nothing) || set_field!(cond, :value, bytes(value))
    (timestamp < 0)     || set_field!(cond, :timestamp, Int64(timestamp))
    isempty(iterators)  || set_field!(cond, :iterators, iterators)

    mutate(upd, row, cond)
end

update{T}(upd::BatchUpdates{T}, row, col_family, col_qualifier, col_visibility, value, timestamp::Integer=time_ms()) = mutate(upd, row, _col_update(col_family, col_qualifier, col_visibility, value, timestamp))
delete{T}(upd::BatchUpdates{T}, row, col_family, col_qualifier) = mutate(upd, row, _col_delete(col_family, col_qualifier))
update(cu::ConditionalUpdates, col_family, col_qualifier, col_visibility, value, timestamp::Integer=time_ms()) = mutate(cu, _col_update(col_family, col_qualifier, col_visibility, value, timestamp))
delete(cu::ConditionalUpdates, col_family, col_qualifier) = mutate(cu, _col_delete(col_family, col_qualifier))

type BatchWriter{T <: AbstractWriterOptions}
    session::AccumuloSession
    tablename::String
    writername::String
    options::T
end

function batch_writer(session::AccumuloSession, tablename::String;
            max_memory::Integer=DEFAULT_MAX_MEMORY, latency_ms::Integer=DEFAULT_MAX_LATENCY,
            timeout_ms::Integer=DEFAULT_TIMEOUT, threads::Integer=DEFAULT_MAX_WRITE_THREADS,
            durability::Integer=Durability.DEFAULT)
    opts = thriftbuild(WriterOptions, Dict(:maxMemory => Int64(max_memory),
                :latencyMs => Int64(latency_ms),
                :timeoutMs => Int64(timeout_ms),
                :threads => Int32(threads),
                :durability => Int32(durability)))
    tbl = tablename
    writer = createWriter(client(session), handle(session), tbl, opts)
    BatchWriter(session, tbl, writer, opts)
end

function batch_writer(f::Function, args...; kwargs...)
    w = batch_writer(args...; kwargs...)
    try
        f(w)
    finally
        close(w)
    end
end

function conditional_batch_writer(session::AccumuloSession, tablename::String;
            max_memory::Integer=-1, timeout_ms::Integer=DEFAULT_TIMEOUT,
            threads::Integer=DEFAULT_MAX_WRITE_THREADS, authorizations::SET=[],
            durability::Integer=Durability.DEFAULT)
    opts = ConditionalWriterOptions()
    (max_memory < 0) || set_field!(opts, :maxMemory, Int64(max_memory))
    (timeout_ms < 0) || set_field!(opts, :timeoutMs, Int64(timeout_ms))
    (threads < 0) || set_field!(opts, :threads, Int32(threads))
    if !isempty(authorizations)
        auths = Set([bytes(x) for x in authorizations])
        set_field!(opts, :authorizations, auths)
    end
    (durability == Durability.DEFAULT) || set_field!(opts, :durability, Int32(durability))

    tbl = tablename
    writer = createConditionalWriter(client(session), handle(session), tbl, opts)
    BatchWriter(session, tbl, writer, opts)
end

function conditional_batch_writer(f::Function, args...; kwargs...)
    w = conditional_batch_writer(args...; kwargs...)
    try
        f(w)
    finally
        close(w)
    end
end

close(writer::BatchWriter{WriterOptions}) = closeWriter(client(writer.session), writer.writername)
close(writer::BatchWriter{ConditionalWriterOptions}) = closeConditionalWriter(client(writer.session), writer.writername)

flush(writer::BatchWriter{WriterOptions}) = flush(client(writer.session), writer.writername)
flush(writer::BatchWriter{ConditionalWriterOptions}) = nothing  # there's no flush API for conditional writes. this function just exists for uniformity

update(writer::BatchWriter{WriterOptions}, mutations::BatchUpdates{ColumnUpdates}) = update(client(writer.session), writer.writername, mutations.mutations)
update(writer::BatchWriter{ConditionalWriterOptions}, mutations::BatchUpdates{ConditionalUpdates}) = updateRowsConditionally(client(writer.session), writer.writername, mutations.mutations)

update(session::AccumuloSession, tablename::String, mutations::BatchUpdates{ColumnUpdates}) = updateAndFlush(client(session), handle(session), tablename, mutations.mutations)
function update(session::AccumuloSession, tablename::String, upd::BatchUpdates{ConditionalUpdates})
    results = Dict{Vector{UInt8}, ConditionalStatus}()
    for (row, cu) in upd.mutations
        result = updateRowConditionally(client(session), handle(session), tablename, row, cu)
        results[row] = result
    end
    results
end


text:
Dlist = (2,5,7,12)
@testset "PeriodicCMPS: environments and gauging with bond dimension $D" for D in Dlist
    for T in (Float64, ComplexF64)
        Q = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D)) for j=1:5])
        R = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D))/D for j=1:3])
        Ψ = InfiniteCMPS(Q, R)

        ρL, λL, infoL = leftenv(Ψ; Kmax = 20)
        TL = LeftTransfer(Ψ)
        @test norm(-∂(ρL) + TL(ρL) - 2*λL*ρL) <= 10*infoL.normres[1]
        @test ρL == ρL'

        ρR, λR, infoR = rightenv(Ψ; Kmax = 20)
        TR = RightTransfer(Ψ)
        @test norm(∂(ρR) + TR(ρR) - 2*λR*ρR) <= 10*infoR.normres[1]
        @test λL ≈ λR
        @test ρR == ρR'

        Zf = localdot(ρL, ρR)
        Z = Zf(0)
        @test Zf ≈ FourierSeries([Z])
        @test imag(Z) <= sqrt(infoL.normres[1]*infoR.normres[1])

        Ψn = copy(Ψ)
        ρL2, = leftenv!(Ψn; Kmax = 20)
        @test ρL2 ≈ ρL
        Qn = Ψn.Q
        Rn, = Ψn.Rs
        @test norm(-∂(ρL)+Qn'*ρL2 + ρL*Qn + Rn'*ρL*Rn) <= 1e-9

        Ψn = copy(Ψ)
        ρR2, = rightenv!(Ψn; Kmax = 20)
        @test ρR2 ≈ ρR
        Qn = Ψn.Q
        Rn, = Ψn.Rs
        @test norm(∂(ρR) + Qn*ρR + ρR*Qn' + Rn*ρR*Rn') <= 1e-9

        Ψn = copy(Ψ)
        ρL2, ρR2 = environments!(Ψn; Kmax = 20)
        @test ρL2 ≈ ρL/sqrt(Z)
        @test ρR2 ≈ ρR/sqrt(Z)
        Qn = Ψn.Q
        Rn, = Ψn.Rs
        @test norm(-∂(ρL) + LeftTransfer(Ψn)(ρL)) <= 1e-9
        @test norm(∂(ρR) + RightTransfer(Ψn)(ρR)) <= 1e-9
    end
end

@testset "PeriodicCMS: energy environments for bond dimension $D" for D in Dlist
    v = fit(x->-3+sin(x), FourierSeries; Kmax = 1)
    H = ∫(∂ψ̂'*∂ψ̂ + v*ψ̂'*ψ̂ + 0.3*(ψ̂*ψ̂ + ψ̂'*ψ̂') + 2.5*(ψ̂')^2*ψ̂^2, (-Inf,+Inf))
    for T in (Float64, ComplexF64)
        Q = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D)) for j=1:5])
        R = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D))/D for j=1:3])
        Ψ = InfiniteCMPS(Q, R)

        ρL, ρR = environments!(Ψ; Kmax = 20)
        HL, EL, eL, hL = leftenv(H, (Ψ,ρL,ρR); Kmax = 30)
        HR, ER, eR, hR = rightenv(H, (Ψ,ρL,ρR); Kmax = 30)

        @test eL ≈ eR
        @test EL ≈ ER
        @test norm(∂(HL) - LeftTransfer(Ψ)(HL) - hL) <= 1e-9*norm(HL)
        @test norm(∂(HR) + RightTransfer(Ψ)(HR) + hR) <= 1e-9*norm(HR)
        @test abs(dot(HL,ρR)) <= 1e-9*norm(HL)
        @test abs(dot(ρL,HR)) <= 1e-9*norm(HR)
    end
end

@testset "PeriodicCMS: local gradients with bond dimension $D" for D in Dlist
    for T in (Float64, ComplexF64)
        Q = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D)) for j=1:5])
        R1 = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D))/D for j=1:3])
        R2 = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D))/D for j=1:3])
        Rs = (R1, R2)
        Ψ = InfiniteCMPS(Q, Rs)
        ρL, ρR = environments!(Ψ; Kmax = 20)

        QR1 = Q*R1 - R1*Q + ∂(R1)
        QR2 = Q*R2 - R2*Q + ∂(R2)

        @test CMPSKit.localgradientQ(ψ̂[1], Q, Rs, ρL, ρR) == zero(Q)
        @test CMPSKit.localgradientQ(ψ̂[1]^2, Q, Rs, ρL, ρR) == zero(Q)
        @test CMPSKit.localgradientQ(∂ψ̂[1], Q, Rs, ρL, ρR) == zero(Q)
        @test CMPSKit.localgradientQ(ψ̂[2]', Q, Rs, ρL, ρR) == zero(Q)
        @test CMPSKit.localgradientQ((ψ̂[2]')^2, Q, Rs, ρL, ρR) == zero(Q)
        @test CMPSKit.localgradientQ(∂(ψ̂[2]'), Q, Rs, ρL, ρR) ≈ ρL*ρR*R2' - R2'*ρL*ρR
        @test CMPSKit.localgradientQ((∂ψ̂[1])'*∂ψ̂[2], Q, Rs, ρL, ρR) ≈
                ρL*QR2*ρR*R1' - R1'*ρL*QR2*ρR
        @test CMPSKit.localgradientQ(ψ̂[1]'*ψ̂[2], Q, Rs, ρL, ρR) == zero(Q)
        @test CMPSKit.localgradientQ((ψ̂[1]')^2*ψ̂[1]^2, Q, Rs, ρL, ρR) == zero(Q)

        @test CMPSKit.localgradientRs(ψ̂[1], Q, Rs, ρL, ρR) == (zero(Q), zero(Q))
        @test CMPSKit.localgradientRs(ψ̂[1]^2, Q, Rs, ρL, ρR) == (zero(Q), zero(Q))
        @test CMPSKit.localgradient∂Rs(ψ̂[1]^2, Q, Rs, ρL, ρR) == (zero(Q), zero(Q))

        @test CMPSKit.localgradientRs(∂ψ̂[1], Q, Rs, ρL, ρR) == (zero(Q), zero(Q))
        @test CMPSKit.localgradient∂Rs(∂ψ̂[1], Q, Rs, ρL, ρR) == (zero(Q), zero(Q))

        @test all(isapprox.(CMPSKit.localgradientRs(ψ̂[2]', Q, Rs, ρL, ρR),
                            (zero(Q), ρL*ρR)))
        @test all(isapprox.(CMPSKit.localgradientRs((ψ̂[2]')^2, Q, Rs, ρL, ρR),
                            (zero(Q), R2'*ρL*ρR + ρL*ρR*R2')))

        @test all(isapprox.(CMPSKit.localgradientRs(∂(ψ̂[2]'), Q, Rs, ρL, ρR),
                            (zero(Q), Q'*ρL*ρR - ρL*ρR*Q')))
        @test all(isapprox.(CMPSKit.localgradient∂Rs(∂(ψ̂[2]'), Q, Rs, ρL, ρR),
                            (zero(Q), ρL*ρR)))

        @test all(isapprox.(CMPSKit.localgradientRs((∂ψ̂[1])'*∂ψ̂[2], Q, Rs, ρL, ρR),
                            (Q'*ρL*QR2*ρR - ρL*QR2*ρR*Q', zero(Q))))
        @test all(isapprox.(CMPSKit.localgradient∂Rs((∂ψ̂[1])'*∂ψ̂[2], Q, Rs, ρL, ρR),
                            (ρL*QR2*ρR, zero(Q))))

        @test all(isapprox.(CMPSKit.localgradientRs(ψ̂[1]'*ψ̂[2], Q, Rs, ρL, ρR),
                            (ρL*R2*ρR, zero(Q))))
        @test all(isapprox.(CMPSKit.localgradient∂Rs(ψ̂[1]'*ψ̂[2], Q, Rs, ρL, ρR),
                            (zero(Q), zero(Q))))

        @test all(isapprox.(CMPSKit.localgradientRs((ψ̂[1]')^2*ψ̂[1]^2, Q, Rs, ρL, ρR),
                            (R1'*ρL*R1*R1*ρR + ρL*R1*R1*ρR*R1', zero(Q))))
    end
end

@testset "PeriodicCMPS: global gradients with bond dimension D = $D" for D in Dlist
    α = fit(x->-1 + 0.8*sin(x), FourierSeries; Kmax = 1)
    β = rand()
    γ = rand()
    H = ∫(∂ψ̂'*∂ψ̂ + α*ψ̂'*ψ̂ + β*(ψ̂*ψ̂ + ψ̂'*ψ̂') + γ*(ψ̂')^2*ψ̂^2, (-Inf,+Inf))
    T = ComplexF64
    Q = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D)) for j=1:5])
    R = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D))/D for j=1:3])
    Ψ = InfiniteCMPS(Q, R)
    ρL, ρR = environments!(Ψ; Kmax = 20)

    HL, = leftenv(H, (Ψ, ρL, ρR); Kmax = 20)
    HR, = rightenv(H, (Ψ, ρL, ρR); Kmax = 20)

    gradQ, gradRs = gradient(H, (Ψ, ρL, ρR), HL, HR)

    QL = Ψ.Q
    RL = Ψ.Rs[1]
    QRL = QL*RL-RL*QL
    𝒟R = QRL + ∂(RL)
    gradR = gradRs[1]

    @test gradQ ≈ ρL*𝒟R*ρR*RL' - RL'*ρL*𝒟R*ρR + HL*ρR + ρL*HR

    @test gradR ≈ -differentiate(ρL*𝒟R*ρR) + QL'*ρL*𝒟R*ρR - ρL*𝒟R*ρR*QL' +
                    α*ρL*RL*ρR + β*RL'*ρL*ρR + β*ρL*ρR*RL' + γ*RL'*ρL*RL*RL*ρR + γ*ρL*RL*RL*ρR*RL' + HL*RL*ρR + ρL*RL*HR
end

@testset "PeriodicCMS: test ground state algorithm" begin
    D = 2
    T = ComplexF64
    α = fit(x->-1 + 0.8*sin(x), FourierSeries; Kmax = 1)
    β = 0.
    γ = 1.
    H = ∫(∂ψ̂'*∂ψ̂ + α*ψ̂'*ψ̂ + β*(ψ̂*ψ̂ + ψ̂'*ψ̂') + γ*(ψ̂')^2*ψ̂^2, (-Inf,+Inf))
    Kmax = 10

    eigalg = Arnoldi(; krylovdim = D^2*(2*Kmax+1), tol = 1e-10)
    linalg = GMRES(; krylovdim = D^2*(2*Kmax+1), tol = 1e-10, maxiter = 1)
    optalg = ConjugateGradient(; gradtol = 1e-7, verbosity = 2)

    gradtol = 1e-7
    optalg = LBFGS(30; verbosity = 2, gradtol = gradtol)
    eigalg = Arnoldi(; krylovdim = 64, tol = 1e-10)
    linalg = GMRES(; krylovdim = 64, tol = 1e-10)
    for k = 1:3
        A = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D)) for j=1:5])
        A = (A - A')/2
        R = FourierSeries([exp(-4*(j>>1))*randn(T, (D,D))/D for j=1:3])
        Q = A - 1/2 * R'*R
        Ψ = InfiniteCMPS(Q, R)

        ΨL, ρL, ρR, E, e, normgrad, numfg, history =
            groundstate(H, Ψ;
                        optalg = optalg, eigalg = eigalg, linalg = linalg, Kmax = Kmax)

        @test E ≈ -0.237009267723921
        @test ∫(expval(ψ̂'*ψ̂, ΨL, one(ρR), ρR), (0,1)) ≈ 0.5616439424330847 atol=gradtol
    end
end


text:
#Global Setting
use_gpu = false
use_darrays = false #Only relevant if use_gpu = false

hard_clustering = false
ignore_local = false

global_weight = 1.0
local_weight= 1.0

split_delays = false

# #Model Parameters
# iterations = 10
#
# total_dim = 3
# local_dim = 3
#
# α = 4.0
# γ = 2.0
# global_weight = 1.0
# local_weight= 1.0
#
#
# global_hyper_params = niw_hyperparams(5.0,
#     ones(3)*2.5,
#     10.0,
#     [[0.8662817 0.78323282 0.41225376];[0.78323282 0.74170384 0.50340258];[0.41225376 0.50340258 0.79185577]])
#
# local_hyper_params = niw_hyperparams(1.0,
#     [217.0,510.0],
#     10.0,
#     Matrix{Float64}(I, 2, 2)*0.5)


text:
# ##############################  Gradients  ###################################

# function grad_potential_energy(r::Atom, p::GaN, Type1::Symbol, Type2::Symbol)

#     if (Type1 == :Ga) && (Type2 == :Ga) # Ga-Ga interaction
#         return grad_potential_energy(r.Position, p.lj_Ga_Ga)
#     elseif (Type1 == :N) && (Type2 == :N) # N-N interaction
#         return grad_potential_energy(r.Position, p.lj_N_N)
#     else
#         return grad_potential_energy(r.Position, p.bm_Ga_N)
#     end
# end

# function grad_potential_energy(r::Vector{Atom}, p::GaN)
#     n = length(r)
#     # Initialize
#     rtemp = r[1].Position - r[2].Position
#     pe = [0. * grad_potential_energy(rtemp, p, :Ga, :Ga), 0. * grad_potential_energy(rtemp, p, :N, :N), 0. * grad_potential_energy(rtemp, p, :Ga, :N)]
#     for i = 1:(n-1)
#         ri = r[i].Position
#         ri_Type = r[i].Type
#         for j = (i+1):n
#             rj = r[j].Position
#             rj_Type = r[j].Type
#             rtemp =  ri - rj
#             if ri_Type == r[j].Type
#                 if r[i].Type == :Ga
#                     pe[1] += grad_potential_energy(rtemp, p, ri_Type, rj_Type)
#                 else
#                     pe[2] += grad_potential_energy(rtemp, p, ri_Type, rj_Type)
#                 end
#             else
#                 pe[3] += grad_potential_energy(rtemp, p, ri_Type, rj_Type)
#             end
#         end
#     end
#     return pe
# end

# function grad_force(r::Atom, p::GaN, Type1::Symbol, Type2::Symbol)

#     if (Type1 == :Ga) && (Type2 == :Ga) # Ga-Ga interaction
#         return grad_force(r.Position, p.lj_Ga_Ga)
#     elseif (Type1 == :N) && (Type2 == :N) # N-N interaction
#         return grad_force(r.Position, p.lj_N_N)
#     else
#         return grad_force(r.Position, p.bm_Ga_N)
#     end
# end

# function grad_force(r::Vector{Atom}, p::GaN)
#     # Initialize
#     rtemp = r[1].Position - r[2].Position
#     f = [ [0. * grad_force(rtemp, p, :Ga, :Ga), 0. * grad_force(rtemp, p, :N, :N), 0. * grad_force(rtemp, p, :Ga, :N)] for j = 1:n]
#     for i = 1:(n-1)
#         ri = r[i].Position
#         for j = (i+1):n
#             rj = r[j].Position
#             rtemp = ri - rj
#             if r[i].Type == r[j].Type
#                 if r[i].Type == :Ga
#                     f[i][1] += grad_force(rtemp, p, r[i].Type, r[j].Type)
#                     f[j][1] -= f[i][1]
#                 else
#                     f[i][2] += grad_force(rtemp, p, r[i].Type, r[j].Type)
#                     f[j][2] -= f[i][2]
#                 end
#             else
#                 f[i][3] += grad_force(rtemp, p, r[i].Type, r[j].Type)
#                 f[j][3] -= f[i][3]
#             end
#         end
#     end
#     return f
# end

# function grad_virial(r::Atom, p::GaN, Type1::Symbol, Type2::Symbol)
#     df = grad_force(r, p, Type1, Type2)
#     if Type1 != Type2
#         dfdA = df[:dfdA]
#         dfdrho = df[:dfdρ]
#         return (dvdA = dfdA[1]*r.Position[1] + dfdA[2]*r.Position[2] + dfdA[3]*r.Position[3],
#             dvdρ = dfdrho[1]*r.Position[1] + dfdrho[2]*r.Position[2] + dfdrho[3]*r.Position[3])
#     else
#         dfde = df[:dfdϵ]
#         dfdsig = df[:dfdσ]
#     return (dvdϵ = dfde[1]*r.Position[1] + dfde[2]*r.Position[2] + dfde[3]*r.Position[3],
#             dvdσ = dfdsig[1]*r.Position[1] + dfdsig[2]*r.Position[2] + dfdsig[3]*r.Position[3])
#     end
# end

# function grad_virial(r::Vector{Atom}, p::GaN)
#     n = length(r)
#     # Initialize
#     rtemp = r[1].Position - r[2].Position
#     pe = [0. * grad_virial(rtemp, p, :Ga, :Ga), 0. * grad_virial(rtemp, p, :N, :N), 0. * grad_virial(rtemp, p, :Ga, :N)]
#     for i = 1:(n-1)
#         for j = (i+1):n
#             rtemp = r[i] - r[j]
#             if r[i].Type == r[j].Type
#                 if r[i].Type == :Ga
#                     pe[1] += grad_virial(rtemp, p, r[i].Type, r[j].Type)
#                 else
#                     pe[2] += grad_virial(rtemp, p, r[i].Type, r[j].Type)
#                 end
#             else
#                 pe[3] += grad_virial(rtemp, p, r[i].Type, r[j].Type)
#             end
#         end
#     end
#     return pe
# end


text:
@sumProductRule(:node_type     => Bernoulli,
                :outbound_type => Message{Bernoulli},
                :inbound_types => (Nothing, Message{PointMass}),
                :name          => SPBernoulliOutNP)

@sumProductRule(:node_type     => Bernoulli,
                :outbound_type => Message{Beta},
                :inbound_types => (Message{PointMass}, Nothing),
                :name          => SPBernoulliIn1PN)

@sumProductRule(:node_type     => Bernoulli,
                :outbound_type => Message{Bernoulli},
                :inbound_types => (Nothing, Message{Beta}),
                :name          => SPBernoulliOutNB)

@naiveVariationalRule(:node_type     => Bernoulli,
                      :outbound_type => Message{Bernoulli},
                      :inbound_types => (Nothing, ProbabilityDistribution),
                      :name          => VBBernoulliOut)

@naiveVariationalRule(:node_type     => Bernoulli,
                      :outbound_type => Message{Beta},
                      :inbound_types => (ProbabilityDistribution, Nothing),
                      :name          => VBBernoulliIn1)


text:
##Ok let's go and write an iterator for Kelleher's algorithm
struct part
  n
end

import Base
#import IterativeSolvers
import SpecialFunctions

#initializer
function Base.iterate(iter::part)
  a=zeros(Int64,iter.n)
  k=2
  y=iter.n-1
  if k!=1
    x=a[k-1]+1
    k-=1
    while 2*x <=y
      a[k]=x
      y-=x
      k+=1
    end
    l=k+1
    if x<=y
      a[k]=x
      a[l]=y
      #put!(c,a[1:k+1])
      x+=1
      y-=1
      return a[1:k+1],(a,k,l,x,y)
    end
    a[k]=x+y
    y=x+y-1
    ## header again there is probably an easier way to code this
    p= a[1:k]
    if k!=1
      x=a[k-1]+1
      k-=1
      while 2*x <=y
        a[k]=x
        y-=x
        k+=1
      end
      l=k+1
    else
      k=-1 # sentinel value see below
    end
    return p,(a,k,l,x,y)
  end
end

function Base.iterate(iter::part,state)
  a,k,l,x,y=state
  if k==-1; return nothing; end #minus 1 is a sentinel value
  if x<=y
    a[k]=x
    a[l]=y
    #put!(c,a[1:k+1])
    x+=1
    y-=1
    return a[1:k+1],(a,k,l,x,y)
  else
    a[k]=x+y
    y=x+y-1
    p=a[1:k]# copy next output
    if k!=1
      x=a[k-1]+1
      k-=1
      while 2*x <=y
        a[k]=x
        y-=x
        k+=1
      end
      l=k+1
    else
      k=-1 #set sentinel value for termination
    end
    return p,(a,k,l,x,y)
  end
end

# #small performance test
# #Watch this PyHoltz
# zeit=time()
# for (num,p) in enumerate(part(15))
#   println(num,p)
# end
# zeit=time()-zeit
# print("iterator:",zeit)

function part2mult(p)
  z=sum(p)
  mu=zeros(Int64,z)
  if p !=[0]
    for i in p
      mu[i]+=1
    end
  end
  return mu
end

# #next tiny test
# for (num,p) in enumerate(part(5))
#   println(num,p,part2mult(p))
# end

function multinomcoeff(mu)
  return SpecialFunctions.factorial(float(sum(mu)))/prod(SpecialFunctions.factorial.(float.(mu)))
end

function weigh(mu)
  weight=0
  for (g,mu_g) in enumerate(mu)
    weight+=g*mu_g #TODO: check type
  end
  return weight
end

function generate_MN(N)
  #MN=[Dict{Tuple{Int64,Int64},Array{Array{Int64,1},1}}()]
  MN=[Dict{Array{Int64,1},Array{Array{Int64,1},1}}()]
  #MN=[Dict()]
  for k=1:N
    #MU=Dict{Tuple{Int64,Int64},Array{Array{Int64,1},1}}()
    MU=Dict{Array{Int64,1},Array{Array{Int64,1},1}}()
    #MU=Dict()
    for n=1:k
      key=[0,n]
      mu=[0]
      mu=part2mult(mu)
      #mu=Array{Int64,1}[]
      if key in keys(MU)
        push!(MU[key],mu)
      else
        MU[key]=[mu]
      end
    end

    for m = 1:k
      for mu in part(m)
        if mu == [k]
          continue
        end
        mu=part2mult(mu)
        for n=0:k-m
          key=[sum(mu),n]
          if key in keys(MU)
            push!(MU[key],mu)
          else
            MU[key]=[mu]
          end
        end
      end
    end
    push!(MN,MU)
  end
  return MN
end

function tuple2idx(tpl,ord)
  m,n=tpl
  return sum(ord-m+1:ord)+n+m
end

function mult2str(mu)
  txt=""
  for mu_g in mu
    txt*=string(mu_g)*" "
  end
  if length(txt)!=0
    txt=txt[1:end-1]
    txt*="\n"
  end
 return txt
end

function generate_multi_indices_at_order(k;to_disk=false,compressed=false)
  if to_disk
    pack="../src/NLEVP/compressed_perturbation_data/" #This file location is relative to the deps directory as this is the one where build pkg is run.
    dir="$k/"
  else
    Mu=[ Array{Int16,1}[] for i in 1:tuple2idx((k,0),k) ] #TODO specify type
  end
  if to_disk && !ispath(pack*dir)
    mkpath(pack*dir)
  elseif to_disk && ispath(pack*dir)
    return nothing
  end
  for n=1:k
    key=(0,n)
    p=[0]
    mu=part2mult(p)
    out=compressed ? p :  mu
    println("n and k : $n and $k")
    println("key: $key")
    if to_disk
      fname="$(key[1])_$(key[2])"
      println("fname: $fname")
      open(pack*dir*fname,"w") do file
        write(file,mult2str(out))
      end
    else
      idx=tuple2idx(key,k)
      push!(Mu[idx],out)
    end
  end

  for m = 1:k
    for p in part(m)
      if p == [k]
        continue
      end
      mu=part2mult(p)
      out=compressed ? p :  mu
      for n=0:k-m
        key=(sum(mu),n)

        if to_disk
          fname="$(key[1])_$(key[2])"
          open(pack*dir*fname,"a") do file
            write(file,mult2str(out))
          end
        else
          idx=tuple2idx(key,k)
          push!(Mu[idx],out)
        end

      end
    end
  end

  if !to_disk
    return Mu
  end
end

function efficient_MN(N)
  MN=Array{Array{Array{Array{Int16,1},1},1}}(undef,N)
  #MN=[]
  for k =1:N
  #push!(MN,generate_multi_indices_at_order(k))
  MN[k]=generate_multi_indices_at_order(k)
  end
  return MN
end

function disk_MN(N)
  for k=1:N
    generate_multi_indices_at_order(k,to_disk=true,compressed=true)
  end
end
##
# zeit=time()
# MN=generate_MN(50);
# zeit=time()-zeit
# println("zeit: ", zeit)
##
# function perturb(L,N,v0,v0Adj)
#   #normalize
#   v0/=sqrt(v0'*v0)
#   v0Adj/=v0Adj'*L(1,0)*v0
#   λ=Array{ComplexF64}(undef,N+1)
#   v=Array{Array{Complex{Float64}},1}(undef,N+1)
#   v[1]=v0
#   dim=size(v0)[1]
#   L00=L(0,0)
#   sparse = SparseArrays.issparse(L00)
#   L00=SparseArrays.lu(L(0,0),check=false) #L(0,0) is singular but in most of the cases we are lucky and a LU factorization exists
#   if sparse && L00.status!=0 #lu failed... TODO: implement LU check for dense arrays
#      L00=SparseArrays.qr(L(0,0)) #we try a QR factorization
#   end
#   #TODO consider using qr for all cases
#
#
#   for k=1:N
#     r=zeros(ComplexF64,dim)
#     for n = 1:k
#       r+=L(0,n)*v[k-n+1] #+1 to start indexing with zero
#     end
#     for m =1:k
#       for mu in part(m)
#         if mu ==[k]
#           continue
#         end
#         mu=part2mult(mu)
#         for n=0:k-m
#           coeff=1
#           for (g, mu_g) in enumerate(mu)
#             coeff*=λ[g+1]^mu_g #+1 because indexing starts with zero
#           end
#           r+=L(sum(mu),n)*v[k-n-m+1]*multinomcoeff(mu)*coeff
#         end
#       end
#     end
#     λ[k+1]=-v0Adj'*r/(v0Adj'*L(1,0)*v0)
#     #v[k+1]= IterativeSolvers.gmres(L(0,0),-(r+λ[k+1]*L(1,0)*v0)) This is crap!!!
#     #v[k+1]=L(0,0)\-(r+λ[k+1]*L(1,0)*v0) #This works but LU factorization should be done outside of the loop
#     v[k+1]=L00\-(r+λ[k+1]*L(1,0)*v0) #yeah !
#     v[k+1]-=(v0'*v[k+1])*v0 #if v[k+1] is orthogonal to v0 it features minimum norm
#     #
#     #println("status: $(L00.status)")
#     #res=L(0,0)*v[k+1]+(r+λ[k+1]*L(1,0)*v0)
#     #res=LinearAlgebra.norm(res)
#   end
#   return λ,v
# end

function perturb(L,N,v0,v0Adj)
  #normalize
  v0/=sqrt(v0'*v0)
  v0Adj/=v0Adj'*L(1,0)*v0
  λ=Array{ComplexF64}(undef,N+1)
  v=Array{Array{Complex{Float64}},1}(undef,N+1)
  v[1]=v0
  dim=size(v0)[1]
  L00=L(0,0)
  sparse = SparseArrays.issparse(L00)
  L00=SparseArrays.lu(L(0,0),check=false) #L(0,0) is singular but in most of the cases we are lucky and a LU factorization exists
  if sparse && L00.status!=0 #lu failed... TODO: implement LU check for dense arrays
     L00=SparseArrays.qr(L(0,0)) #we try a QR factorization
  end
  #TODO consider using qr for all cases

  for k=1:N
    r=zeros(ComplexF64,dim)
    for n = 1:k
      r+=L(0,n)*v[k-n+1] #+1 to start indexing with zero
    end
    for m =1:k
      for mu in part(m)
        if mu ==[k]
          continue
        end
        mu=part2mult(mu)
        for n=0:k-m
          coeff=1
          for (g, mu_g) in enumerate(mu)
            coeff*=λ[g+1]^mu_g #+1 because indexing starts with zero
          end
          r+=L(sum(mu),n)*v[k-n-m+1]*multinomcoeff(mu)*coeff
        end
      end
    end
    λ[k+1]=-v0Adj'*r/(v0Adj'*L(1,0)*v0)
    #v[k+1]= IterativeSolvers.gmres(L(0,0),-(r+λ[k+1]*L(1,0)*v0)) This is crap!!!
    #v[k+1]=L(0,0)\-(r+λ[k+1]*L(1,0)*v0) #This works but LU factorization should be done outside of the loop
    v[k+1]=L00\-(r+λ[k+1]*L(1,0)*v0) #yeah !
    v[k+1]-=(v0'*v[k+1])*v0 #if v[k+1] is orthogonal to v0 it features minimum norm
    #
    #println("status: $(L00.status)")
    #res=L(0,0)*v[k+1]+(r+λ[k+1]*L(1,0)*v0)
    #res=LinearAlgebra.norm(res)
  end
  return λ,v
end

#TODO: consider stronger compression
function perturb_disk(L,N,v0,v0Adj)
  #normalize
  zeit=time()
  v0/=sqrt(v0'*v0)
  v0Adj/=v0Adj'*L(1,0)*v0
  λ=Array{ComplexF64}(undef,N+1)
  v=Array{Array{Complex{Float64}},1}(undef,N+1)
  v[1]=v0
  dim=size(v0)[1]
  L00=L(0,0)
  sparse = SparseArrays.issparse(L00)
  L00=SparseArrays.lu(L(0,0),check=false) #L(0,0) is singular but in most of the cases were are lucky and a LU factorization exists
  if sparse && L00.status!=0 #lu failed... TODO: implement LU check for dense arrays
     L00=SparseArrays.qr(L(0,0)) #we try a QR factorization
  end
  #TODO consider using qr for all cases
  pack=(@__DIR__)*"/compressed_perturbation_data/"
  for k=1:N
    dir="$k/"
    r=zeros(ComplexF64,dim)
    for m=0:k
      for n=0:k-m
        if m==n==0 || k==m==1
          continue
        end
        fname="$(m)_$(n)"
        w=zeros(ComplexF64,dim)
        open(pack*dir*fname,"r") do file
          for str in eachline(file)
            #str=readline(file)
            mu=part2mult(parse.(Int64,split(str)))

            coeff=1
            for (g,mu_g) in enumerate(mu)
              coeff*=λ[g+1]^mu_g #+1 because indexing starts with zero
            end
            w+=v[k-n-weigh(mu)+1]*multinomcoeff(mu)*coeff
            end
        end
        r+=L(m,n)*w
      end
    end
    #println("##########")
    #println("r: $r")
    #v[k+1]= IterativeSolvers.gmres(L(0,0),-(r+λ[k+1]*L(1,0)*v0)) This is crap!!!
    λ[k+1]=-v0Adj'*r/(v0Adj'*L(1,0)*v0)
    #println("λ: $(λ[k+1])")
    #v[k+1]=L(0,0)\-(r+λ[k+1]*L(1,0)*v0) #This works but LU factorization should be done outside of the loop
    #println("rhs: $(r+λ[k+1]*L(1,0)*v0)")
    v[k+1]=L00\-(r+λ[k+1]*L(1,0)*v0) #yeah !
    v[k+1]-=(v0'*v[k+1])*v0 #if v[k+1] is orthogonal to v0 it features minimum norm
    #v[k+1]./=exp(1im*angle(v[k+1][1]))
    #res=L(0,0)*v[k+1]+(r+λ[k+1]*L(1,0)*v0)
    #res=LinearAlgebra.norm(res)

    #normalization
    c=0.0+0.0im
    #println("order: $k, time: $(time()-zeit)")
    for l=1:k-1
        c-=.5*v[l+1]'*v[k-l+1]
        #c+=v[l+1]'*v[k-l+1]

    end
    v[k+1]+=c*v[1]
    #v[k+1]-=v[1]*(c/2.0+v[1]'*v[k+1])

  end
  return λ,v
end

# function perturb_fast(L,N,v0,v0Adj)
#   #normalize
#   v0/=sqrt(v0'*v0)
#   v0Adj/=v0Adj'*L(1,0)*v0
#   λ=Array{ComplexF64}(undef,N+1)
#   v=Array{Array{Complex{Float64}},1}(undef,N+1)
#   v[1]=v0
#   dim=size(v0)[1]
#   L00=L(0,0)
#   sparse = SparseArrays.issparse(L00)
#   L00=SparseArrays.lu(L(0,0),check=false) #L(0,0) is singular but in most of the cases were are lucky and a LU factorization exists
#   if sparse && L00.status!=0 #lu failed... TODO: implement LU check for dense arrays
#      L00=SparseArrays.qr(L(0,0)) #we try a QR factorization
#   end
#   #TODO consider using qr for all cases
#
#
#   for k=1:N
#     r=zeros(ComplexF64,dim)
#
#     for (m,n) in MN[k+1] #MN is ein dictionary mit schlüsseln (m,n)
#       w=zeros(ComplexF64,dim)
#       for mu in MN[k+1][(m,n)]
#         coeff=1
#         for (g,mu_g) in enumerate(mu)
#           coeff*=λ[g+1]^mu_g #+1 because indexing starts with zero
#           w+=v[k-n-weigh(mu)+1]*multinomcoeff(mu)*coeff
#         end
#       end
#       r+=L(m,n)*w
#     end
#     #v[k+1]= IterativeSolvers.gmres(L(0,0),-(r+λ[k+1]*L(1,0)*v0)) This is crap!!!
#
#     λ[k+1]=-v0Adj'*r/(v0Adj'*L(1,0)*v0)
#     #v[k+1]=L(0,0)\-(r+λ[k+1]*L(1,0)*v0) #This works but LU factorization should be done outside of the loop
#     v[k+1]=L00\-(r+λ[k+1]*L(1,0)*v0) #yeah !
#     res=L(0,0)*v[k+1]+(r+λ[k+1]*L(1,0)*v0)
#   end
#   return λ,v
# end

function perturb_norm(L,N,v0,v0Adj)
  #normalize
  zeit=time()
  Y=-L.terms[end].coeff #TODO: check for __aux__
  Y=convert(SparseArrays.SparseMatrixCSC{ComplexF64,Int},Y)
  v0/=sqrt(v0'*Y*v0)
  v0Adj=SparseArrays.lu(Y)\v0Adj #TODO: this step is highly redundant
  v0Adj/=v0Adj'*Y*L(1,0)*v0
  λ=Array{ComplexF64}(undef,N+1)
  v=Array{Array{Complex{Float64}},1}(undef,N+1)
  v[1]=v0
  dim=size(v0)[1]
  L00=L(0,0)
  sparse = SparseArrays.issparse(L00)
  L00=SparseArrays.lu(L(0,0),check=false) #L(0,0) is singular but in most of the cases were are lucky and a LU factorization exists
  if sparse && L00.status!=0 #lu failed... TODO: implement LU check for dense arrays
     L00=SparseArrays.qr(L(0,0)) #we try a QR factorization
  end
  #TODO consider using qr for all cases
  pack=(@__DIR__)*"/compressed_perturbation_data/"
  for k=1:N
    dir="$k/"
    r=zeros(ComplexF64,dim)
    for m=0:k
      for n=0:k-m
        if m==n==0 || k==m==1
          continue
        end
        fname="$(m)_$(n)"
        w=zeros(ComplexF64,dim)
        open(pack*dir*fname,"r") do file
          for str in eachline(file)
            #str=readline(file)
            mu=part2mult(parse.(Int64,split(str)))

            coeff=1
            for (g,mu_g) in enumerate(mu)
              coeff*=λ[g+1]^mu_g #+1 because indexing starts with zero
            end
            w+=v[k-n-weigh(mu)+1]*multinomcoeff(mu)*coeff
            end
        end
        r+=L(m,n)*w
      end
    end
    #println("##########")
    #println("r: $r")
    #v[k+1]= IterativeSolvers.gmres(L(0,0),-(r+λ[k+1]*L(1,0)*v0)) This is crap!!!
    λ[k+1]=-v0Adj'*Y*r/(v0Adj'*Y*L(1,0)*v0)
    #println("λ: $(λ[k+1])")
    #v[k+1]=L(0,0)\-(r+λ[k+1]*L(1,0)*v0) #This works but LU factorization should be done outside of the loop
    #println("rhs: $(r+λ[k+1]*L(1,0)*v0)")
    v[k+1]=L00\-(r+λ[k+1]*L(1,0)*v0) #yeah !
    v[k+1]-=(v0'*Y*v[k+1])*v0 #if v[k+1] is orthogonal to v0 it features minimum norm
    #v[k+1]./=exp(1im*angle(v[k+1][1]))
    #res=L(0,0)*v[k+1]+(r+λ[k+1]*L(1,0)*v0)
    #res=LinearAlgebra.norm(res)

    #normalization
    c=0.0+0.0im
    println("order: $k, time: $(time()-zeit)")
    for l=1:k-1
        c-=.5*v[l+1]'*Y*v[k-l+1]
        #c+=v[l+1]'*v[k-l+1]

    end
    v[k+1]+=c*v[1]
    #v[k+1]-=v[1]*(c/2.0+v[1]'*v[k+1])

  end
  return λ,v
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

abstract type AbstractRule end

struct Rule <: AbstractRule
    label::String
    when::String
    then::String
end

function Rule(rule::Proto.Rule)
    label = rule.label
    when = rule.when
    then = rule.then

    return Rule(label, when, then)
end

function Rule(rules::Vector{Proto.Rule})
    return Rule.(rules)
end

"""
    Remote

Wrapper type that encapsulates a Rule and transaction.
"""
struct RemoteRule{D <: AbstractRule, T <: AbstractCoreTransaction}
    rule::D
    transaction::T
end

"""
    as_remote(x, t)

Create a `Remote`(@ref) object for a concept `x` with a
transaction `t`.
"""
function as_remote(x::D, t::T) where {
    D <: AbstractRule, T <: AbstractCoreTransaction
}
    return RemoteRule{D, T}(x, t)
end

function Base.show(io::IO, rule::Rule)
    println(io, "Rule:")
    println(io, "label: $(rule.label)")
    println(io, "when: $(rule.when)")
    println(io, "then: $(rule.then)")
end

is_remote(rule::Rule) = false
is_remote(rule::RemoteRule) = true

function set_label(remote::RemoteRule, new_label::AbstractString)
    dbResult = execute(remote.rule.transaction, RuleRequestBuilder.rule_set_label_req(remote.label, new_label))
    remote.rule.label = new_label
    return remote
end

"""
    delete(remote::RemoteRule)
To delete a type in the database pack the type with the transaction to a RemoteRule.
"""
function delete(remote::RemoteRule)
    execute(remote.transaction, RuleRequestBuilder.rule_delete_req(remote.rule.label))
end

function is_deleted(remote::RemoteRule)
    dbResult = execute(remote.transaction, LogicManagerRequestBuilder.get_rule_req(remote.rule.label))
    return !hasproperty(dbResult[1].logic_manager_res.get_rule_res,:res)
end


text:
module PowerFlowAnalysis

using JuMP
using Ipopt

using PowerModels

# package code goes here
include("core/common.jl")

include("core/import.jl")
include("core/pfsys.jl")
include("core/lurie.jl")

include("solvability/brouwer.jl")
include("solvability/luriecert.jl")

end # module

max_stars_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_issues_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_forks_repo_name: hunghtd/PowerFlowAnalysis-tutorial


text:
# This file was generated, do not modify it. # hide
abs(x) < eps()

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
function num_vecjac!(
    du,
    f,
    x,
    v,
    cache1 = similar(v),
    cache2 = similar(v);
    compute_f0 = true,
)
    if DiffEqBase.numargs(f) != 2
        du .= num_jacvec(f, x, v)
        return du
    end
    compute_f0 && (f(cache1, x))
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    vv = reshape(v, size(x))
    for i = 1:length(x)
        x[i] += ϵ
        f(cache2, x)
        x[i] -= ϵ
        du[i] = (((cache2 .- cache1) ./ ϵ)'*vv)[1]
    end
    return du
end

function num_vecjac(f, x, v, f0 = nothing)
    vv = reshape(v, axes(x))
    f0 === nothing ? _f0 = f(x) : _f0 = f0
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    du = similar(x)
    for i = 1:length(x)
        x[i] += ϵ
        f0 = f(x)
        x[i] -= ϵ
        du[i] = (((f0 .- _f0) ./ ϵ)'*vv)[1]
    end
    return vec(du)
end


text:
using sb, LinearAlgebra, Statistics

# observable storage structure
struct storages
    time :: Vector{Float64}
    energy :: Vector{Float64}
    function storages()
        new(Vector{Float64}(),Vector{Float64}())
    end
end

# storage callback
function store!(storage, t, ops, sbm)
    append!(storage.time,t)
    append!(storage.energy,sb.ehenergy(ops,sbm))
end

function integrate()
    storage = storages()
    # Create a model
    sbm = SBModel(OhmicSBParams(No = 100, α = 0.1))
    # Create containers starting from a localized state
    ops = EhrenfestOps([[1.0 0.0]; [0.0 0.0]],sbm)
    # run dynamics
    RunEhrenfest!(sbm,ops,2000,0.001,store!, storage)
    return storage.time,storage.energy
end

ts, energies = integrate()

@test mean(energies.-energies[1]) < 0.003


text:
# ---
# title: 383. Ransom Note
# id: problem383
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/ransom-note/description/>
# hidden: true
# ---
#
# Given an arbitrary ransom note string and another string containing letters
# from all the magazines, write a function that will return true if the ransom
# note can be constructed from the magazines ; otherwise, it will return false.
#
# Each letter in the magazine string can only be used once in your ransom note.
#
#
#
# **Example 1:**
#
#
#
#     Input: ransomNote = "a", magazine = "b"
#     Output: false
#
#
# **Example 2:**
#
#
#
#     Input: ransomNote = "aa", magazine = "ab"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: ransomNote = "aa", magazine = "aab"
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * You may assume that both strings contain only lowercase letters.
#
#
## @lc code=start
using LeetCode
using DataStructures

function can_construct(ransom_note::String, magazine::String)::Bool
    words = DefaultDict{Char,Int}(0)
    for i in magazine
        words[i] += 1
    end
    for i in ransom_note
        words[i] == 0 && return false
        words[i] -= 1
    end
    return true
end
## @lc code=end


text:
module Aurora

using Reexport

using LaTeXStrings
using LinearAlgebra
using MultivariateStats
using NearestNeighbors
using RecipesBase

using Statistics
@reexport using StatsBase

abstract type AbstractAurora end
abstract type AbstractFittedAurora end
StatsBase.predict(fitted_aurora::AbstractFittedAurora) = fitted_aurora.μs

include("replicated_sample.jl")
include("auroral.jl")
include("auroraknn.jl")

export ReplicatedSample,
       Auroral,
       CoeyCunningham,
       AuroraKNN

end


text:
# on Julia prompt:
# julia> STDIN
# TTY(open, 0 bytes waiting)
# julia> STDOUT
# TTY(open, 0 bytes waiting)

read(STDIN, Char)
write(STDOUT, "Julia") #> Julia5
readbytes(STDIN,3)
readline(STDIN)
# input:  # Julia
# output: # "Julia\r\n" (Windows) - "Julia\n" (Linux)

names(IOStream)
# 4-element Array{Symbol,1}:
#  :handle
#  :ios
#  :name
#  :mark

IOStream.types
# (Ptr{None},Array{Uint8,1},String,Int64)

# stream can be STDIN or any other input stream:
for line in eachline(stream)
    print("Found $line")
end

while !eof(stream)
    x = read(stream, Char)
    println("Found: $x")
end

# files:
# reading:
fname = "example.dat"
f1 = open(fname)
# IOStream(<file example.dat>)

data = readlines(f1)
# 3-element Array{Union(ASCIIString,UTF8String),1}:
#  "this is line 1.\r\n"
#  "this is line 2.\r\n"
#  "this is line 3."

for line in data
	println(line)
	# process line
end
close(f1)

open(fname) do file
    process(file)
end

function process()
	# ...
end

open(fname) do file
    for line in eachline(file)
        print(line) # process line
    end
end

# writing:
fname = "example2.dat"
f2 = open(fname, "w")
write(f2, "I write myself to a file\n") #> 24
println(f2, "even with println!")
close(f2)

# looping over all files in a folder:
for file in readdir()
	# process file
end

max_stars_repo_path: Module 1/Chapter08/Old/io.jl

max_issues_repo_path: Module 1/Chapter08/Old/io.jl

max_forks_repo_path: Module 1/Chapter08/Old/io.jl


text:
export train_local!

function train_local!{T<:Number}(net::HopfieldNetwork, patterns::Array{T,2})
	nn,np = size(patterns)

	@assert false "unimplemented"
	# calc synaptic weight (Local iterative learning method)
	train(net,patterns)
	w = zeros(T,size(net.weight))
	for mu=1:np
		for p=1:nn
			for q=1:nn
				w[p,q] = patterns[p,mu] * conj(patterns[q,mu]) / nn
			end
		end
		net.weight += d*w
	end
end


text:
# Write functions as pure math
cartesian(r, θ, ϕ) = r*cos(θ)*sin(ϕ), r*sin(θ)*sin(ϕ), r*cos(ϕ)

f(x) = 2x^2+5x-10
∇f(x) = 4x+5
∇²f(x) = 4

# Set notation are defined operators
10 ∈ [12,13,14,15] # returns false

max_stars_repo_name: JamieMair/intro-to-julia-talk

max_issues_repo_name: JamieMair/intro-to-julia-talk

max_forks_repo_name: JamieMair/intro-to-julia-talk


text:
# ------------------------------------------------------------------------------------------
# # Julia for Data Science - Data
# Prepared by [@nassarhuda](https://github.com/nassarhuda)! 😃
#
# In the next few notebooks, we will discuss why *Julia* is the tool you want to use for
# your data science applications.
#
# We will cover the following:
# 1. Reading and writing files
# 1. DataFrames
# 1. RDatasets
# 1. FileIO
# 1. File types
#
#
# ### Data: Build a strong relationship with your data.
# Every data science task has one main ingredient, the _data_! Most likely, you want to use
# your data to learn something new. But before the _new_ part, what about the data you
# already have? Let's make sure you can **read** it, **store** it, and **understand** it
# before you start using it.
#
# Julia makes this step really easy with data structures and packages to process the data,
# as well as existing functions that are readily usable on your data.
#
# The goal of this first part is get you acquainted with some Julia's tools to manage your
# data.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# **Reading and writing to files is really easy in Julia.**

#
# To see this, let's download a csv file from github that we can work with.
#
# Note: `download` depends on external tools such as curl, wget or fetch. So you must have
# one of these.
# ------------------------------------------------------------------------------------------

P = download("https://raw.githubusercontent.com/nassarhuda/easy_data/master/programming_languages.csv","programminglanguages.csv")

# ------------------------------------------------------------------------------------------
# We can use shell commands like `ls` in Julia by preceding them with a semicolon.
# ------------------------------------------------------------------------------------------

;ls

# ------------------------------------------------------------------------------------------
# And there's the *.csv file we downloaded!
#
# By default, `readdlm` will fill an array with the data stored in the input .csv file. If
# we set the keyword argument `header` to `true`, we'll get a second output array for the
# headers.
# ------------------------------------------------------------------------------------------

using DelimitedFiles
P,H = readdlm("programminglanguages.csv",header=true)

P

# ------------------------------------------------------------------------------------------
# You can use different delimiters with the function `readdlm`.
#
# To write to files, we can use `writedlm`.

#
# Let's write this same data to a file with a different delimiter.
# ------------------------------------------------------------------------------------------

writedlm("programming_languages_data.txt", P, '-')

# ------------------------------------------------------------------------------------------
# We can now check that this worked using a shell command to glance at the file,
# ------------------------------------------------------------------------------------------

;head -10 programming_languages_data.txt

# ------------------------------------------------------------------------------------------
# and also check that we can use `readdlm` to read our new text file correctly.
# ------------------------------------------------------------------------------------------

P_new_delim = readdlm("programming_languages_data.txt", '-');
P == P_new_delim

# ------------------------------------------------------------------------------------------
# ### DataFrames!
# *Shout out to R fans!*
# One other way to play around with data in Julia is to use a DataFrame.
#
# This requires loading the `DataFrames` package.
#
# Run this command to install all the packages used in the "Julia for Data Science" project
# -- (those packages are listed in this file: [`Project.toml`](/edit/introductory-
# tutorials/broader-topics-and-ecosystem/intro-to-julia-for-data-science/Project.toml)):
# ------------------------------------------------------------------------------------------

] instantiate

using DataFrames
df = DataFrame(year = P[:,1], language = P[:,2])

# ------------------------------------------------------------------------------------------
# You can access columns by header name, or column index.
#
# In this case, `df[1]` is equivalent to `df.year` or `df[!, :year]`.
#
# Note that if we want to index columns by header name, we precede the header name with a
# colon. In Julia, this means that the header names are treated as *symbols*.
# ------------------------------------------------------------------------------------------

df.year

# ------------------------------------------------------------------------------------------
# **`DataFrames` provides some handy features when dealing with data**
#
# First, it uses julia's "missing" type.
# ------------------------------------------------------------------------------------------

a = missing
typeof(a)

# ------------------------------------------------------------------------------------------
# Let's see what happens when we try to add a "missing" type to a number
# ------------------------------------------------------------------------------------------

a + 1

# ------------------------------------------------------------------------------------------
# ### RDatasets
#
# We can use RDatasets to play around with pre-existing datasets
# ------------------------------------------------------------------------------------------

using RDatasets
iris = dataset("datasets", "iris")

# ------------------------------------------------------------------------------------------
# Note that data loaded with `dataset` is stored as a DataFrame. 😃
# ------------------------------------------------------------------------------------------

typeof(iris)

# ------------------------------------------------------------------------------------------
# `DataFrames` provides the `describe` can give you quick statistics about each column in
# your dataframe
# ------------------------------------------------------------------------------------------

describe(iris)

# ------------------------------------------------------------------------------------------
# You can create your own dataframe quickly as follows
# ------------------------------------------------------------------------------------------

foods = ["apple", "cucumber", "tomato", "banana"]
calories = [missing,47,22,105]
typeof(calories)

using Statistics
mean(calories)

# ------------------------------------------------------------------------------------------
# `missing` ruins everything! 😑
# ------------------------------------------------------------------------------------------

mean(skipmissing(calories))

# ------------------------------------------------------------------------------------------
# In fact, `describe' will drop these values too
# ------------------------------------------------------------------------------------------

describe(DataFrame(c=calories))

# ------------------------------------------------------------------------------------------
# Note that `typeof(calories)` is `Vector{Union{Missing, Int64},1}`
#
# If we want to replace all `missing` values a default value, we can do it like this:
# ------------------------------------------------------------------------------------------

newcalories = replace(calories, missing => 0)

# ------------------------------------------------------------------------------------------
# Now let's create a `DataFrame` that shows foods and their calories from two `DataArray`s!
# ------------------------------------------------------------------------------------------

dataframe_calories = DataFrame(item=foods,calories=calories)

# ------------------------------------------------------------------------------------------
# Let's generate a second `DataFrame` that shows foods and their prices.
# ------------------------------------------------------------------------------------------

prices = [0.85,1.6,0.8,0.6,]

dataframe_prices = DataFrame(item=foods,price=prices)

# ------------------------------------------------------------------------------------------
# We can also `join` these two dataframes together because they share a common column,
# `item`.
# ------------------------------------------------------------------------------------------

DF = join(dataframe_calories,dataframe_prices,on=:item)

# ------------------------------------------------------------------------------------------
# Note that we used the keyword argument `on` to say that we wanted to join these dataframes
# together based on the `item` column.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### FileIO
# ------------------------------------------------------------------------------------------

using FileIO
julialogo = download("https://avatars0.githubusercontent.com/u/743164?s=200&v=4","julialogo.png")

# ------------------------------------------------------------------------------------------
# Again, let's check that this download worked!
# ------------------------------------------------------------------------------------------

;ls

# ------------------------------------------------------------------------------------------
# Next, let's load the Julia logo, stored as a .png file
# ------------------------------------------------------------------------------------------

X1 = load("julialogo.png")

# ------------------------------------------------------------------------------------------
# We see below that Julia stores this logo as an array of colors.
# ------------------------------------------------------------------------------------------

@show typeof(X1);
@show size(X1);

# ------------------------------------------------------------------------------------------
# ### File types
# In Julia, many file types are supported so you do not have to transfer a file you got from
# another language to a text file before you read it.
#
# *Some packages that achieve this:*
# MAT CSV NPZ JLD FASTAIO
#
#
# Let's try using MAT to write a file that stores a matrix.
# ------------------------------------------------------------------------------------------

using MAT

A = rand(5,5)
matfile = matopen("densematrix.mat", "w")
write(matfile, "A", A)
close(matfile)

# ------------------------------------------------------------------------------------------
# Now try opening densematrix.mat with MATLAB!
# ------------------------------------------------------------------------------------------

newfile = matopen("densematrix.mat")
read(newfile,"A")

names(newfile)

close(newfile)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using ReservoirComputing
using ParameterizedFunctions
using OrdinaryDiffEq
using Plots

#lorenz system parameters
u0 = [1.0,0.0,0.0]
tspan = (0.0,1000.0)
p = [10.0,28.0,8/3]

#define lorenz system
function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

#solve system
prob = ODEProblem(lorenz, u0, tspan, p)
sol = solve(prob, ABM54(), dt=0.02)
v = sol.u
data = Matrix(hcat(v...))

#model parameters
shift = 300
approx_res_size = 300
radius = 1.2
activation = tanh
degree = 6
sigma = 0.1
train_len = 5000
predict_len = 1250
beta = 0.0
alpha = 1.0
nla_type = NLAT2()
extended_states = false

#get data
train = data[:, shift:shift+train_len-1]
test = data[:, shift+train_len:shift+train_len+predict_len-1]

#create echo state network
esn = ESN(approx_res_size,
    train,
    degree,
    radius,
    activation = activation,
    sigma = sigma,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

#training and prediction
W_out = ESNtrain(esn, beta)
output = ESNpredict(esn, predict_len, W_out)

#plots and images
comp = plot(transpose(output),layout=(3,1), label="predicted")
comp = plot!(transpose(test),layout=(3,1), label="actual")
savefig(comp, "lorenz_coord")
sf = plot(transpose(output)[:,1], transpose(output)[:,2], transpose(output)[:,3], label="predicted")
sf2 = plot(transpose(test)[:,1], transpose(test)[:,2], transpose(test)[:,3], label="actual")
final = plot(sf, sf2, layout=(1,2), label = ["predicted" "actual"])
savefig(final, "lorenz_attractor")


text:
function genrandfile(fname::AbstractString, nblks::Int)
    println("genrandfile starting at $(myid())")
    dfs = HDFSClient(TeraSort.DFSHOST, TeraSort.DFSPORT, UserGroupInformation(TeraSort.DFSUSER))
    open(dfs, "$fname-$(myid())", "w") do fptr
        for blk in 1:nblks
            for rec in 1:RECSPERBLK
                write(fptr, randstring(RECSZ-1), '\n')
            end
        end
    end
    println("genrandfile created $nblks blocks at $(myid())")
    nblks
end

function distgenrandfile(fname::AbstractString, nblkspernode::Int)
    wrkrs = workers()
    nwrkrs = length(wrkrs)
    refs = Array(RemoteRef, nwrkrs)
    idx = 1
    for nid in wrkrs
        refs[idx] = @spawnat nid TeraSort.genrandfile(fname, nblkspernode)
        idx += 1
    end
    nblks = 0
    for idx in 1:nwrkrs
        nblks += fetch(refs[idx])
    end
    nblks
end

function teragen()
    nblks,t,b,g = @timed distgenrandfile(joinpath(INDIR, "randstr"), 2)
    println("created $nblks blocks at 2 blocks per node")
    println("time:$t, memused:$b, gctime:$g")
end


text:
#=
    using the machinery NamedTupleTools

    prototypes
      NamedTuple{(:a, :b, :c),T} where T<:Tuple

    nt.fieldname
=#

#=
    using the machinery of Hyperscript supporting tag elaboration
=#

#=
   sorts of Content
=#

@enum AkoContent noContent=1 akoTextual akoImage akoTable akoListed akoForm akoInteger akoFloat akoNumber

const akoContent = Dict(:nothing=>noContent, :Nothing=>noContent,
    :text=>akoTextual, :string=>akoTextual, :AbstractString=>akoTextual,
    :img=>akoImage, :image=>akoImage, :form=>akoForm, :table=>akoTable, :ol=>akoListed, :ul=>akoListed,
    :Int8=>akoInteger, :Int16=>akoInteger, :Int64=>akoInteger, :Int64=>akoInteger, :Integer=>akoInteger,
    :Float16=>akoFloat, :Float32=>akoFloat, :Float64=>akoFloat, :AbstractFloat=>akoFloat, :Number=>akoNumber)

AkoContent(x::Symbol) = get(akoContent, x, nothing)

mutable struct Content{T} <: AbstractContent
    ako::AkoContent
    value::T

    function Content(value::T) where {T}
        ako = AkoContent(Symbol(T))
        if isnothing(ako)
            ako = AkoContent(Symbol(supertype(T)))
        end
        isnothing(ako) && println("\n  >>>  missing AkoContent for ",T," <: ", supertype(T), "\n")
        return new{T}(ako, value)
    end
end

const MaybeContent{T} = Union{Nothing, Content{T}}
iscontent(x::MaybeContent{T}) where {T} = !isnothing(x)
iscontent(x::Content{T}) where {T} = true
iscontent(x) = false

#=
   sorts of Context
=#

@enum AkoContext begin
    noContext=1 akoBackground akoForeground
    akoEntry akoSequence akoTabular akoHierarchy
    akoColor akoFont akoWord akoText
    akoBorder akoMargin akoPadding
    akoBox akoColumn akoGrid akoFlex
    akoOutline akoTransform
end

const akoContext = Dict(:nothing=>noContext, :Nothing=>noContext,
    :body=>akoHierarchy, :head=>akoHierarchy,
    :table=>akoTabular, :ol=>akoSequence, :ul=>akoSequence, :form=>akoEntry,
    # https://www.w3.org/TR/css-color/
    # https://www.w3.org/TR/css-color-adjust-1/
    # https://www.w3.org/TR/css-color-4/
    :color=>akoColor, Symbol("alpha-value")=>akoColor, Symbol("rendering-intent")=>akoColor,
    :lab=>akoColor, :lch=>akoColor, :hwb=>akoColor, :hsl=>akoColor, :rgb=>akoColor, :hue=>akoColor,
    # https://www.w3.org/TR/css-text/
    Symbol("font-family")=>akoFont, Symbol("font-size")=>akoFont, Symbol("font-weight")=>akoFont,
    Symbol("font-variant")=>akoFont, Symbol("font-kerning")=>akoFont, Symbol("font-stretch")=>akoFont,
    Symbol("font-size-adjust")=>akoFont, Symbol("line-height")=>akoFont, Symbol("letter-spacing")=>akoFont,
    Symbol("full-size")=>akoFont, :small=>akoFont, Symbol("full-width")=>akoFont,
    # https://www.w3.org/TR/css-text/
    Symbol("text")=>akoText, Symbol("text-transform")=>akoText, Symbol("text-indent")=>akoText,
    Symbol("hanging-punctuation")=>akoText, Symbol("white-space")=>akoText,
    Symbol("tab-size")=>akoText, Symbol("line-break")=>akoText, :hyphens=>akoText,
    Symbol("overflow-wrap")=>akoText, Symbol("text-align")=>akoText,
    Symbol("text-align-all")=>akoText, Symbol("text-align-last")=>akoText, Symbol("text-justify")=>akoText,
    Symbol("text-justify-all")=>akoText, Symbol("inter-character")=>akoText,
    # https://www.w3.org/TR/css-text-decor-3/
    #
    :word=>akoWord, Symbol("word-break")=>akoWord, Symbol("word-spacing")=>akoWord,  Symbol("word-wrap")=>akoWord,
    Symbol("keep-all")=>akoWord,  Symbol("inter-word")=>akoText,
    # https://www.w3.org/TR/css-lists-3/
    :outline=>akoOutline,
    :padding=>akoPadding, :margin=>akoMargin, :border=>akoBorder,
    :box=>akoBox, :column=>akoColumn, :grid=>akoGrid, :flex=>akoFlex,
    # https://www.w3.org/TR/css-images-3/
    :transform=>akoTransform
    # https://www.w3.org/TR/css-contain-1/ "This CSS module describes the contain property, which indicates that the element’s subtree is independent of the rest of the page."
    # https://www.w3.org/TR/css-contain-2/
    # https://www.w3.org/TR/css-values-3/
    # https://www.w3.org/TR/css-sizing-3/ "This module extends the CSS sizing properties with keywords that represent content-based "intrinsic" sizes and context-based "extrinsic" sizes, allowing CSS to more easily describe boxes that fit their content or fit into a particular layout context."
    )

AkoContext(x::Symbol) = get(akoContext, x, nothing)

mutable struct Context{T} <: AbstractContext
    ako::AkoContext
    value::T

    function Context(value::T) where {T}
        ako = AkoContext(Symbol(T))
        if isnothing(ako)
            ako = AkoContext(Symbol(supertype(T)))
        end
        isnothing(ako) && println("\n  >>>  missing AkoContext for ",T," <: ", supertype(T), "\n")
        return new{T}(ako, value)
    end
end

const MaybeContext{T} = Union{Nothing, Context{T}}
iscontext(x::MaybeContext{T}) where {T} = !isnothing(x)
iscontext(x::Context{T}) where {T} = true
iscontext(x) = false


text:
"""
    Duality(optimizer)

Duality uses Lagrangian relaxation to compute over-approximated bounds for a network

# Problem requirement
1. Network: any depth, any activation function that is monotone
2. Input: hyperrectangle
3. Output: halfspace

# Return
`BasicResult`

# Method
Lagrangian relaxation.
Default `optimizer` is `GLPKSolverMIP()`.

# Property
Sound but not complete.

# Reference
[K. Dvijotham, R. Stanforth, S. Gowal, T. Mann, and P. Kohli,
"A Dual Approach to Scalable Verification of Deep Networks,"
*ArXiv Preprint ArXiv:1803.06567*, 2018.](https://arxiv.org/abs/1803.06567)
"""
@with_kw struct Duality{O<:AbstractMathProgSolver}
    optimizer::O = GLPKSolverMIP()
end

# can pass keyword args to optimizer
# Duality(optimizer::DataType = GLPKSolverMIP; kwargs...) =  Duality(optimizer(kwargs...))

function solve(solver::Duality, problem::Problem)
    model = Model(solver = solver.optimizer)
    c, d = tosimplehrep(problem.output)
    λ = init_multipliers(model, problem.network)
    μ = init_multipliers(model, problem.network)
    o = dual_value(solver, problem, model, λ, μ)
    @constraint(model, last(λ) .== -c)
    status = solve(model, suppress_warnings = true)
    return interpret_result(solver, status, o - d[1])
end

# False if o > 0, True if o <= 0
function interpret_result(solver::Duality, status, o)
    status != :Optimal && return BasicResult(:Unknown)
    getvalue(o) <= 0.0 && return BasicResult(:SAT)
    return BasicResult(:UNSAT)
end

function dual_value(solver::Duality,
                    problem::Problem,
                    model::Model,
                    λ::Vector{Vector{Variable}},
                    μ::Vector{Vector{Variable}})
    bounds = get_bounds(problem)
    layers = problem.network.layers
    # input layer
    o = input_layer_value(layers[1], μ[1], bounds[1])
    o += activation_value(layers[1], μ[1], λ[1], bounds[1])
    # other layers
    for i in 2:length(layers)
        o += layer_value(layers[i], μ[i], λ[i-1], bounds[i])
        o += activation_value(layers[i], μ[i], λ[i], bounds[i])
    end
    @objective(model, Min, o)
    return o
end

function activation_value(layer::Layer,
                          μᵢ::Vector{Variable},
                          λᵢ::Vector{Variable},
                          bound::Hyperrectangle)
    o = zero(eltype(μᵢ))
    b_hat = linear_transformation(layer, bound)
    l_hat, u_hat = low(b_hat), high(b_hat)
    l, u = layer.activation(l_hat), layer.activation(u_hat)

    o += sum(symbolic_max.(μᵢ.*l_hat, μᵢ.*u_hat))
    o += sum(symbolic_max.(λᵢ.*l, λᵢ.*u))

    # also possible
    # o += sum(@. symbolic_max(μᵢ*l_hat, μᵢ*u_hat))
    # o += sum(@. symbolic_max(λᵢ*l, λᵢ*u))

    # original
    # for j in 1:length(l)
    #     o += symbolic_max(μᵢ[j]*l_hat[j], μᵢ[j]*u_hat[j])
    #     o += symbolic_max(λᵢ[j]*l[j], λᵢ[j]*u[j])
    # end
    return o
end

function layer_value(layer::Layer,
                     μᵢ::Vector{Variable},
                     λᵢ::Vector{Variable},
                     bound::Hyperrectangle)
    (W, b) = (layer.weights, layer.bias)
    o = λᵢ' * bound.center - μᵢ' * (W * bound.center + b)
    # instead of for-loop:
    o += sum(symbolic_abs.(λᵢ .- W'*μᵢ) .* bound.radius) # TODO check that this is equivalent to before
    return o
end

function input_layer_value(layer::Layer,
                           μᵢ::Vector{Variable},
                           input::Hyperrectangle)
    W, b = layer.weights, layer.bias
    o = -μᵢ' * (W*input.center .+ b)
    o += sum(symbolic_abs.(μᵢ'*W) .* input.radius)   # TODO check that this is equivalent to before
    return o
end


text:
using HMCUtilities
using Test
using LinearAlgebra

include("utility.jl")
include("constraints.jl")


text:
using Statistics
using NCDatasets
using Formatting

using ArgParse
using JSON

include("LinearRegression.jl")
include("nanop.jl")

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Ocean data file. New variable will be appended."
            arg_type = String
            required = true

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

        "--output-file"
            help = "Output 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

        "--varname"
            help = "Variable name. If it is 2D then it should be (lon, lat, time), if it is 3D then it should be (lon, lat, lev/depth, time)"
            arg_type = String
            required = true

        "--detrend"
            help = "Whether to detrend the data or not."
            action = :store_true

    end

    return parse_args(ARGS, s)
end

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

output_file = parsed["output-file"]

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

    dims = size(ds[parsed["varname"]])

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

    if length(dims) == 3 # 2D + time
        rng = (:,:,beg_t:end_t)
        z_exists = false
    elseif length(dims) == 4 # 3D + time
        rng = (:,:,:,beg_t:end_t)
        z_exists = true
    else
        ErrorException("Unknown dimension") |> throw
    end

    global data  = replace(ds[parsed["varname"]][rng...], missing=>NaN)

    if z_exists
        global (Nx, Ny, Nz, Nt) = size(data)
    else
        global (Nx, Ny, Nt) = size(data)
        Nz = 1
        data = reshape(data, Nx, Ny, Nz, Nt)
    end

    if mod(Nt, 12) != 0
        ErrorException("Time record is not multiple of 12") |> throw
    end

    global nyears = Int64(Nt / 12)
end

Dataset(parsed["domain-file"], "r") do ds
    global mask = replace(ds["mask"], missing=>NaN)
end

# MM   = Monthly Mean
# A    = Anomaly
# YYC  = Year-to-Year Correlation

data = nanmean(data, dims=(1,))[1, :, :, :]

data_MM    = zeros(Float64, Ny, Nz, 12)
data_MA    = zeros(Float64, Ny, Nz, Nt)
data_MAVAR = zeros(Float64, Ny, Nz, 12)

x = collect(Float64, 1:Nt)
for j=1:Ny, k=1:Nz

    d = view(data, j, k, :)
    if parsed["detrend"]
        d = detrend(x, d)
    end

    data_MM[j, k, :] = mean( reshape(d, 12, :), dims=2 )[:, 1]

#    println("size: ", size(data_MA))

    data_MA[j, k, :] = d - repeat( data_MM[j, k, :], outer=nyears)

    for m = 1:12
        d_yy = view(data_MA, j, k, m:12:(m+nyears*12-1))
        data_MAVAR[j, k, m] = std(d_yy)
    end

end

Dataset(output_file, "c") do ds

    defDim(ds, "months", 12)
    defDim(ds, "time", Inf)
    defDim(ds, "Ny", Ny)
    defDim(ds, "Nz", Nz)

    for (varname, vardata, vardim, attrib) in [
        (format("{:s}_MM",    parsed["varname"]), data_MM,    ("Ny", "Nz", "months"), Dict()),
        (format("{:s}_MA",    parsed["varname"]), data_MA,    ("Ny", "Nz", "time"),   Dict()),
        (format("{:s}_MAVAR", parsed["varname"]), data_MAVAR, ("Ny", "Nz", "months"), Dict()),
    ]

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

        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:
using DiffEqJump, DiffEqBase, OrdinaryDiffEq, Statistics
using Test

rate = (u,p,t) -> u
affect! = function (integrator)
  integrator.u += 1
end
jump = ConstantRateJump(rate,affect!)

rate = (u,p,t) -> 0.5u
affect! = function (integrator)
  integrator.u -= 1
end
jump2 = ConstantRateJump(rate,affect!)

prob = DiscreteProblem(1.0,(0.0,3.0))
jump_prob = JumpProblem(prob,Direct(),jump)

sol = solve(jump_prob,FunctionMap())

# using Plots; plot(sol)

prob = DiscreteProblem(10.0,(0.0,3.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)

sol = solve(jump_prob,FunctionMap())

# plot(sol)

nums = Int[]
@time for i in 1:10000
  local jump_prob = JumpProblem(prob,Direct(),jump,jump2)
  local sol = solve(jump_prob,FunctionMap())
  push!(nums,sol[end])
end

@test mean(nums) - 45 < 1

prob = DiscreteProblem(1.0,(0.0,3.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)

sol = solve(jump_prob,FunctionMap())

nums = Int[]
@time for i in 1:10000
  local jump_prob = JumpProblem(prob,Direct(),jump,jump2)
  local sol = solve(jump_prob,FunctionMap())
  push!(nums,sol[2])
end

@test sum(nums .== 0)/10000 - 0.33 < 0.02


text:
module SpringCollab2020BlockJellySpawnTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/BlockJellySpawnTrigger" BlockJellySpawnTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight)

const placements = Ahorn.PlacementDict(
    "Block Jelly Spawn (Spring Collab 2020)" => Ahorn.EntityPlacement(
        BlockJellySpawnTrigger,
        "rectangle"
    )
)

end


text:
"""
    `showmap(run)` plotting of map for a given data set.
"""
function showmap(run::Run;
    minzoom::Integer = 0, maxzoom::Integer = 21,
    mapservice::AbstractString="google", maptype::AbstractString="roadmap")
    mapcoor = mapprojection(run, minzoom = minzoom, maxzoom = maxzoom)
    mapfile = downloadmap(mapcoor, mapservice = mapservice, maptype = maptype)
    runmap = permuteddimsview(rawview(channelview(load(mapfile))), (2, 3, 1)) / 256
    xylim = [mapcoor.lon_left, mapcoor.lon_right, mapcoor.lat_bottom, mapcoor.lat_top]
    xticks([])
    yticks([])
    imshow(runmap, origin="upper", extent=xylim, aspect=1.5)
end

"""
    `plotrun(run::Run)` plot running data.
"""
function plotrun(run::Run; s::AbstractFloat = 1.0)
    kmh = speed_kmh(run)
    kmh = kmh / maximum(kmh)
    scatter(run.lon[2:end], run.lat[2:end], c = kmh, cmap = ColorMap(:coolwarm), s = s)
end


text:
struct TreeField
    width::UInt32
    height::UInt32
    data::Array{Bool, 1}
end

function fromStdin()::TreeField
    height = 0
    width = 0
    data = Bool[]
    for line in eachline()
        height += 1
        width = length(line)
        for i in firstindex(line):lastindex(line)
            push!(data, line[i] != '.')
        end
    end
    TreeField(width, height, data)
end

function getIndex(field::TreeField, x, y)::Bool
    field.data[(y - 1) * field.width + mod1(x, field.width)]
end

function main()
    field = fromStdin()
    x = 1
    y = 1
    trees = 0
    while y <= field.height
        if getIndex(field, x, y)
            trees += 1
        end
        y += 1
        x += 3
    end

    println(trees)
end

main()


text:
# performance testing

using Formatting

# performance of format parsing

fp1() = FormatExpr("{1} + {2} + {3}")

fp1()
@time for i=1:10000; fp1(); end

fp2() = FormatExpr("abc {1:*<5d} efg {2:*>8.4e} hij")

fp2()
@time for i=1:10000; fp2(); end

# performance of string formatting

const fe1 = fp1()
const fe2 = fp2()

sf1(x, y, z) = format(fe1, x, y, z)
sf1(10, 20, 30)
@time for i=1:10000; sf1(10, 20, 30); end

sf2(x, y) = format(fe2, x, y)
sf2(10, 20)
@time for i=1:10000; sf2(10, 20); end


text:
mutable struct Man
    name::String

    println("Initialized!")
end

function hello(self::Man)
    println("Hello ", self.name, "!")
end

function goodbye(self::Man)
    println("Good-bye ", self.name, "!")
end

function main()
    m = Man("David")
    hello(m)
    goodbye(m)
end

main()


text:
using Pkg, DocumentationGenerator
docspath = normpath(joinpath(@__DIR__, ".."))
packages = sort!(readdir(joinpath(docspath, "build")), by=x->lowercase(x))

badge_regex = r"(\(|^| )(https://[^\)\n]*?svg[^\)\n]*?)((\))|($)|( ))"ms
known_badges = [
    "appveyor",
    "codecov",
    "coveralls",
    "travis",
    "circle",
    "gitlab"
]

function badge_value(badge)
    try
        svg = read(download(badge), String)
        if occursin("passing", svg)
            return "passes"
        elseif occursin("failing", svg)
            return "failing"
        else occursin("coverage", svg)
            m = match(r">(\d\d)%<", svg)
            return m === nothing ? "unknown" : parse(Int, m[1])
        end
    catch err
        @error "error while retrieving CI badge" error=err
        return "unknown"
    end
end

function extract_badges(root)
    src = joinpath(root, "_packagesource")
    readme = joinpath(src, "README.md")
    if isfile(readme)
        readme_str = read(readme, String)
        urls = getindex.(eachmatch(badge_regex, readme_str), 2)
        return Dict(map(known_badges) do badge
            res = findall(url-> occursin(badge, url), urls)
            badge => if isempty(res)
                []
            else
                map(res) do idx
                    Dict("url" => urls[idx], "value" => badge_value(urls[idx]))
                end
            end
        end)
    else
        println(basename(root), " doesn't contain readme")
        return Dict{String, Any}(Pair.(known_badges, ([],)))
    end
end

function extract_info_dict(pkg)
    results = Dict{String, Any}()
    build = joinpath(@__DIR__, "..", "build", pkg)
    root = joinpath(build, readdir(build)[1])
    results["badges"] = extract_badges(root)
    if isfile(root, "meta.toml")
        toml = Pkg.TOML.parsefile(joinpath(root, "meta.toml"))
        merge!(results, toml)
    end
    results
end

package_info = Dict{String, Any}()
for pkg in packages
    package_info[pkg] = extract_info_dict(pkg)
end

using JSON

open(joinpath(@__DIR__, "..", "build", "packages_meta.json"), "w") do io
    JSON.print(io, package_info)
end


text:
"""
    perturb_input(st_snow, prec, tair, itime)

Perturb input data for snow model
"""
function perturb_input(st_snow, prec, tair, itime)

    n = Uniform(0.5, 1.5)
    prec_noise = rand(n, 1)

    n = Normal(0.0, 2)
    tair_noise = rand(n, 1)

    # Assign inputs to snow model

    for izone in eachindex(st_snow.prec)

        st_snow.prec[izone] = prec[izone, itime] * prec_noise[1]
        st_snow.tair[izone] = tair[izone, itime] + tair_noise[1]

    end

end

"""
    particle_filter(st_snow, st_hydro, prec, tair, epot, q_obs, nens)

Run the particle filter for any combination of snow and hydrological routing model.
"""
function particle_filter(st_snow, st_hydro, prec, tair, epot, q_obs, nens; test_forecast = false)

    srand(1)

    # Number of elevation bands (rows) and time steps (cols)

    nzones = size(prec, 1)
    ntimes = size(prec, 2)

    # Initilize state variables

    st_snow  = [deepcopy(st_snow) for i in 1:nens]
    st_hydro = [deepcopy(st_hydro) for i in 1:nens]

    # Initilize particles

    wk = ones(nens) / nens

    # Allocate arrays

    q_res  = zeros(ntimes, 3)
    q_sim  = zeros(1, nens)

    # Arrays for forecast

    if test_forecast
        q_sim_forecast = []
        q_obs_forecast = []
    end

    # Start time loop

    for itime = 1:ntimes

        for ipart = 1:nens

            perturb_input(st_snow[ipart], prec, tair, itime)

            run_timestep(st_snow[ipart])

            get_input(st_snow[ipart], st_hydro[ipart], epot, itime)

            run_timestep(st_hydro[ipart])

            q_sim[ipart] = st_hydro[ipart].q_sim

        end

        # Run particle filter

        if true

            for ipart = 1:nens

                wk[ipart] = pdf(Normal(q_obs[itime], max(0.1 * q_obs[itime], 0.1)), q_sim[ipart]) * wk[ipart]

            end

            if sum(wk) > 0.0
                wk = wk / sum(wk)
            else
                wk = ones(nens) / nens
            end

            # Perform resampling

            Neff = 1 / sum(wk.^2)

            if round(Int64, Neff) < round(Int64, nens * 0.8)

                println("Resampled at step: $itime")

                indx = resample(wk)

                st_snow  = [deepcopy(st_snow[i]) for i in indx]
                st_hydro = [deepcopy(st_hydro[i]) for i in indx]

                wk = ones(nens) / nens

            end

        end

        # Store results

        q_res[itime, 1] = sum(wk[:] .* q_sim[:])
        q_res[itime, 2] = minimum(q_sim)
        q_res[itime, 3] = maximum(q_sim)

        # Run forecast

        if test_forecast

            tmp_sim, tmp_obs = run_forecast(st_snow, st_hydro, prec, tair, epot, q_obs, itime, nens, 7)

            if isempty(q_sim_forecast)
                q_sim_forecast = tmp_sim
                q_obs_forecast = tmp_obs
            else
                q_sim_forecast = hcat(q_sim_forecast, tmp_sim)
                q_obs_forecast = hcat(q_obs_forecast, tmp_obs)
            end

        end

    end

    # Return results

    if test_forecast
        q_sim_forecast = transpose(q_sim_forecast)
        q_obs_forecast = transpose(q_obs_forecast)
        return q_res, q_sim_forecast, q_obs_forecast
    else
        return q_res
    end

end

"""
    enkf_filter(st_snow, st_hydro, prec, tair, epot, q_obs, nens)

Run the ensemble Kalman filter for any combination of snow and hydrological
routing model.
"""
function enkf_filter(st_snow, st_hydro, prec, tair, epot, q_obs, nens; test_forecast = false)

    srand(1)

    # Number of elevation bands (rows) and time steps (cols)

    nzones = size(prec, 1)
    ntimes = size(prec, 2)

    # Initilize state variables

    st_snow  = [deepcopy(st_snow) for i in 1:nens]
    st_hydro = [deepcopy(st_hydro) for i in 1:nens]

    # Allocate arrays

    q_res  = zeros(ntimes, 3)
    q_sim  = zeros(1, nens)

    # Arrays for forecast

    if test_forecast
        q_sim_forecast = []
        q_obs_forecast = []
    end

    # Start time loop

    for itime = 1:ntimes

        # Perturb inputs and run models

        for iens = 1:nens

            perturb_input(st_snow[iens], prec, tair, itime)

            run_timestep(st_snow[iens])

            get_input(st_snow[iens], st_hydro[iens], epot, itime)

            run_timestep(st_hydro[iens])

            q_sim[iens] = st_hydro[iens].q_sim

        end

        # Run filter part

        if q_obs[itime] >= 0

            # Perturb observations

            sigma = max(0.1 * q_obs[itime], 0.1)

            obs_ens = q_obs[itime] + sigma * randn(Float64, 1, nens)

            # Update states of hydrological model

            enkf_hydro(st_hydro, obs_ens, q_sim)

            # Update states of snow model

            enkf_snow(st_snow, obs_ens, q_sim)

            # Update simulated discharge

            q_sim  = enkf(q_sim, obs_ens, q_sim)

            q_sim[q_sim .< 0] = 0.

        end

        # Store results

        q_res[itime, 1] = mean(q_sim)
        q_res[itime, 2] = minimum(q_sim)
        q_res[itime, 3] = maximum(q_sim)

        # Run forecast

        if test_forecast

            tmp_sim, tmp_obs = run_forecast(st_snow, st_hydro, prec, tair, epot, q_obs, itime, nens, 7)

            if isempty(q_sim_forecast)
                q_sim_forecast = tmp_sim
                q_obs_forecast = tmp_obs
            else
                q_sim_forecast = hcat(q_sim_forecast, tmp_sim)
                q_obs_forecast = hcat(q_obs_forecast, tmp_obs)
            end

        end

    end

    # Return results

    if test_forecast
        q_sim_forecast = transpose(q_sim_forecast)
        q_obs_forecast = transpose(q_obs_forecast)
        return q_res, q_sim_forecast, q_obs_forecast
    else
        return q_res
    end

end

"""
    run_forecast(st_snow, st_hydro, prec, tair, epot, q_obs, itime, nens, ndays)

Run forecast
"""
function run_forecast(st_snow, st_hydro, prec, tair, epot, q_obs, itime, nens, ndays)

    # Run forecast for ndays

    q_sim_forecast = [NaN for i in 1:ndays]
    q_obs_forecast = [NaN for i in 1:ndays]

    if itime < (length(q_obs) - ndays - 1)

        st_snow  = [deepcopy(st_snow[i]) for i in 1:nens]
        st_hydro = [deepcopy(st_hydro[i]) for i in 1:nens]

        q_sim = zeros(nens)

        for itmp in (itime+1):(itime+ndays)

            # Perturb inputs and run models

            for iens = 1:nens

                perturb_input(st_snow[iens], prec, tair, itmp)

                run_timestep(st_snow[iens])

                get_input(st_snow[iens], st_hydro[iens], epot, itmp)

                run_timestep(st_hydro[iens])

                q_sim[iens] = st_hydro[iens].q_sim

            end

            q_sim_forecast[itmp-itime] = mean(q_sim)
            q_obs_forecast[itmp-itime] = q_obs[itmp]

        end

    end

    return q_sim_forecast, q_obs_forecast

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann