text:
## GridSpec function with default gridName argument:

GridSpec() = GridSpec("LLC90")

## GridSpec function with gridName argument:

"""
    GridSpec(gridName)

Return a `gmcgrid` specification that provides grid files `path`,
`class`, `nFaces`, `ioSize`, `facesSize`, `ioPrec`, & a `read` function
(not yet) using hard-coded values for `"LLC90"`, `"CS32"`, `"LL360"` (for now).
"""
function GridSpec(gridName,gridParentDir="./")

if gridName=="LLC90";
    grDir=gridParentDir*"GRID_LLC90/";
    nFaces=5;
    grTopo="llc";
    ioSize=[90 1170];
    facesSize=[(90, 270), (90, 270), (90, 90), (270, 90), (270, 90)]
    ioPrec=Float64;
elseif gridName=="CS32";
    grDir=gridParentDir*"GRID_CS32/";
    nFaces=6;
    grTopo="cs";
    ioSize=[32 192];
    facesSize=[(32, 32), (32, 32), (32, 32), (32, 32), (32, 32), (32, 32)]
    ioPrec=Float32;
elseif gridName=="LL360";
    grDir=gridParentDir*"GRID_LL360/";
    nFaces=1;
    grTopo="ll";
    ioSize=[360 160];
    facesSize=[(360, 160)]
    ioPrec=Float32;
else;
    error("unknown gridName case");
end;

mygrid=gcmgrid(grDir,grTopo,nFaces,facesSize, ioSize, ioPrec, read, write)

return mygrid

end

## GridLoad function

"""
    GridLoad(mygrid::gcmgrid)

Return a `Dict` of grid variables read from files located in `mygrid.path` (see `?GridSpec`).

Based on the MITgcm naming convention, grid variables are:

- XC, XG, YC, YG, AngleCS, AngleSN, hFacC, hFacS, hFacW, Depth.
- RAC, RAW, RAS, RAZ, DXC, DXG, DYC, DYG.
- DRC, DRF, RC, RF (one-dimensional)
"""
function GridLoad(mygrid::gcmgrid)

    GridVariables=Dict()

    list0=("XC","XG","YC","YG","AngleCS","AngleSN","RAC","RAW","RAS","RAZ",
    "DXC","DXG","DYC","DYG","Depth")
    for ii=1:length(list0)
        tmp1=mygrid.read(mygrid.path*list0[ii]*".data",MeshArray(mygrid,mygrid.ioPrec))
        tmp2=Symbol(list0[ii])
        @eval (($tmp2) = ($tmp1))
        GridVariables[list0[ii]]=tmp1
    end

    mygrid.ioPrec==Float64 ? reclen=8 : reclen=4

    list0=("DRC","DRF","RC","RF")
    for ii=1:length(list0)
        fil=mygrid.path*list0[ii]*".data"
        tmp1=stat(fil)
        n3=Int64(tmp1.size/reclen)

        fid = open(fil)
        tmp1 = Array{mygrid.ioPrec,1}(undef,n3)
        read!(fid,tmp1)
        tmp1 = hton.(tmp1)

        tmp2=Symbol(list0[ii])
        @eval (($tmp2) = ($tmp1))
        GridVariables[list0[ii]]=tmp1
    end

    list0=("hFacC","hFacS","hFacW")
    n3=length(GridVariables["RC"])
    for ii=1:length(list0)
        tmp1=mygrid.read(mygrid.path*list0[ii]*".data",MeshArray(mygrid,mygrid.ioPrec,n3))
        tmp2=Symbol(list0[ii])
        @eval (($tmp2) = ($tmp1))
        GridVariables[list0[ii]]=tmp1
    end

    return GridVariables

end

"""
    GridOfOnes(grTp,nF,nP)

Define all-ones grid variables instead of using `GridSpec` & `GridLoad`.
"""
function GridOfOnes(grTp,nF,nP)

    grDir=""
    grTopo=grTp
    nFaces=nF
    ioSize=[nP nP*nF]
    facesSize=Array{NTuple{2, Int},1}(undef,nFaces)
    facesSize[:].=[(nP,nP)]
    ioPrec=Float32

    mygrid=gcmgrid(grDir,grTopo,nFaces,facesSize, ioSize, ioPrec, read, write)

    GridVariables=Dict()
    list0=("XC","XG","YC","YG","RAC","RAZ","DXC","DXG","DYC","DYG","hFacC","hFacS","hFacW","Depth");
    for ii=1:length(list0);
        tmp1=fill(1.,nP,nP*nF);
        tmp1=mygrid.read(tmp1,MeshArray(mygrid,Float64));
        tmp2=Symbol(list0[ii]);
        @eval (($tmp2) = ($tmp1))
        GridVariables[list0[ii]]=tmp1
    end

    return GridVariables

end

"""
    findtiles(ni::Int,nj::Int,mygrid::gcmgrid)
    findtiles(ni::Int,nj::Int,grid::String="llc90",gridParentDir="./")

Return a `MeshArray` map of tile indices for tile size `ni,nj`
"""
function findtiles(ni::Int,nj::Int,mygrid::gcmgrid)
    mytiles = Dict()

    GridVariables=GridLoad(mygrid)

    mytiles["nFaces"]=mygrid.nFaces;
    #mytiles.fileFormat=mygrid.fileFormat;
    mytiles["ioSize"]=mygrid.ioSize;

    XC=GridVariables["XC"];
    YC=GridVariables["YC"];
    XC11=copy(XC); YC11=copy(XC);
    XCNINJ=copy(XC); YCNINJ=copy(XC);
    iTile=copy(XC); jTile=copy(XC); tileNo=copy(XC);
    tileCount=0;
    for iF=1:XC11.grid.nFaces
        #global tileCount,XC,YC,XC11,YC11,iTile,jTile,tileNo
        face_XC=XC.f[iF]; face_YC=YC.f[iF];
    #ordering convention that was used in first generation nctile files:
    #    for ii=1:size(face_XC,1)/ni;
    #        for jj=1:size(face_XC,2)/nj;
    #ordering convention that is consistent with MITgcm/pkg/exch2:
        for jj=Int.(1:size(face_XC,2)/nj);
            for ii=Int.(1:size(face_XC,1)/ni);
                tileCount=tileCount+1;
                tmp_i=(1:ni).+ni*(ii-1)
                tmp_j=(1:nj).+nj*(jj-1)
                tmp_XC=face_XC[tmp_i,tmp_j]
                tmp_YC=face_YC[tmp_i,tmp_j]
                XC11.f[iF][tmp_i,tmp_j].=tmp_XC[1,1]
                YC11.f[iF][tmp_i,tmp_j].=tmp_YC[1,1]
                XCNINJ.f[iF][tmp_i,tmp_j].=tmp_XC[end,end]
                YCNINJ.f[iF][tmp_i,tmp_j].=tmp_YC[end,end]
                iTile.f[iF][tmp_i,tmp_j]=collect(1:ni)*ones(Int,1,nj)
                jTile.f[iF][tmp_i,tmp_j]=ones(Int,ni,1)*collect(1:nj)'
                tileNo.f[iF][tmp_i,tmp_j]=tileCount*ones(Int,ni,nj)
            end
        end
    end

    mytiles["XC"] = XC;
    mytiles["YC"] = YC;
    mytiles["XC11"] = XC11;
    mytiles["YC11"] = YC11;
    mytiles["XCNINJ"] = XCNINJ;
    mytiles["YCNINJ"] = YCNINJ;
    mytiles["iTile"] = iTile;
    mytiles["jTile"] = jTile;
    mytiles["tileNo"] = tileNo;

    return mytiles

end

findtiles(ni::Int,nj::Int,grid::String="llc90",gridParentDir="./") = findtiles(ni,nj,GridSpec(grid,gridParentDir))


text:
include("WaveOptics.jl")

function psf(;
    objective_mag,
    NA,
    λ,
    medium_index,
    f_tubelens,
    camera_pixel_size,
    psf_shape,
    zstep,
    oversampling,
)
    camera_num_pixels = max(psf_shape[1], psf_shape[2])
    z_num_pixels = psf_shape[3]

    # Simulate a PSF as a square image
    NA = NA > medium_index ? medium_index : NA
    recipe = WaveOptics.WideFieldOpticalRecipe(
        Float64(objective_mag),
        Float64(NA),
        Float64(λ),
        Float64(medium_index),
        Float64(f_tubelens),
    )

    sim_size_px = camera_num_pixels * oversampling
    sim_size_m = camera_num_pixels * camera_pixel_size / objective_mag

    model = WaveOptics.DefocusModel(recipe, sim_size_m, sim_size_px)

    # Compute the 3D PSF at the sampling resolution
    zmax = zstep * (z_num_pixels - 1) * 0.5
    zrange = range(-zmax, step = zstep, length = z_num_pixels)

    psf = WaveOptics.psf_3d(model, zrange)
    psf_sensor = WaveOptics.downsample_to_sensor_resolution(psf, oversampling)

    # Crop the PSF to the specified size (psf_shape)
    d = abs(psf_shape[1] - psf_shape[2])
    p1 = div(d, 2)
    p2 = d - p1
    if psf_shape[1] != camera_num_pixels
        psf_sensor = psf_sensor[p1+1:end-p2, :, :]
    elseif psf_shape[2] != camera_num_pixels
        psf_sensor = psf_sensor[:, p1+1:end-p2, :]
    end

    factor = sum(psf_sensor[:, :, div(z_num_pixels, 2)])

    return psf_sensor ./ factor
end

ext: jl


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using CompilerTools
using CompilerTools.LambdaHandling
using Base.Test

## Tests for CompilerTools.AliasAnalysis

function test_alias_1(x::Int64, y::Int64, z::Int64)
    A = rand(x, y)
    D_arr = rand(y, z)
    C = zeros(x, z)
    B = D_arr

    for i = 1:y
        if( i%2 == 0 )
            for j = 1:z
                D_arr[i,j] /= 2
            end
        end
    end

    for i = 1:x
        for j = 1:y
            A[i,j] += 1;
        end
        for j = 1:z
            for k = 1:y
                C[i,j] += A[i,k] * B[k,j] + 1
            end
        end
    end
end

ast = code_typed(test_alias_1, (Int64,Int64,Int64,))[1]
linfo, body = lambdaToLambdaVarInfo(ast)
#println("ast = ", ast)
#println("linfo = ", linfo)

#cfg_2 = CompilerTools.CFGs.from_lambda(ast) :: CompilerTools.CFGs.CFG

#CompilerTools.AliasAnalysis.set_debug_level(3)
#CompilerTools.CFGs.set_debug_level(4)
#CompilerTools.LivenessAnalysis.set_debug_level(5)
#CompilerTools.LambdaHandling.set_debug_level(5)

lives = CompilerTools.LivenessAnalysis.from_lambda(linfo, body)
#println("lives = ")
#println(lives)
handled = CompilerTools.AliasAnalysis.from_lambda(linfo, body, lives)
handled = map(x -> lookupVariableName(x, linfo), handled)
#println("handled = ", handled)

@test (in(:A, handled))
@test (in(:D_arr, handled) == false)
@test (in(:B, handled) == false)
#TODO: C is not detected now.
#@test (in(:C, handled))


text:
function _val(m::nmod, p)
  if m == 0
    return inf
  end
  return valuation(lift(m), p)
end

function _min_val(M, p)
  L = Union{PosInf, Int}[_val(a, p) for a in M]
  return minimum(L)
end

@doc Markdown.doc"""
    _last_block_index(G::Union{nmod_mat, fmpz_mod_mat}, p) -> Int, Int, Int

Return the starting index of the last modular block, as well as its valuation
and the valuation of the second to last modular block.

Assumes that $G$ is a block diagonal matrix.
"""
function _last_block_index(G::Union{nmod_mat, fmpz_mod_mat}, p)
  n = nrows(G)
  val = _min_val(G[n,:], p)
  val_current = val
  for k in 1:n-1
    val_current = _min_val(G[n-k,:], p)
    if val != val_current
      return n-k+1, val, val_current
    end
  end
  return 1, val, val_current
end

@doc Markdown.doc"""
    _hensel_qf(Z::T, G::T, F::T, a, b, p) where T <: Union{nmod_mat, fmpz_mod_mat}

The real worker for `hensel_qf`. Input is
  - $Z, G, F$, symmetric `n \times n` matrices.
  - $a, b$, integers with $a<b$.
We require that the triple $(Z, G, F)$ is $a$-adapted.

Return a matrix `F_l` such that $(Z, G, F_l)$ is $b$-adapted in particular $F
\equiv F_l \mod p^b$.
"""
function _hensel_qf(Z::T, G::T, F::T, a, b, p) where {T <: Union{nmod_mat, fmpz_mod_mat}}
  #@req _min_val(Z-F*G*F',p)>=a,"input must be adapted"
  i, s1, s2 = _last_block_index(G, p)
  R = base_ring(Z)
  Z = change_base_ring(R, divexact(lift(Z), p^s1))
  G = change_base_ring(R, divexact(lift(G), p^s1))
  s = s2 - s1
  if s == 0
    @assert i == 1
    s = b-a
  end
  even = p == 2
  Zn = Z[i:end,i:end] - F[i:end,1:i-1]*G[1:i-1,1:i-1]*F[i:end,1:i-1]'
  Gn = G[i:end,i:end]
  if even
    F[i:end,i:end] = _hensel_qf_modular_even(Zn, Gn, F[i:end,i:end], a, b)
  else
    F[i:end,i:end] = _hensel_qf_modular_odd(Zn, Gn, F[i:end,i:end], a, b)
  end
  K = inv(G[i:end,i:end]*F[i:end,i:end]')
  if i == 1
    return F
  end
  while a < b
    # an extra line that recomputes the upper block diagonal
    # if the input really is adapted this should not be necessary
    # but in any case it does not hurt
    F[1:i-1,i:end] = (Z[1:i-1,i:end] - F[1:i-1,1:i-1]*G[1:i-1,1:i-1]*F[i:end,1:i-1]') * K
    Zn = Z[1:i-1,1:i-1] - F[1:i-1,i:end]*G[i:end,i:end]*F[1:i-1,i:end]'
    F[1:i-1,1:i-1] = _hensel_qf(Zn, G[1:i-1,1:i-1], F[1:i-1,1:i-1], a, a+s, p)
    F[1:i-1,i:end] = (Z[1:i-1,i:end] - F[1:i-1,1:i-1]*G[1:i-1,1:i-1]*F[i:end,1:i-1]') * K
    a = a + s
  end
  return F
end

@doc Markdown.doc"""
    _hensel_qf_modular_odd(Z::T, G::T, F::T, a, b)
                                        where T <: Union{nmod_mat, fmpz_mod_mat}

Helper function for `_hensel_qf`. Matrices $Z$ and $G$ are assumed to be
modular symmetric matrices. We require that the triple $(Z,G,F)$ is
`a`-adapted.
"""
function _hensel_qf_modular_odd(Z::T, G::T, F::T, a, b) where {T <: Union{nmod_mat, fmpz_mod_mat}}
  while a < b
    Y = divexact(Z - F*G*F', 2)
    F = F + Y*inv(G*F')
    a = 2 * a
  end
  return F
end

function _solve_X(Y::Union{nmod_mat, fmpz_mod_mat}, b, g)
  F = FiniteField(2)
  Y = change_base_ring(F, lift(Y))
  b = [F(lift(i)) for i in b]
  g = [F(lift(i)) for i in g]
  return _solve_X(Y, b, g)
end

function _solve_X_ker(Y::Union{nmod_mat, fmpz_mod_mat}, b, g)
  F = FiniteField(2)
  Y = change_base_ring(F, lift(Y))
  b = [F(lift(i)) for i in b]
  g = [F(lift(i)) for i in g]
  return _solve_X_ker(Y, b, g)
end

function _solve_X_get_A_and_c(Y::gfp_mat, b, g)
  k = base_ring(Y)
  Y = matrix(k, nrows(Y), ncols(Y), [k(lift(a)) for a in Y])'

  @req issymmetric(Y) "Y must be symmetric"
  @req ncols(Y) == nrows(Y) "Y must be a square matrix"
  n = ncols(Y)

  equations = Vector{elem_type(k)}[]
  for i in 1:n
    R = zero_matrix(k, n, n)
    R[i,:] = g
    R[i,i] += 1
    eqn = reshape(collect(R'), :)
    push!(eqn, b[i])
    push!(equations, eqn)
  end
  # equations to a symmetric matrix
  for i in 1:n
    for j in i+1:n
      R = zero_matrix(k, n, n)
      R[i, j] = 1
      R[j, i] = 1
      eq = reshape(collect(R'), :)
      push!(eq, Y[i, j])
      push!(equations, eq)
    end
  end
  r = length(equations)
  l = length(equations[1])
  equations = elem_type(k)[i for i in Iterators.flatten(equations)]
  A = matrix(k,r,l, equations)
  c = A[:, l]
  A = A[:, 1:end-1]
  return A, c
end

@doc Markdown.doc"""
    _solve_X(Y::gfp_mat, b, g, ker=false) -> gfp_mat

Solve a certain linear equation modulo $2$. This is a helper function for
`_hensel_qf_modular_even`. We find the solution $X$ such that

$$Y = X + X^t$$

and

$$b_i = X_{ii} + \sum_{j=1}^n X_{ij}g_j \quad i \in \{1, \dots, n\}.$$
"""
function _solve_X(Y::gfp_mat, b, g)
  k = base_ring(Y)
  n = ncols(Y)
  # A*Xcoeff == c
  A, c = _solve_X_get_A_and_c(Y, b, g)
  fl, Xcoeff = can_solve_with_solution(A, c, side=:right)
  @assert fl
  X = matrix(k, n, n, reshape(collect(Xcoeff'), :))
  # confirm the computation
  @hassert :Lattice Y == X + X'
  for i in 1:n
    @hassert :Lattice b[i] == X[i,i] + sum([X[i,j]*g[j] for j in 1:n])
  end
  return X
end

function _solve_X_ker(Y::gfp_mat, b, g)
  # A*Xcoeff == c
  k = base_ring(Y)
  n = ncols(Y)
  A, c = _solve_X_get_A_and_c(Y, b, g)
  fl, Xcoeff = can_solve_with_solution(A, c, side=:right)
  if ker
    Ker = dense_matrix_type(k)[]
    r, K = right_kernel(A)
    for i in 1:r
      X = matrix(k, n, n, elem_type(k)[a for a in K[:,i]])
      push!(Ker, X)
    end
    return Ker
  end
end

@doc Markdown.doc"""
    hensel_qf(G::T, F::T, a, b, p) where T <: Union{nmod_mat, fmpz_mod_mat}

Lift `F` modulo `p^n` satisfying `G == F * G * F'`.

- `G` -- a block diagonal matrix of the form `[G0*p^n0, G1*p^n1,...,Gk*p^nk]`
  with integers `nk < .... < n1 < n0` and `Gi` unimodular and symmetric.
- `F` -- invertible `p`-adic matrix
  such that `(G, G, F)` is `a`-adapted
- `a` -- integer the starting precision
- `b`-- integer the target precision

Return `Fk` such that
- `Fk`- the lift of `F` such that
  `Z == F * G * F'` modulo `p^n` with `n = prec`
"""
function hensel_qf(G::T, F::T, a, b, p) where {T <: Union{nmod_mat, fmpz_mod_mat}}
  # Input checks
  @req isunit(det(F)) "F must be invertible"
  @req ncols(G)== ncols(F) && nrows(G) == nrows(F) "G, F must have the same size"
  @req base_ring(G) == base_ring(F) "not the same basering"
  @req issymmetric(G) "G must be symmetric"
  R = base_ring(G)
  #n = modulus(R)
  #@req(b > n,"Desired precision is higher than base ring precision")
  for k in 1:ncols(G)-1
    n1 = _min_val(G[k,:], p)
    n2 = _min_val(G[k+1,:], p)
    @req n1 >= n2 "block valuations must be descending"
  end
  @req _min_val(F*G*F'-G, p) >= a "F must satisfy G == F * G * F' mod p^a"
  if p == 2 & a == 1
    @req _min_val(diagonal(F*G*F'-G), p) >= a+1 "input is not adapted"
  end
  if ncols(F) == 0
    return F
  end
  # the real worker
  F = _hensel_qf(G, G, F, a, b, p) #works inplace
  return F
end

@doc Markdown.doc"""
    _block_indices_vals(G:::Union{nmod_mat, fmpz_mod_mat}, p)
                                                     -> Vector{Int}, Vector{Int}

Return a list of indices and a list of valuation of the homogeneous blocks.

The matrix `G` is assumed to be a symmetric `p`-adic block diagonal matrix with
modular blocks which have descending valuations.
"""
function _block_indices_vals(G::Union{nmod_mat, fmpz_mod_mat}, p)
  indices = []
  valuations = []
  while ncols(G) != 0
    i, val, _ = _last_block_index(G, p)
    push!(indices, i)
    push!(valuations, val)
    G = G[1:i-1,1:i-1]
  end
  reverse!(indices)
  reverse!(valuations)
  return indices, valuations
end

#=
    Helper function for `_hensel_qf`.

    Deals with the case that `G` is modular and `p=2`.

    INPUT:

    - `Z` -- symmetric `p`-adic `n \times n` matrix
    - `G` -- symmetric `p`-adic `n \times n` matrix
    - `F` -- symmetric `p`-adic `n \times n` matrix
    - `a` -- integer
    - `b` -- integer

    We require that the triple `(Z, G, F)` is `a`-adapted.

    OUTPUT:

    - `Fl` such that `(Z, G, Fl)` is `b`-adapted
    - raises a `ValueError` if `F` cannot be lifted
=#
function _hensel_qf_modular_even(Z::T, G::T, F::T, a, b) where {T <: Union{nmod_mat, fmpz_mod_mat}}
  n = ncols(Z)
  @req a != 0 "a must be a non-zero integer"
  if a == 1
    R = base_ring(Z)
    v = _min_val(G, 2)::Int
    G = divexact(G, 2^v)
    Z = divexact(Z, 2^v)
    Y = Z - F*G*F'
    X = _solve_X(divexact(Y, 2), [divexact(y, 4) for y in diagonal(Y)], diagonal(inv(G)))
    X = 2 * change_base_ring(R, lift(X))::T
    F = F + X*inv(G*F')
    a = 2
  end
  while a < b
    Y = Z - F*G*F'
    for i in 1:n
      #Y[i,i+1:end] = 0
      for j in i+1:ncols(Y)
        Y[i, j] = 0
      end
      Y[i,i] = divexact(Y[i,i], 2)
    end
    F = F + Y*inv(G*F')
    a = 2*a - 1
  end
  # confirm computation
  @hassert :Lattice _min_val(Z-F*G*F', 2) >= b
  @hassert :Lattice _min_val(diagonal(Z-F*G*F'),2) >= b + 1
  return F
end


text:
using QPDAS, Random, Test, LinearAlgebra
import QPDAS: deleterowcol!, addrowcol!, CholeskySpecial

Random.seed!(3141592)

m, n = 30, 100
i = [6,10,25,9,8]
# i = [6,8,9,10,25]
# i = [25,10,9,8,6]
A = randn(m,n)
M = A*A'
b = randn(m)
# Get smaller matrix
idx = BitArray(undef, m)
idx .= true
for j in i
    idx[j] = false
end
MS = M[idx, idx]
bs = b[idx]
# Get reference answer
as = MS\bs

# Test
F = CholeskySpecial(cholesky(M))
for j in i
    deleterowcol!(F, j)
end

a = F\b
norm(a[idx]-as)
@test a[idx] ≈ as
@test a[.!idx] ≈ b[.!idx]

# Test reverse, in some other order
i = [10,6,8,25,9]
for j in i
    addrowcol!(F, j)
end

@test cholesky(M)\b ≈ F.F\b
@test cholesky(M).U ≈ F.U

for j in randperm!(i)
    deleterowcol!(F, j)
end
for j in randperm!(i)
    addrowcol!(F, j)
end

@test cholesky(M)\b ≈ F.F\b
@test cholesky(M)\b ≈ F\b
@test cholesky(M).U ≈ F.U

@test size(F) == size(M)
@test size(F,1) == size(M,1)

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
module TestEnumerationArray

using FactCheck
using MultidimensionalTables

facts("EnumerationArray tests") do
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"])).elems --> Int[1 0 2 2;0 2 1 3]
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]).elems --> Int[2 0 3 3;0 3 2 1]
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]))[2,2].value --> "b"
  @fact EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"])[2,2].value --> "b"
  @fact tbltool.wrap_array(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]))[1:3]) --> tbltool.wrap_array(EnumerationArray(@nalift(["a",NA,NA])))
  @fact tbltool.wrap_array(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"])[1:3]) --> tbltool.wrap_array(EnumerationArray(@nalift(["a",NA,NA]), ["c","a","b"]))
  @fact size(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"])) --> (2, 4)
  @fact tbltool.wrap_array(reshape(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]), 4, 2)) --> tbltool.wrap_array(transpose(EnumerationArray(@nalift(["a" NA  NA "b";"b" "a" "b" "c"]), ["c","a","b"])))
  @fact tbltool.wrap_array(transpose(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]))) --> tbltool.wrap_array(EnumerationArray(@nalift(["a" NA;NA "b";"b" "a";"b" "c"])))
  @fact tbltool.wrap_array(transpose(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]))) --> tbltool.wrap_array(permutedims(EnumerationArray(@nalift(["a" NA "b" "b";NA "b" "a" "c"]), ["c","a","b"]), (2,1)))
  @fact tbltool.wrap_array(reverse(enumeration(["hello", "hello", "hi"]))) --> tbltool.wrap_array(enumeration(["hi", "hello", "hello"]))
end

end


text:
include("LatticeVars.jl")

function main()
    LatVecs = zeros(3,3)
    LatVecs[1,:] = [5.0, 5.0, 0.0]
    LatVecs[2,:] = [5.0, 0.0, 5.0]
    LatVecs[3,:] = [0.0, 6.0, 5.0]

    lattice_vars = LatticeVars( LatVecs )
    avec = lattice_vars.avec
    bvec = lattice_vars.bvec
    res = (avec*bvec')/(2*pi)
    display(res); println()

    println("Pass here")
end

main()


text:
@testset "quadraturepoints" begin
    for T in testtypes
        @test typeof(quadraturepoints(Triangle{T})) == QuadPts2D{T}
        @test typeof(quadraturepoints(Tetrahedron{T})) == QuadPts3D{T}
    end
end


text:
"""
Approximates p(query|evidence) with N weighted samples using likelihood
weighted sampling
"""
@with_kw struct LikelihoodWeightingInference <: InferenceMethod
    nsamples::Int = 500
end

"""
Approximates p(query|evidence) with `nsamples` likelihood weighted samples.

Since this uses a Factor, it is only efficient if the number of samples
is (signifcantly) greater than the number of possible instantiations for the
query variables
"""
function infer(im::LikelihoodWeightingInference, inf::InferenceState{BN}) where {BN<:DiscreteBayesNet}
    bn = inf.pgm
    nodes = names(bn)
    query = inf.query
    evidence = inf.evidence

    ϕ = Factor(query, map(n -> ncategories(bn, n), query))

    # if nodes are evidence
    evidence_mask = [haskey(evidence, s) for s in nodes]

    sample = Assignment()
    # add the evidence to the sample
    merge!(sample, evidence)

    # manual index into ϕ.potential since categorical implies Base.OneTo
    q_ind = similar(query, Int)

    for i in 1:im.nsamples
        w = 1.0

        for (cpd, is_ev) in zip(bn.cpds, evidence_mask)
            nn = name(cpd)
            d = cpd(sample)

            if !is_ev
                sample[nn] = rand(d)
            else
                w *= pdf(d, sample[nn])
            end
        end

        # pick out the query variables to index by
        for (i, q) in enumerate(query)
            q_ind[i] = sample[q]
        end

        ϕ.potential[q_ind...] += w
    end

    normalize!(ϕ)
    return ϕ
end


text:
using JuliaLambdaExample
using Test

@testset "JuliaLambdaExample.jl" begin
    result = JuliaLambdaExample.handle_event("1000000", String[])
    @test result isa String
    @test_nowarn parse(Float64, result)
    @test isapprox(parse(Float64, result), π, atol = 0.1)
end

max_stars_repo_name: ramboldio/aws-lambda-container-julia

max_issues_repo_name: ramboldio/aws-lambda-container-julia

max_forks_repo_name: ramboldio/aws-lambda-container-julia


text: my_f(x,y) = 2x+3y


text:
# Copyright Wen-Chih Chen
# the code is for computing basic DEA models based on orientation (input or output)
# and returns-to-scale (crs or vrs)

using JuMP
using Gurobi # Gurobi is used as the LP solver
using FrontierEfficiencyAnalysis # for large scale problems and/or with size-limited solver

function basicDEAmodel()
  #----------------------------------------
  # data section
  data = readcsv("example.csv") # input User's (.csv) data path
  # each row is for a DMU; each column is an input or output
  # in the example the first two columns are the inputs while the others are the outputs
  inputs = data[:, 1:2]
  outputs = data[:, 3:end]

  # ---------------------------------------
  # to determine the efficiencies of all DMUs
  results = zeros(size(inputs)[1])
  getAllEfficiencies!(results, inputs, outputs) # get all DMU's efficiencies of input-oriented crs model
  ## or use the following to set the options of orientation, returns-to-scale and method
  # getAllEfficiencies!(results, inputs, outputs, orientation = "input", rts = "crs", method = 0)
  ## orientation = "input" (default) or "output"
  ## rts = "crs" (default) or "vrs"
  ## method = 0 (default) or 1 (by FrontierEfficiencyAnalysis for large scale problems and/or with size-limited solver)
  ## e.g. for output-oriented vsr model and solve by FrontierEfficiencyAnalysis do:
  # getAllEfficiencies!(results, inputs, outputs, orientation = "output", rts = "vrs", method = 1)
  println(results)

  # ---------------------------------------
  # to determine the efficiency of DMU k
  k = 2
  score = getEfficiency(inputs, outputs, k) # get DMU k's efficiency of input-oriented crs model
  ## or use the following to set the options of orientation, returns-to-scale and method
  # score = getEfficiency(inputs, outputs, k, orientation = "input", rts = "crs", method = 0)
  ## orientation = "input" (default) or "output"
  ## rts = "crs" (default) or "vrs"
  ## method = 0 (default) or 1 (by FrontierEfficiencyAnalysis for large scale problems and/or with size-limited solver)
  ## e.g. for output-oriented vsr model and solve by FrontierEfficiencyAnalysis do:
  # score = getEfficiency(inputs, outputs, k, orientation = "output", rts = "vrs", method = 1)
  println(score)

end

# compute DEA efficiencies for all DMUs in the data set and store results in scores
function getAllEfficiencies!(scores, inputData, outputData; orientation = "input", rts = "crs", method = 0)
  # parameters
  # scores: efficiency results
  # inputData: input data matrix (DMU size) x (input size)
  # outputData: output data matrix (DMU size) x (output size)
  # orientation: input- or output- oriented model (value= "input" or "output")
  # rts: returns to scale setting (value = "crs" or "vrs")
  # # method: how to solve LP (value = 0 (default) or 1 (by FrontierEfficiencyAnalysis.jl))

  if !(orientation == "input" || orientation == "output")
    error("orientation not supported")
  end

  if !(rts == "crs" || rts == "vrs")
    error("returns to scale not supported")
  end

  for t = 1:size(inputData)[1]
    scores[t] = getEfficiency(inputData, outputData, t, orientation = orientation, rts = rts, method = method)
  end
end

# compute DEA efficiency for DMU k
function getEfficiency(inputs, outputs, k; orientation = "input", rts = "crs", method = 0)
  # parameters
  # inputs: input data matrix (DMU size) x (input size)
  # outputs: output data matrix (DMU size) x (output size)
  # k: index of DMU to be evaluated
  # orientation: input- or output- oriented model (value= "input" or "output")
  # rts: returns to scale setting (value = "crs" or "vrs")
  # method: how to solve LP (value = 0 (default) or 1 (by FrontierEfficiencyAnalysis.jl))
  scale, inputNo = size(inputs)
  outputNo = size(outputs)[2]

  model = Model(solver = GurobiSolver(OutputFlag=0)) # Gurobi is used as the LP solver here. Users can choose their favorite solver.
  @variable(model, Lambda[1:scale] >= 0)
  @variable(model, Theta)

  if orientation == "input" # input-oriented model
    @objective(model, Min, Theta)
    @constraint(model, inputCon[i=1:inputNo], sum(Lambda[r]*inputs[r,i] for r = 1:scale) <= Theta*inputs[k,i])
    @constraint(model, outputCon[j=1:outputNo], sum(Lambda[r]*outputs[r,j] for r = 1:scale) >= outputs[k,j])
  elseif orientation == "output" # output-oriented model
    @objective(model, Max, Theta)
    @constraint(model, inputCon[i=1:inputNo], sum(Lambda[r]*inputs[r,i] for r = 1:scale) <= inputs[k,i])
    @constraint(model, outputCon[j=1:outputNo], sum(Lambda[r]*outputs[r,j] for r = 1:scale) >= Theta*outputs[k,j])
  else
    error("orientation not support")
  end

  if rts == "vrs"
    @constraint(model, sum(Lambda[r] for r = 1:scale) == 1)
  elseif rts == "crs"
  else
    error("returns to scale not supported")
  end

  if method == 0 # defalut
    solve(model)
  else  # using FrontierEfficiencyAnalysis for large scale problem or with size limited solver
    solveDEA(model)
  end

  return getvalue(Theta)
end

basicDEAmodel()

max_stars_repo_name: wen-chih/code-for-DEA


text:
using Images
function voronoi(w, h, n_centroids)
    dist = (point,vector) -> sqrt.((point[1].-vector[:,1]).^2 .+ (point[2].-vector[:,2]).^2)
    dots = [rand(1:h, n_centroids) rand(1:w, n_centroids) rand(RGB{N0f8}, n_centroids)]
    img  = zeros(RGB{N0f8}, h, w)
    for x in 1:h, y in 1:w
        distances = dist([x,y],dots) # distance
        nn = findmin(distances)[2]
        img[x,y]  = dots[nn,:][3]
    end
    return img
end
img = voronoi(800, 600, 200)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function KS_solve_scf( LF::LF3dGrid, Gv::GvectorsT,
                      ∇2, precH,
                      V_ionic, Focc, Ncols::Int64;
                      rho0 = nothing,
                      Potentials0 = nothing,
                      α_t=3e-5, Niter=1000,
                      verbose=false )
  #
  Npoints = LF.Nx * LF.Ny * LF.Nz
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h

  # starting rho
  if rho0 == nothing
    srand(1234)
    v = rand( Npoints, Ncols )
    v = orthonormalize( LF, v )
    rho = calc_rho( Focc, v )
  else
    srand(1234)
    v = rand( Npoints, Ncols ) # for initial guess in
    rho = copy( rho )
  end

  # starting Potentials
  if Potentials0 == nothing
    # initialize potential: Ionic, Hartree, XC
    V_Hartree = solve_poisson_FFT( Gv, rho )
    V_xc = excVWN( rho ) + rho .* excpVWN( rho )
    Potentials = PotentialsT( V_ionic, V_Hartree, V_xc )
  else
    Potentials = copy(Potentials0)
  end

  Etot_old = 0.0
  rho_old = copy(rho)
  Energies = EnergiesT(0.0, 0.0, 0.0, 0.0, 0.0)
  λ = zeros(Ncols)
  for iter = 1:50
    λ, v = diag_lobpcg( LF, ∇2, precH, Potentials, v, verbose_last=true )
    v = v[:,:]/sqrt(ΔV)

    # Calculate energies and update potetials
    Energies = calc_Energies( LF, ∇2, Potentials, rho, Focc, v )
    #
    Etot = Energies.Total
    diffE = abs(Etot-Etot_old)
    Etot_old = Etot
    print_Energies( Energies )
    if diffE < 1e-7
      break
    end
    #
    rho_new = calc_rho( Focc, v )
    diffRho = norm(rho_new-rho)
    @printf("%8d %18.10f %18.10e %18.10e\n", iter, Etot, diffE, diffRho )
    rho = 0.5*rho_new[:] + 0.5*rho[:]
    integRho = sum(rho)*ΔV
    @printf("integRho = %18.10f\n", integRho)
    #
    V_Hartree = solve_poisson_FFT( Gv, rho )
    V_xc = excVWN( rho ) + rho .* excpVWN( rho )
    Potentials = PotentialsT( V_ionic, V_Hartree, V_xc )
  end

  return Energies, v, Potentials
end


text:
###### Operations for JuMPExpressionMatrix ######

function _remove_undef!(expression_array::T) where T <: JuMPExpressionMatrix

    for j in 1:size(expression_array)[2]

        for i in 1:size(expression_array)[1]

            if !isassigned(expression_array, i, j)
                 expression_array[i,j] = zero(eltype(expression_array))
            else
                continue
            end

        end

    end

    return

end

###### Operations for JuMPDenseAxisExpressionMatrix ######

function _remove_undef!(axis_expression_array::JuMP.Containers.DenseAxisArray)
    _remove_undef!(axis_expression_array.data)
end

function _add_to_expression!(expression_array::T,
                             ix::Int64,
                             jx::Int64,
                             var::JV,
                             sign::Int64) where {T, JV <: JuMP.AbstractVariableRef}

    if isassigned(expression_array,  ix, jx)
        JuMP.add_to_expression!(expression_array[ix,jx], sign, var)
    else
        expression_array[ix,jx] = sign*var
    end

    return

end

function _add_to_expression!(expression_array::T,
                             ix::Int64,
                             jx::Int64,
                             var::JV) where {T, JV <: JuMP.AbstractVariableRef}

    if isassigned(expression_array, ix, jx)
        JuMP.add_to_expression!(expression_array[ix,jx], 1.0, var)
    else
        expression_array[ix,jx] = 1.0*var
    end

    return

end

function _add_to_expression!(expression_array::T,
                             ix::Int64,
                             jx::Int64,
                             value::Float64) where T

    if isassigned(expression_array, ix, jx)
        expression_array[ix,jx].constant +=  value
    else
        expression_array[ix,jx] = zero(eltype(expression_array)) + value
    end

    return

end

function _add_to_expression!(expression_array::T,
                            ix::Int64,
                            jx::Int64,
                            parameter::PJ.ParameterRef) where T

    if isassigned(expression_array, ix, jx)
        JuMP.add_to_expression!(expression_array[ix,jx], 1.0, parameter);
    else
        expression_array[ix,jx] = zero(eltype(expression_array)) + parameter;
    end

    return

end


text:
using Random

function _benchmark_opt(tfopt, adcmeopt)
    reset_default_graph()
    n = 100
    Random.seed!(233)
    θ = Variable(ae_init([2,10,10,2]))
    y = fc(rand(10,2), [10,10,2], θ)
    loss = sum(y^2)
    opt = tfopt().minimize(loss)

    global sess = Session(); init(sess)
    qs = zeros(10)
    for i = 1:10
        qs[i] = run(sess, loss)
        run(sess, opt)
    end

# ------------------ ADCME ----------------------
    init(sess)
    opt = adcmeopt()
    uo = UnconstrainedOptimizer(sess, loss, vars=[θ])
    x0 = getInit(uo)
    fs = zeros(10)
    for i = 1:10
        fs[i], df = getLossAndGrad(uo, x0)
        Δ = apply!(opt, x0, df)
        x0 -= Δ
    end
    flag =  maximum(abs.(fs-qs)./abs.(qs))<1e-2
    @info maximum(abs.(fs-qs)./abs.(qs))
    println(fs)
    println(qs)
    return flag
end

function _print_opt(adcmeopt)
    reset_default_graph()
    n = 100
    Random.seed!(233)
    θ = Variable(ae_init([2,10,10,2]))
    y = fc(rand(10,2), [10,10,2], θ)
    loss = sum(y^2)
    global sess = Session(); init(sess)
    opt = adcmeopt()
    uo = UnconstrainedOptimizer(sess, loss, vars=[θ])
    x0 = getInit(uo)
    fs = zeros(10)
    for i = 1:10
        fs[i], df = getLossAndGrad(uo, x0)
        Δ = apply!(opt, x0, df)
        x0 -= Δ
    end
    println("$adcmeopt: ")
    println(fs)
end

@testset "Optimizers" begin
    @test _benchmark_opt(GradientDescentOptimizer, Descent)
    @test _benchmark_opt(AdamOptimizer, ADAM)
    @test _benchmark_opt(AdagradOptimizer, ADAGrad)
    @test _benchmark_opt(AdadeltaOptimizer, ADADelta)
    _print_opt(RMSProp)
    _print_opt(AMSGrad)
    _print_opt(NADAM)
    _print_opt(Momentum)
    _print_opt(Nesterov)
    _print_opt(RADAM)
    _print_opt(AdaMax)
end


text:
module VertexEliminationOrders

using Base.Threads
using Random

# Sub-modules to support dfs and heuristics.
include("mmdqueues.jl"); include("graphs.jl")
using .MMDQueues, .Graphs
export Graph

# Utility functions for loading and analysing graphs.
include("utils.jl")

include("treewidth_heuristics.jl")
include("dfs_utils.jl")
include("dfs.jl")

end


text: #ITensors.hasind(s::Index) = x -> hasind(x, s)


text:
####### STATE SPACE ##############################

function POMDPs.states(pomdp::SingleOCPOMDP)
    env = pomdp.env
    rl = env.params.roadway_length
    cw = env.params.crosswalk_width
    lw = env.params.lane_width
    cl = env.params.crosswalk_length
    x_ped = 0.5*(env.params.roadway_length - env.params.crosswalk_width + 1)
    Y = LinRange(-cl/4, cl/4, Int(floor(cl/2/pomdp.pos_res)) + 1)
    V_ped = LinRange(0, env.params.ped_max_speed, Int(floor(env.params.ped_max_speed/pomdp.vel_res)) + 1)
    X = LinRange(5., rl/2 + 2*cw , Int(floor((rl/2 + 2*cw -5)/pomdp.pos_res)) + 1)
    V = LinRange(0., env.params.speed_limit, Int(floor(env.params.speed_limit/pomdp.vel_res)) + 1)
    y_ego = 0. #XXX might need to be a parameter
    space = SingleOCState[]
    for x in X
        for v in V
            for y in Y
                for v_ped in V_ped
                    ego = VehicleState(VecSE2(x, y_ego, 0.), env.roadway, v)
                    ped = VehicleState(VecSE2(x_ped, y, pi/2), env.roadway, v_ped)
                    push!(space, SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped))
                end
            end
        end
    end
    for x in X
        for v in V
            ego = VehicleState(VecSE2(x, y_ego, 0.), env.roadway, v)
            ped =  get_off_the_grid(pomdp)
            push!(space, SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped))
        end
    end
    push!(space, TERMINAL_STATE)
    return space
end

function POMDPs.stateindex(pomdp::SingleOCPOMDP, s::SingleOCState)
    if isterminal(pomdp, s)
        return n_states(pomdp)
    end
    env = pomdp.env
    rl = env.params.roadway_length
    cw = env.params.crosswalk_width
    lw = env.params.lane_width
    cl = env.params.crosswalk_length
    size_x = Int(floor((rl/2 + 2*cw -5)/pomdp.pos_res) + 1)
    size_v = Int(floor(env.params.speed_limit/pomdp.vel_res) + 1)
    size_y = Int(floor(cl/2/pomdp.pos_res) + 1)
    size_v_ped = Int(floor(env.params.ped_max_speed/pomdp.vel_res) + 1)

    x = s.ego.posG.x
    v = s.ego.v
    y = s.ped.posG.y
    v_ped = s.ped.v

    y_0 = -env.params.crosswalk_length/4
    x_0 = 5
    x_ind = ceil(Int, (x - x_0)/pomdp.pos_res) + 1
    v_ind = ceil(Int, v/pomdp.vel_res) + 1
    y_ind = ceil(Int, (y - y_0)/pomdp.pos_res) + 1

    if !off_the_grid(pomdp, s.ped)
        v_ped_ind = ceil(Int, v_ped/pomdp.vel_res) + 1
        i = LinearIndices((size_v_ped, size_y, size_v, size_x))[v_ped_ind, y_ind, v_ind, x_ind]
    else
        off_grid_ind = LinearIndices((size_v, size_x))[v_ind, x_ind]
        i = size_x*size_v*size_y*size_v_ped + off_grid_ind
    end
    return i
end

function n_states(pomdp::SingleOCPOMDP)
    env = pomdp.env
    rl = env.params.roadway_length
    cw = env.params.crosswalk_width
    lw = env.params.lane_width
    cl = env.params.crosswalk_length
    size_x = floor(Int, (rl/2 + 2*cw -5)/pomdp.pos_res) + 1
    size_v = floor(Int, env.params.speed_limit/pomdp.vel_res) + 1
    size_y = floor(Int, cl/2/pomdp.pos_res) + 1
    size_v_ped = floor(Int, env.params.ped_max_speed/pomdp.vel_res) + 1
    return size_x*size_v*(size_y*size_v_ped + 1) + 1
end

### OBSERVATION SPACE

function POMDPs.observations(pomdp::SingleOCPOMDP)
    return states(pomdp)
end

function POMDPs.obsindex(pomdp::SingleOCPOMDP, o::SingleOCObs)
    return stateindex(pomdp, o)
end

function n_observations(pomdp::SingleOCPOMDP)
    return n_states(pomdp)
end

#### ACTION SPACE
function POMDPs.actions(pomdp::SingleOCPOMDP)
    return [SingleOCAction(-2*pomdp.max_acc), SingleOCAction(-pomdp.max_acc), SingleOCAction(0.0), SingleOCAction(pomdp.max_acc)]
end

function POMDPs.actionindex(pomdp::SingleOCPOMDP, a::SingleOCAction)
    if a.acc < -2.
        return 1
    elseif -2. <= a.acc < 0.
        return 2
    elseif a.acc == 0.
        return 3
    else
        return 4
    end
end

######################### HELPERS ##################################################################

"""
    off_the_grid(pomdp::SingleOCPOMDP, ped::VehicleState)
Check if the current state of the pedestrian is in the grid
"""
function off_the_grid(pomdp::SingleOCPOMDP, ped::VehicleState)
    env = pomdp.env
    obstacle = env.obstacles[1]
    obs_center = get_center(obstacle)
    return ped.posG ≈ obs_center
end

"""
    get_off_the_grid(pomdp::SingleOCPOMDP)
return the off the grid pedestrian state
"""
function get_off_the_grid(pomdp::SingleOCPOMDP)
    env = pomdp.env
    obstacle = env.obstacles[1]
    obs_center = get_center(obstacle)
    return VehicleState(VecSE2(obs_center.x, obs_center.y, 0.), pomdp.env.roadway, 0.0)
end

"""
    get_X_grid(pomdp::SingleOCPOMDP)
Helper function to return the X discretization
"""
function get_X_grid(pomdp::SingleOCPOMDP)
    env = pomdp.env
    rl = env.params.roadway_length
    cw = env.params.crosswalk_width
    lw = env.params.lane_width
    cl = env.params.crosswalk_length
    return LinRange(5., rl/2 + 2*cw , Int(floor((rl/2 + 2*cw -5)/pomdp.pos_res)) + 1)
end

"""
    get_V_grid(pomdp::SingleOCPOMDP)
Helper function to return the ego velSingleOCity discretization
"""
function get_V_grid(pomdp::SingleOCPOMDP)
    env = pomdp.env
    V = LinRange(0., env.params.speed_limit, Int(floor(env.params.speed_limit/pomdp.vel_res)) + 1)
end

"""
    get_Y_grid(pomdp::SingleOCPOMDP)
Helper function to return the Y discretization
"""
function get_Y_grid(pomdp::SingleOCPOMDP)
    cl = pomdp.env.params.crosswalk_length
    return LinRange(-cl/4, cl/4, Int(floor(cl/2/pomdp.pos_res)) + 1)
end

"""
    get_V_ped_grid(pomdp::SingleOCPOMDP)
Helper funcion to return the VelSingleOCity discretization for the pedestrian
"""
function get_V_ped_grid(pomdp::SingleOCPOMDP)
    env = pomdp.env
    V_ped = LinRange(0, env.params.ped_max_speed, Int(floor(env.params.ped_max_speed/pomdp.vel_res)) + 1)
end

"""
    get_y_index(pomdp::SingleOCPOMDP, y::Float64)
return the index of y in the Y grid
"""
function get_y_index(pomdp::SingleOCPOMDP, y::Float64)
    y_0 = -pomdp.env.params.crosswalk_length/4
    y_ind = Int(ceil((y - y_0)/pomdp.pos_res)) + 1
    return y_ind
end

"""
    get_v_ped_index(pomdp::SingleOCPOMDP, v_ped::Float64)
return the index of v_ped in the V_ped grid
"""
function get_v_ped_index(pomdp::SingleOCPOMDP, v_ped::Float64)
    v_ped_ind = Int(ceil(v_ped/pomdp.vel_res)) + 1
    return v_ped_ind
end


text:
#
# Computes the electrostatic interaction of a pair of atoms, given the distance
#

function qpair(d,q1,q2)

  Coulomb = 332.05382e0
  qpair = Coulomb*q1*q2 / d

  return qpair
end


text:
const fn_primitives = [:sample, :observe, :transform, :plate, :track]

@doc raw"""
    track(t, address, value)

Tracks the `value` in the trace `t` at `address`. Does not affect the log
probability of the trace.
"""
track(t, address, value) = transform(t, address, x -> x, value)
export track

get_fname_from_expr(f) = f.args[1].args[1]  # function name, the function has to be defined function ... end

function expand_trace_statements(__ex, primitives)
    if @capture(__ex, function_name_(xs__))
        function_name in primitives ? :($function_name(t, $(xs...))) : __ex
    else
        __ex
    end
end

function expand_statements(__ex, primitives)
    if @capture(__ex, varname_ = function_name_(xs__))
        if function_name in primitives
            if function_name != :plate
                :($varname = $function_name(t, $(string(varname)), $(xs...)))
            else
                # plate(t, op, address, ...) <- plate(op,...)
                :($varname = $function_name(t, $(xs[1]), $(string(varname)), $(xs[2:end]...)))
            end
        else
            __ex
        end
    else
        __ex
    end
end

@doc raw"""
    macro model(f)

Create a valid `CrimsonSkyline` model from a dsl that omits passthrough of the `Trace` data structure
and obviates address creation. For example, the model
```
function f(t::Trace, data::Vector{Float64})
    loc = sample(t, "loc", Normal())
    log_scale = sample(t, "log_scale", Normal())
    scale = transform(t, "scale", x -> exp(x), log_scale)
    plate(t, observe, "data", Normal(loc, scale), data)
end
```
can be created instead by writing
```
@model function simple_model(data::Vector{Float64})
    loc = sample(Normal())
    log_scale = sample(Normal())
    scale = track(exp(log_scale))
    data = plate(observe, Normal(loc, scale), data)
end
```
The resulting function `f` is exactly equal to the first definition.
The LHS of assignments *must* be parsed as `Symbols`, e.g., `my_variable = ...` is always safe;
`Expr`s like `out[i] = ...` won't _necessarily_ parse correctly.
The address of each site is equal to the string representation of the LHS of the assignment.
For example, `log_scale = sample(Normal())` means that there will be a site called `"log_scale"`
in the trace.
"""
macro model(f)
    fname = get_fname_from_expr(f)
    f = prewalk(x -> expand_statements(x, fn_primitives), f)
    f = prewalk(x -> expand_trace_statements(x, [fname]), f)
    return quote
        $(esc(fname)) = $f
    end
end
export @model


text:
###
# Methods to transform of heta-compiler generated files into Platform content
###

const SUPPORTED_VERSIONS = ["0.6.5", "0.6.6"]

# transformation of tuple to Platform
function Platform(
    models::NamedTuple,
    scenarios::Tuple,
    version::String
)
    # TODO: semver approach might be better
    @assert length(indexin(version, SUPPORTED_VERSIONS)) > 0 "Heta compiler of the version \"$version\" is not supported."

    print("Loading platform... ")
    model_pairs = [pair[1] => Model(pair[2]...) for pair in pairs(models)]

    platform = Platform(
        OrderedDict{Symbol,Model}(model_pairs),
        OrderedDict{Symbol,Scenario}()
    )
    println("OK!")

    return platform
end

# transformation of tuple to Model
function Model(
    init_func::Function,
    ode_func::Function,
    time_events::NamedTuple,
    c_events::NamedTuple,
    stop_events::NamedTuple,
    saving_generator::Function,
    constants_num::NamedTuple,
    events_active::NamedTuple,
    records_output::NamedTuple
)
    events = Pair[]
    ## FIXME : remove event name from heta-compiler
    for (name,event_tuple) in pairs(time_events) # time events
        evt = name => TimeEvent(event_tuple[1], event_tuple[2], event_tuple[4])
        push!(events, evt)
    end
    for (name,event_tuple) in pairs(c_events) # c events
        evt = name => CEvent(event_tuple[1], event_tuple[2], event_tuple[4])
        push!(events, evt)
    end
    for (name,event_tuple) in pairs(stop_events) # stop events
        evt = name => StopEvent(event_tuple[1], event_tuple[3])
        push!(events, evt)
    end

    #observable_pairs = filter((p) -> p[2], pairs(records_output)) # from records_output
    #observables = Symbol[p[1] for p in observable_pairs]
    records_output_ = collect(Pair{Symbol,Bool}, pairs(records_output))
    events_active_ = collect(Pair{Symbol,Bool}, pairs(events_active))

    ### fake run, disabled because it slows model loading down
    #=
    _u0, _p0 = init_func(constants_num)
    constants = LVector(constants_num)
    _params = Params(constants, _p0)
    prob = ODEProblem(ode_func, _u0, (0.,1.), _params)

    # check if default alg can solve the prob
    integrator = init(prob, DEFAULT_ALG)
    step!(integrator)
    ret = check_error(integrator)
    ret != :Success && @warn "Default algorithm returned $ret status. Consider using a different algorithm."
    =#
    model = Model(
        init_func,            # init_func
        ode_func,             # ode_func
        NamedTuple(events),   # events :: Changed to NamedTuple
        saving_generator,     # saving_generator
        records_output_,
        constants_num,        # constants :: Changed to NamedTuple
        NamedTuple(events_active_) # events_active :: Changed to NamedTuple
    )

    return model
end


text:
function randpoly(R::Union{ FqNmodPolyRing, FqPolyRing }, dmax::Int)
  r = R()
  F = base_ring(R)
  d = rand(0:dmax)
  for i = 0:d
    setcoeff!(r, i, rand(F))
  end
  return r
end

@testset "Extensions of finite fields" begin
  Fp = GF(2)
  Fpx, x = Fp["x"]
  f = x^3 + x + 1
  Fq = FqNmodFiniteField(f, :$)
  a = gen(Fq)
  Fqy, y = Fq["y"]
  g = y^3 + a*y + a^2
  G, FqytoG = Hecke.field_extension(g)
  @test characteristic(G) == 2
  @test degree(G) == 9
  GtoFqy = pseudo_inv(FqytoG)

  # Check if FqytoG does something reasonable with constant polynomials
  @test iszero(FqytoG(Fqy(1))^8 - FqytoG(Fqy(1)))
  @test iszero(FqytoG(Fqy(a))^8 - FqytoG(Fqy(a)))
  @test iszero(FqytoG(Fqy(a^2))^8 - FqytoG(Fqy(a^2)))

  # Check if FqytoG is linear
  for i = 1:5
    h1 = randpoly(Fqy, 5)
    h2 = randpoly(Fqy, 5)
    @test mod(h1, g) == GtoFqy(FqytoG(h1))
    @test mod(h2, g) == GtoFqy(FqytoG(h2))
    @test FqytoG(h1 + h2) == FqytoG(h1) + FqytoG(h2)
    @test FqytoG(h1*h2) == FqytoG(h1)*FqytoG(h2)
  end
  @test iszero(FqytoG(Fqy()))

  # Now everything with FqFiniteFields
  Fp = GF(fmpz(2))
  Fpx, x = Fp["x"]
  f = x^3 + x + 1
  Fq = FqFiniteField(f, :$)
  a = gen(Fq)
  Fqy, y = Fq["y"]
  g = y^3 + a*y + a^2
  G, FqytoG = Hecke.field_extension(g)
  @test characteristic(G) == 2
  @test degree(G) == 9
  GtoFqy = pseudo_inv(FqytoG)

  # Check if FqytoG does something reasonable with constant polynomials
  @test iszero(FqytoG(Fqy(1))^8 - FqytoG(Fqy(1)))
  @test iszero(FqytoG(Fqy(a))^8 - FqytoG(Fqy(a)))
  @test iszero(FqytoG(Fqy(a^2))^8 - FqytoG(Fqy(a^2)))

  # Check if FqytoG is linear
  for i = 1:5
    h1 = randpoly(Fqy, 5)
    h2 = randpoly(Fqy, 5)
    @test mod(h1, g) == GtoFqy(FqytoG(h1))
    @test mod(h2, g) == GtoFqy(FqytoG(h2))
    @test FqytoG(h1 + h2) == FqytoG(h1) + FqytoG(h2)
    @test FqytoG(h1*h2) == FqytoG(h1)*FqytoG(h2)
    @test iszero(FqytoG(Fqy()))
  end
end


text:
using LearnBenchmark

d = 5
N = 5000
bound_ratio = .01

offset = zeros(d,N)
offset[1,:] .+= 1

dataset = [randn(d,N), randn(d,N) + offset]

L = d+1

@info "Ensemble Learner Results"
@info bayeserror(dataset,L,bound_ratio)
#println()
@info "Base Learner Results"
@info baselearner(dataset,3, bound_ratio)


text:
using Mmap
using Printf

# Test by reading it back in
s = open("zset")   # default is read-only
n = read(s, Int)
println("size of zset is $n")
zs = Mmap.mmap(s, Array{ComplexF64,1}, n)
close(s)
f =  zs[1]
l = zs[n]
@printf("First = %f, %f  last = %f, %f\n", f.re, f.im, l.re, l.im)

max_stars_repo_name: rabarar/juliaSnippets

max_issues_repo_name: rabarar/juliaSnippets

max_forks_repo_name: rabarar/juliaSnippets


text:
#=
Test Environment
 - Julia : v1.3.1
 - Flux  : v0.10.1
=#

# Classifies MNIST digits with a convolution network.
# Writes out saved model to the file "mnist_conv.bson".
# Demonstrates basic model construction, training, saving,
# conditional early-exits, and learning rate scheduling.
#
# This model, while simple, should hit around 99% test
# accuracy after training for approximately 20 epochs.

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition
using Printf, BSON
using Dates
using CUDAnative: device!
using CuArrays
using Random
using Dates

model_file = joinpath(dirname(@__FILE__),"conv_gpu_minibatch.bson")

epochs = 100
batch_size = 128
gpu_device = 0

# set using GPU device
device!(gpu_device)
CuArrays.allowscalar(false)

# Bundle images together with labels and groups into minibatch
function make_minibatch(imgs,labels,batch_size)
    len = length(imgs)
    sz = size(imgs[1])
    data_set =
    [(cat([reshape(Float32.(imgs[i]),sz...,1,1) for i in idx]...,dims=4),
      float.(onehotbatch(labels[idx],0:9)) ) for idx in partition(1:len,batch_size) ]
    return data_set
end

# Train data load
train_labels = MNIST.labels()
train_imgs = MNIST.images()
# Make train data to minibatch
train_set = make_minibatch(train_imgs,train_labels,batch_size)

# Test data load
test_labels = MNIST.labels(:test)
test_imgs = MNIST.images(:test)
test_set = make_minibatch(test_imgs,test_labels,batch_size)

#=
 Define our model. We will use a simple convolutional architecture with
 three iterations of Conv -> ReLu -> MaxPool, followed by a final Dense
 layer that feeds into a softmax probability output.
=#
@info "Construncting model..."
model = Chain(
  # First convolution, operating upon a 28x28 image
  Conv((3,3), 1=>16, pad=(1,1), relu),
  MaxPool((2,2)),

  # Second convolution, operating upon a 14x14 image
  Conv((3,3), 16=>32, pad=(1,1), relu),
  MaxPool((2,2)),

  # Third convolution, operating upon a 7x7 image
  Conv((3,3), 32=>32, pad=(1,1), relu),
  MaxPool((2,2)),

  # Reshape 3d tensor into a 2d one, at this point it should be (3,3,32,N)
  # which is where we get the 288 in the `Dense` layer below:
  x -> reshape(x, :, size(x,4)),
  Dense(288,10),

  # Finally, softmax to get nice probabilities
  softmax,
)

m = model |> gpu

#=
`loss()` calculates the crossentropy loss between our prediction `ŷ`
 (calculated from `m(x)`) and the ground truth `y`. We augment the data
 a bit, adding gaussian random noise to our image to make it more robust.
 =#
function loss(x,y)
  ŷ = m(x)
  return crossentropy(ŷ,y)
end

function accuracy(data_set)
  batch_size = size(data_set[1][1])[end]
  l = length(data_set)*batch_size
  s = 0f0
  for (x,y) in data_set
    s += sum((onecold(m(x|>gpu) |> cpu) .== onecold(y|>cpu)))
  end
  return s/l
end

# Make sure our is nicely precompiled befor starting our training loop
m(train_set[1][1] |> gpu)

# Train our model with the given training set using the ADAM optimizer and
# printing out performance aganin the test set as we go.
opt = ADAM(0.001)

@info "Beginning training loop..."
best_acc = 0.0
last_improvement = 0

for epoch_idx in 1:epochs
  global best_acc, last_improvement
  suffle_idxs = collect(1:length(train_set))
  shuffle!(suffle_idxs)
  for idx in suffle_idxs
    (x,y) = train_set[idx]
    # We augment `x` a little bit here, adding in random noise
    x = (x .+ 0.1f0*randn(eltype(x),size(x))) |> gpu
    y = y|> gpu
    Flux.train!(loss, params(m), [(x, y)],opt)
  end
  acc = accuracy(test_set)
  @info(@sprintf("[%d]: Test accuracy: %.4f",epoch_idx,acc))

  # If our accuracy is good enough, quit out.
  if acc >= 0.999
    @info " -> Early-exiting: We reached our target accuracy of 99.9%"
    break
  end

  # If this is the best accuracy we've seen so far, save the model out
  if acc >= best_acc
    @info " -> New best accuracy! saving model out to $(model_file)"
    model = m |> cpu
    acc = acc |> cpu
    BSON.@save model_file model epoch_idx acc
    best_acc = acc
    last_improvement = epoch_idx
  end

  #If we haven't seen improvement in 5 epochs, drop out learing rate:
  if epoch_idx - last_improvement >= 5 && opt.eta > 1e-6
    opt.eta /= 10.0
    @warn " -> Haven't improved in a while, dropping learning rate to $(opt.eta)!"

    # After dropping learing rate, give it a few epochs to improve
    last_improvement = epoch_idx
  end

  if epoch_idx - last_improvement >= 10
    @warn " -> We're calling this converged."
    break
  end
end

max_stars_repo_name: mrchaos/model-zoo

max_issues_repo_name: mrchaos/model-zoo

max_forks_repo_name: mrchaos/model-zoo


text:
module Twiddle

"""
    repeatbyte{T<:Unsigned}(::Type{T}, byte::UInt8)

Repeats the bitpattern of the `byte` throughout the an Unsigned integer of
type `T`.

This is useful for bit-twiddling code that is to work on any word size.

Let's use a very simple example to show why this is useful.
Let's say you have a function in which a byte is masked with a bitmask
0x33.

```julia
f(x::UInt8) = x & 0x33
```

Suppose you wanted this to work for any word size, so 8 bit, 16 bit, 32 bit, 64
bit, and 128 bit words? One way is to write `f` as a set of seperate methods:

```julia
f(x::UInt8) = x & 0x33
f(x::UInt16) = x & 0x3333
f(x::UInt32) = x & 0x33333333
f(x::UInt64) = x & 0x3333333333333333
f(x::UInt128) = x & 0x33333333333333333333333333333333
```

This is clearly not ideal, you would be writing out the same method multiplet
times, but with a differently sized literal for the mask each time.

Alternatively, you could write one parametric function with repeatbyte:

```julia
f2{T<:Unsigned}(x::T) = x & Twiddle.repeatbyte(T, 0x33)
```

You might expect this to be less efficient - `repeatbyte` uses several operations
to generate the values 0x33, 0x3333 and so on, whereas in the seperate methods,
those literal values are hard coded. However, thanks to constant folding during
compilation, those operations are done once at compilation time and so the
native instructions generated are identical.
"""
@inline function repeatbyte(::Type{T}, byte::UInt8) where {T<:Unsigned}
    return div(typemax(T), 0xff) * byte
end

"""
    mask{T<:Unsigned}(::Type{T}, n::Integer)

Creates a bit mask for given number of bits `n`.

The mask starts from the least significant bit, and end at bit `n`.

e.g:

```jldoctest
julia> Twiddle.mask(UInt64, 8)
0x00000000000000ff
```
"""
@inline mask(::Type{T}, n::Integer) where {T<:Unsigned} = (T(1) << n) - 0x1
@inline mask(n::Integer) = mask(UInt64, n)

"""
    swapbits{T<:Unsigned}(x::T, i::Integer, j::Integer)

Swap the i'th and j'th bits in an unsigned integer.
Note this uses zero based indexes for `i` and `j`.

E.g. to swap the LSB and MSB of a byte: 1001 1000 (0x98) -> 0001 1001 (0x19)

```@example
swapbits(0x98, 0, 7)
```
"""
@inline function swapbits(x::T, i::Integer, j::Integer) where {T<:Unsigned}
    ibit = (x >> i) & T(1)
    jbit = (x >> j) & T(1)
    ixj = ibit ⊻ jbit
    ixj = (ixj << i) | (ixj << j)
    return x ⊻ ixj
end

include("nibbles.jl")
include("bitpairs.jl")

end


text:
"""
Class for variational Gaussian Processes models (non-sparse)

```julia
MCGP(X::AbstractArray{T1,N1},y::AbstractArray{T2,N2},kernel::Union{Kernel,AbstractVector{<:Kernel}},
    likelihood::LikelihoodType,inference::InferenceType;
    verbose::Int=0,optimizer::Union{Bool,Optimizer,Nothing}=Adam(α=0.01),atfrequency::Integer=1,
    mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(),
    IndependentPriors::Bool=true,ArrayType::UnionAll=Vector)
```

Argument list :

**Mandatory arguments**

 - `X` : input features, should be a matrix N×D where N is the number of observation and D the number of dimension
 - `y` : input labels, can be either a vector of labels for multiclass and single output or a matrix for multi-outputs (note that only one likelihood can be applied)
 - `kernel` : covariance function, can be either a single kernel or a collection of kernels for multiclass and multi-outputs models
 - `likelihood` : likelihood of the model, currently implemented : Gaussian, Bernoulli (with logistic link), Multiclass (softmax or logistic-softmax) see [`Likelihood Types`](@ref likelihood_user)
 - `inference` : inference for the model, can be analytic, numerical or by sampling, check the model documentation to know what is available for your likelihood see the [`Compatibility Table`](@ref compat_table)

**Keyword arguments**

 - `verbose` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
- `optimizer` : Optimizer for kernel hyperparameters (to be selected from [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl))
- `atfrequency` : Choose how many variational parameters iterations are between hyperparameters optimization
- `mean` : PriorMean object, check the documentation on it [`MeanPrior`](@ref meanprior)
 - `IndependentPriors` : Flag for setting independent or shared parameters among latent GPs
 - `ArrayType` : Option for using different type of array for storage (allow for GPU usage)
"""
mutable struct MCGP{T<:Real,TLikelihood<:Likelihood{T},TInference<:Inference{T},TGP<:Abstract_GP{T},N} <: AbstractGP{T,TLikelihood,TInference,TGP,N}
    X::Matrix{T} #Feature vectors
    y::LatentArray #Output (-1,1 for classification, real for regression, matrix for multiclass)
    nSamples::Int64 # Number of data points
    nDim::Int64 # Number of covariates per data point
    nFeatures::Int64 # Number of features of the GP (equal to number of points)
    nLatent::Int64 # Number pf latent GPs
    f::NTuple{N,_MCGP} # Vector of latent GPs
    likelihood::TLikelihood
    inference::TInference
    verbose::Int64 #Level of printing information
    atfrequency::Int64
    Trained::Bool
end

function MCGP(X::AbstractArray{T},y::AbstractVector,kernel::Kernel,
            likelihood::Union{TLikelihood,Distribution},inference::TInference;
            verbose::Int=0,optimizer::Union{Bool,Optimizer,Nothing}=Adam(α=0.01),atfrequency::Integer=1,
            mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(), variance::Real=1.0,
            ArrayType::UnionAll=Vector) where {T<:Real,TLikelihood<:Likelihood,TInference<:SamplingInference}

            X,y,nLatent,likelihood = check_data!(X,y,likelihood)
            @assert check_implementation(:MCGP,likelihood,inference) "The $likelihood is not compatible or implemented with the $inference"

            nFeatures = nSamples = size(X,1); nDim = size(X,2);
            if isa(optimizer,Bool)
                optimizer = optimizer ? Adam(α=0.01) : nothing
            end

            if typeof(mean) <: Real
                mean = ConstantMean(mean)
            elseif typeof(mean) <: AbstractVector{<:Real}
                mean = EmpiricalMean(mean)
            end

            latentf = ntuple(_->_MCGP{T}(nFeatures,kernel,mean,variance),nLatent)

            likelihood = init_likelihood(likelihood,inference,nLatent,nSamples,nFeatures)
            inference = tuple_inference(inference,nLatent,nSamples,nSamples)
            inference.xview = view(X,:,:)
            inference.yview = view_y(likelihood,y,1:nSamples)
            MCGP{T,TLikelihood,typeof(inference),_MCGP{T},nLatent}(X,y,
                    nFeatures, nDim, nFeatures, nLatent,
                    latentf,likelihood,inference,
                    verbose,atfrequency,false)
end

function Base.show(io::IO,model::MCGP{T,<:Likelihood,<:Inference}) where {T}
    print(io,"Monte Carlo Gaussian Process with a $(model.likelihood) sampled via $(model.inference) ")
end

get_f(model::MCGP) = getproperty.(model.f,:f)
get_y(model::MCGP) = model.inference.yview
get_Z(model::MCGP) = [model.inference.xview]


text:
abstract type AbstractFit end

# # modified from `Figure` in PyPlot
"""
    Fit <: AbstractFit

struct for fit with individual parameters.
"""
mutable struct Fit <: AbstractFit
    o::PyObject
end

"""
    ArrayFit <: AbstractFit

struct for fit with parameters collected in an `Array`.
"""
mutable struct ArrayFit <: AbstractFit
    o::PyObject
end

# modified from `Figure` in PyPlot
PyObject(f::T) where {T<:AbstractFit} = getfield(f, :o)
convert(::Type{T}, o::PyObject) where {T<:AbstractFit} = T(o)
==(f::T, g::T) where {T<:AbstractFit} = PyObject(f) == PyObject(g)
==(f::AbstractFit, g::PyObject) = PyObject(f) == g
==(f::PyObject, g::AbstractFit) = f == PyObject(g)
hash(f::AbstractFit) = hash(PyObject(f))
pycall(f::AbstractFit, args...; kws...) = pycall(PyObject(f), args...; kws...)
(f::AbstractFit)(args...; kws...) = pycall(PyObject(f), PyAny, args...; kws...)
Base.Docs.doc(f::AbstractFit) = Base.Docs.doc(PyObject(f))
#
Base.setproperty!(f::AbstractFit, s::Symbol, x) = setproperty!(PyObject(f), s, x)
Base.setproperty!(f::AbstractFit, s::AbstractString, x) = setproperty!(PyObject(f), s, x)
hasproperty(f::AbstractFit, s::Symbol) = hasproperty(PyObject(f), s)
hasproperty(f::AbstractFit, s::AbstractString) = hasproperty(PyObject(f), s)
Base.propertynames(f::AbstractFit) = propertynames(PyObject(f))
haskey(f::AbstractFit, x) = haskey(PyObject(f), x)

# let the matrix property to be given as an array
function Base.getproperty(f::AbstractFit, s::Symbol)
    # use PyObject(f) instead of f to prevent StackOverflowError
    s === :matrix ? pycall(PyObject(f).matrix, PyAny) : getproperty(PyObject(f), s)
end
function Base.getproperty(f::AbstractFit, s::AbstractString)
    s === "matrix" ? pycall(PyObject(f).matrix, PyAny) : getproperty(PyObject(f), s)
end


text:
using MyGraph

G = DictGraph()

buildconnectedgraph!(collect(1:6), (a,b) -> first(abs.(rand(Int16, 1))), G)

D1 = dijkstra(G;source = 1)
D2 = dijkstra(G;source = 2)
D3 = dijkstra(G;source = 3)
D4 = dijkstra(G;source = 4)
D5 = dijkstra(G;source = 5)
D6 = dijkstra(G;source = 6)


text:
#module Scaling

export NiceScale, scaled

floorx(x) = x<0 ? ceil(x) : floor(x)
ceilx(x) = x<0 ? floor(x) : ceil(x)

mutable struct NiceScale
  minPoint::Float32
  maxPoint::Float32
  maxTicks::Float32
  tickSpacing::Float32
  range::Float32
  niceMin::Float32
  niceMax::Float32
  xfact::Float32
  xoffset::Float32

function NiceScale(min::Real,max::Real; padding=0.05,steps=10.)
     mn = min-padding*(max-min)
     mx = max+padding*(max-min)
     if mx==mn  mn=mn-1; mx=mx+1; end
     a=new(mn,mx,steps,0.,0.,0.,0.)
     a.range = niceNum(a.maxPoint-a.minPoint)
     a.tickSpacing = niceNum( a.range / (a.maxTicks-1),true)
     a.niceMin = ceil(a.minPoint / a.tickSpacing ) * a.tickSpacing
     a.niceMax = floor(a.maxPoint / a.tickSpacing ) * a.tickSpacing
     a.xoffset = (a.minPoint+a.maxPoint)/2
     a.xfact   = (a.maxPoint-a.minPoint)
     return a
end

end # NiceScale

Base.iterate(n::NiceScale,v=1) = v>n.niceMax ? nothing : (v,v+n.tickSpacing)
# Base.start(n::NiceScale) = n.niceMin
# Base.next(n::NiceScale,v) = (v,v+n.tickSpacing)
# Base.done(n::NiceScale,v) = v>n.niceMax
Base.length(n::NiceScale) = Int64((n.niceMax-n.niceMin)/n.tickSpacing)+1

function niceNum(range::Real,round::Bool=false)
   exponent = floor(log10(range))
   fraction = range / 10^exponent

   if round
     if fraction <= 1.5 niceFraction = 1
     elseif fraction <= 3 niceFraction = 2
     elseif fraction <= 7 niceFraction = 5
     else niceFraction = 5
     end
   else
     if fraction <= 1 niceFraction = 1
     elseif fraction <= 2 niceFraction = 2
     elseif fraction <= 5 niceFraction = 5
     else niceFraction = 10
     end
   end

   return niceFraction * 10.0^exponent

end

scaled(x::Real,n::NiceScale) = (x-n.xoffset)/n.xfact

#end # Scaling

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
"""
    abstract type AbstractCallback

Supertype of [`SafeCallback`](#)/`Callback`. When implementing callbacks,
you should subtype [`SafeCallback`](#) instead.
"""
abstract type AbstractCallback end

abstract type SafeCallback <: AbstractCallback end

"""
    abstract type Callback

Supertype of all callbacks. Callbacks add custom functionality to
the training loop by hooking into different [`Events.Event`](#)s

Any `Callback` can be used by passing it to [`Learner`](#). See
`subtypes(FluxTraining.Callback)` for implementations.

## Extending

See [Custom callbacks](/documents/docs/callbacks/custom.md) for a less
succinct tutorial format.

1. Create a `struct MyCallback` that subtypes `FluxTraining.Callback`.
2. Add event handlers by implementing methods for
   [`on`](#)`(event, phase, callback, learner)`. Methods should always
   dispatch on your callback, and may dispatch on specific [`Phases.Phase`](#)s
   and [`Events.Event`](#)s.

   For example, to implement an event handler that runs at the end of every
   step during training: `on(::StepEnd, ::AbstractTrainingPhase, ::MyCallback, learner)`.
3. Define what state the callback accesses and/or modifies by implementing
   [`stateaccess`](#)`(::MyCallback)`. While `learner` is always passed as an
   argument to `on` event handlers, by default a callback can not read or write
   to its fields. See [`stateaccess`](#) for more detail.

   If a callback needs to write some state that other callbacks should be able
   to access, it can store it in `learner.cbstate` if you add a permission in
   `stateaccess`.
4. If the callback needs some one-time initialization, you can implement [`init!`](#)
   which will be run at least once before any step is run.
"""
const Callback = SafeCallback

abstract type UnsafeCallback <: AbstractCallback end

# TODO: implement proper merging of permissions

"""
    stateaccess(callback)

Return a named tuple determining what learner state `callback`
can access. The default is `(;)`, the empty named tuple, meaning
no state can be accessed. Implementations of `stateaccess` should
always return the least permissions possible.

## Extending

For example, the [`ToGPU`](#) callback needs to write both the model and
the batch data, so its `stateaccess` implementation is:

```julia
stateaccess(::ToGPU) = (
    model = Write(),
    params = Write(),
    step = (xs = Write(), ys = Write()),
)
```

When defining `stateaccess`, be careful that you do return a `NamedTuple`.
`(x = Read(),)` is one but `(x = Read())` (without the comma) is parsed as
an assignment with value `Read()`.
"""
stateaccess(::Callback) = (;)
runafter(::AbstractCallback) = ()

"""
    abstract type ConflictResolution

A conflict resolution strategy for resolving write/write conflicts
of two callbacks.

See [`resolveconflict`](#).
"""
abstract type ConflictResolution end

"""
    abstract type NotDefined <: ConflictResolution

The default implementation of [`resolveconflict`](#). If a conflict
is detected, this ensures an error message is printed.
"""
struct NotDefined <: ConflictResolution end

"""
    abstract type Unresolvable <: ConflictResolution

Return from [`resolveconflict`](#) to indicate that two callbacks
are incompatible and cannot be used together.
"""
struct Unresolvable <: ConflictResolution end

"""
    abstract type RunFirst <: ConflictResolution

Return `RunFirst(cb1/cb2)` from [`resolveconflict`](#)`(cb1, cb2)` to indicate
that one of the callbacks should always run before the other.
"""
struct RunFirst <: ConflictResolution cb end

"""
    abstract type NoConflict <: ConflictResolution

Return from [`resolveconflict`](#) to indicate that, while the callbacks
modify the same state, they can be used together without any problems.
"""
struct NoConflict <: ConflictResolution end

function _resolveconflict(cb1, cb2)
    r = resolveconflict(cb1, cb2)
    if r === NotDefined()
        return resolveconflict(cb2, cb1)
    else
        return r
    end
end

"""
    resolveconflict(cb1, cb2)

Define a conflict resolution strategy for resolving a write/write conflict
between two callbacks.

The default is [`NotDefined()`], which will result in an error and a message
to implement this method.

To implement, dispatch on the callback types that you which to resolve (in any
order) and return one of the following:

- [`Unresolvable`](#)`()` if the callbacks must not be used together
- [`RunFirst`](#)`(cb)` if one of the callbacks needs to run first; or
- [`NoConflict`](#)`()` if the callbacks may run together in any order
"""
resolveconflict(::AbstractCallback, ::AbstractCallback) = NotDefined()

"""
    abstract type FitException

Abstract types for exceptions that can be thrown during fitting,
to change its control flow.

See [`CancelStepException`](#), [`CancelEpochException`](#), [`CancelFittingException`](#).
"""
abstract type FitException <: Exception end

"""
    CancelStepException(message)

Throw during fitting to cancel the currently running step.
This prematurely ends the current step without throwing an error.
Must be thrown inside the context of [`runstep`](#).

## Examples

```julia
runepoch(learner, phase) do _
    for (xs, ys) in batches
        runstep(learner, phase, (; xs, ys)) do _, state
            # training logic...
            if isnan(state.loss).
                throw(CancelStepException("Skipping NaN loss"))
            end
        end
    end
end
```
"""
struct CancelStepException <: FitException
    msg::String
end

"""
    CancelEpochException(message)

Throw during fitting to cancel the currently running epoch.
This prematurely ends the current epoch without throwing an error.
Must be thrown inside the context of [`runepoch`](#).

## Examples

```julia
runepoch(learner, phase) do _
    for batch in batches
        step!(learner, phase, batch)
        if learner.step.loss < 1.
            throw(CancelEpochException("Reached target loss"))
        end
    end
end
```
"""
struct CancelEpochException <: FitException
    msg::String
end

"""
    CancelFittingException(msg)

Throw during fitting to cancel it.
"""
struct CancelFittingException <: FitException
    msg::String
end

# Callback hook

"""
    on(event::Event, phase::Phase, callback::AbstractCallback, learner)

Handle `event` with [`Callback`](#) `callback`.
By default, this event handler does nothing for a callback.

To see events which an `AbstractCallback` handles, use

```julia
methods(Training.on, (Any, Any, MyCallbackType, Any)
```

## Extending

You can add event handlers to [`Callback`](#)s by implementing a method for `on`.
See also [`Callback`](#) and [custom callbacks](documents/docs/callbacks/custom.md).

A method of `on` should *always* dispatch on the callback type, i.e.
`on(event, phase, cb::MyCallback, learner)`. It may also dispatch on specific
[`Event`](#)s and [`Phase`](#). It should not dispatch on a specific type for
`learner`.
"""
on(::Event, phase, ::Callback, learner) = return

_on(e, p, cb, learner) = on(e, p, cb, learner)
function _on(e, p, cb::SafeCallback, learner)
    perms = Zygote.ignore() do
        stateaccess(cb)
    end
    on(e, p, cb, protect(learner, perms))
end

"""
    init!(callback, learner)

Initialize a callback. Default is to do nothing.

## Extending

To extend for a callback, implement `init!(cb::MyCallback, learner)`.
`init!` can set up internal state of a callback that depends on `learner`
and can also initialize shared callback state in `learner.cbstate`.
Just like `on` event handlers, the state access permissions must be correctly
defined using [`stateaccess`](#) to do so.

`init!` must also be idempotent, i.e. running it twice on the same `Learner`
should have the same effect as runnning it once.
"""
init!(::Callback, learner) = return


text:
# LANGUAGE: Julia
# AUTHOR: Mahendra Yadav
println("Hello, World!")

#  Execution
# julia hello.jl


text:
module LBPartialBase

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export cons_lb_partial_base

function cons_lb_partial_base(m)

    @expression(m, pb_mass_on_top_of_i[i=1:length(boxes), j=1:length(boxes); i!=j],
        mass_density[i] * m[:pbs_box_ij][i, j] / area_of_box[i]
    )

    @expression(m, pb_sum_of_mass_on_top_of_i[j=1:length(boxes)],
        sum([m[:pb_mass_on_top_of_i][i, j] for i in 1:length(boxes) if i!=j])
    )

    @constraint(m, c_max_allowed_weight[i=1:length(boxes)-1],
        m[:pb_sum_of_mass_on_top_of_i][i] <= max_allowed_weight[i]
    )

    return m
end

end


text:
# Sparisty-promoting LS-RTM of the 2D Marmousi model with on-the-fly Fourier transforms
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using Statistics, Random, LinearAlgebra
using JUDI, SegyIO, HDF5, JOLI, PyPlot

# Load migration velocity model
if ~isfile("marmousi_model.h5")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_model.h5`)
end
n, d, o, m0 = read(h5open("marmousi_model.h5", "r"), "n", "d", "o", "m0")

# Set up model structure
model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)

# Load data
if ~isfile("marmousi_2D.segy")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_2D.segy`)
end
block = segy_read("marmousi_2D.segy")
d_lin = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceDepth")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.03)    # 30 Hz wavelet
q = judiVector(src_geometry, wavelet)

# Set up info structure
ntComp = get_computational_nt(q.geometry, d_lin.geometry, model0)  # no. of computational time steps
info = Info(prod(model0.n), d_lin.nsrc, ntComp)

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

# Setup operators
opt = Options(optimal_checkpointing=false)
M = judiModeling(info, model0, q.geometry, d_lin.geometry; options=opt)
J = judiJacobian(M, q)

# Right-hand preconditioners (model topmute)
Mr = judiTopmute(model0.n, 52, 10)

# Generate distribution of frequecies from source wavelet
q_dist = generate_distribution(q)

# Linearized Bregman parameters
x = zeros(Float32, info.n)
z = zeros(Float32, info.n)
batchsize = 10
niter = 20
fval = zeros(Float32, niter)
J.options.frequencies = Array{Any}(undef, d_lin.nsrc)

# Soft thresholding functions and Curvelet transform
soft_thresholding(x::Array{Float64}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float64, lambda), 0.0)
soft_thresholding(x::Array{Float32}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float32, lambda), 0f0)
C = joCurvelet2D(model0.n[1], model0.n[2]; zero_finest = true, DDT = Float32, RDT = Float64)
lambda = []
t = 1f-4

# Main loop
for j = 1: niter
    println("Iteration: ", j)

    # Select batch and set up left-hand preconditioner
    i = randperm(d_lin.nsrc)[1:batchsize]
    Ml = judiMarineTopmute2D(30, d_lin[i].geometry)    # data topmute

    # Select randomly selected batch of frequencies for each shot
    for k = 1: d_lin.nsrc
        J.options.frequencies[k] = select_frequencies(q_dist; fmin = 0.004f0, fmax = 0.05f0, nf=10)
    end

    # Compute residual and gradient
    r = Ml*J[i]*Mr*x - Ml*d_lin[i]
    g = adjoint(Mr)*adjoint(J[i])*adjoint(Ml)*r

    # Step size and update variable
    fval[j] = .5f0*norm(r)^2
    j == 1 && (global lambda = 0.5*norm(C*t*g, Inf))   # estimate thresholding parameter in 1st iteration

    # Update variables and save snapshot
    global z -= t*g
    global x = adjoint(C)*soft_thresholding(C*z, lambda)
end

# Save final velocity model, function value and history
h5open("lsrtm_marmousi_frequency_result.h5", "w") do file
    write(file, "x", reshape(x, model0.n), "fval", fval)
end

# Plot final image
figure(); imshow(copy(adjoint(reshape(x, model0.n))), extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
title("SPLS-RTM with on on-the-fly DFTs"); xlabel("Lateral position [km]"); ylabel("Depth [km]")

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ManagementPolicySnapShot <: SwaggerModel
    delete::Any # spec type: Union{ Nothing, DateAfterCreation } # spec name: delete

    function ManagementPolicySnapShot(;delete=nothing)
        o = new()
        validate_property(ManagementPolicySnapShot, Symbol("delete"), delete)
        setfield!(o, Symbol("delete"), delete)
        o
    end
end # type ManagementPolicySnapShot

const _property_map_ManagementPolicySnapShot = Dict{Symbol,Symbol}(Symbol("delete")=>Symbol("delete"))
const _property_types_ManagementPolicySnapShot = Dict{Symbol,String}(Symbol("delete")=>"DateAfterCreation")
Base.propertynames(::Type{ ManagementPolicySnapShot }) = collect(keys(_property_map_ManagementPolicySnapShot))
Swagger.property_type(::Type{ ManagementPolicySnapShot }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagementPolicySnapShot[name]))}
Swagger.field_name(::Type{ ManagementPolicySnapShot }, property_name::Symbol) =  _property_map_ManagementPolicySnapShot[property_name]

function check_required(o::ManagementPolicySnapShot)
    true
end

function validate_property(::Type{ ManagementPolicySnapShot }, name::Symbol, val)
end


text:
# TODO: Document this.
macro group_by(src, exprs...)
    Expr(
        :call,
        :GroupBy,
        esc(src),
        Expr(
            :tuple,
            expression_expr.(exprs)...,
        ),
    )
end


text:
using TrueRandom
using Distributions
using Test

@testset "ANU QRNG" begin
    r = ANUQRNG(1024, 768; prefill=false)
    TrueRandom.prefillcache(r)

    @test length(r.cache) == r.batchsize
    @test length(r.cacheprefill) == 0
    @test r.pointer == 1

    @test rand(r, UInt8) isa UInt8
    @test rand(r, UInt8, 2) isa Vector{UInt8}
    @test rand(r, UInt16) isa UInt16
    @test rand(r, UInt16, 2) isa Vector{UInt16}
    @test rand(r, UInt32) isa UInt32
    @test rand(r, UInt32, 2) isa Vector{UInt32}
    @test rand(r, UInt64) isa UInt64
    @test rand(r, UInt64, 2) isa Vector{UInt64}
    @test rand(r, UInt128) isa UInt128
    @test rand(r, UInt128, 2) isa Vector{UInt128}

    @test r.pointer == 48

    @test rand(r, Int8) isa Int8
    @test rand(r, Int8, 2) isa Vector{Int8}
    @test rand(r, Int16) isa Int16
    @test rand(r, Int16, 2) isa Vector{Int16}
    @test rand(r, Int32) isa Int32
    @test rand(r, Int32, 2) isa Vector{Int32}
    @test rand(r, Int64) isa Int64
    @test rand(r, Int64, 2) isa Vector{Int64}
    @test rand(r, Int128) isa Int128
    @test rand(r, Int128, 2) isa Vector{Int128}

    @test r.pointer == 95

    @test rand(r, Float32) isa Float32
    @test rand(r, Float32, 2) isa Vector{Float32}
    @test rand(r, Float64) isa Float64
    @test rand(r, Float64, 2) isa Vector{Float64}

    @test r.pointer == 113
    @test length(r.cacheprefill) == 0

    m = rand(r, Float64, 5, 8)
    @test m isa Matrix{Float64}
    @test 0 <= minimum(m) <= 1
    @test 0 <= maximum(m) <= 1

    p1 = r.pointer
    @test rand(r, UInt16, 1024) isa Any
    @test r.pointer == p1

    @test rand(r, UInt16, 2048) isa Any
    @test r.pointer == p1

    @test randn(r) isa Any
    @test rand(r, Normal()) isa Any
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using Printf
using JSON
using DataStructures
import Base: getindex, time
import GZip

abstract type UCElement end

abstract type Time <: UCElement end

mutable struct Bus <: UCElement
    name::String
    offset::Int
    load::Array{Float64}
    units::Array
    price_sensitive_loads::Array
end

mutable struct CostSegment <: UCElement
    mw::Array{Float64}
    cost::Array{Float64}
end

mutable struct StartupCategory <: UCElement
    delay::Int
    cost::Float64
end

mutable struct Unit <: UCElement
    name::String
    bus::Bus
    max_power::Array{Float64}
    min_power::Array{Float64}
    must_run::Array{Bool}
    min_power_cost::Array{Float64}
    cost_segments::Array{CostSegment}
    min_uptime::Int
    min_downtime::Int
    ramp_up_limit::Float64
    ramp_down_limit::Float64
    startup_limit::Float64
    shutdown_limit::Float64
    initial_status::Union{Int,Nothing}
    initial_power::Union{Float64,Nothing}
    provides_spinning_reserves::Array{Bool}
    startup_categories::Array{StartupCategory}
end # Unit

mutable struct TransmissionLine <: UCElement
    name::String
    offset::Int
    source::Bus
    target::Bus
    reactance::Float64
    susceptance::Float64
    normal_flow_limit::Array{Float64}
    emergency_flow_limit::Array{Float64}
    flow_limit_penalty::Array{Float64}
end

mutable struct Reserves <: UCElement
    spinning::Array{Float64}
end

mutable struct Contingency <: UCElement
    name::String
    lines::Array{TransmissionLine}
    units::Array{Unit}
end

mutable struct PriceSensitiveLoad <: UCElement
    name::String
    bus::Bus
    demand::Array{Float64}
    revenue::Array{Float64}
end

mutable struct UnitCommitmentInstance
    time::Int
    power_balance_penalty::Array{Float64}
    "Penalty for failing to meet reserve requirement."
    shortfall_penalty::Array{Float64}
    units::Array{Unit}
    buses::Array{Bus}
    lines::Array{TransmissionLine}
    reserves::Reserves
    contingencies::Array{Contingency}
    price_sensitive_loads::Array{PriceSensitiveLoad}
end

function Base.show(io::IO, instance::UnitCommitmentInstance)
    print(io, "UnitCommitmentInstance with ")
    print(io, "$(length(instance.units)) units, ")
    print(io, "$(length(instance.buses)) buses, ")
    print(io, "$(length(instance.lines)) lines, ")
    print(io, "$(length(instance.contingencies)) contingencies, ")
    print(io, "$(length(instance.price_sensitive_loads)) price sensitive loads")
end

function read_benchmark(name::AbstractString) :: UnitCommitmentInstance
    basedir = dirname(@__FILE__)
    return UnitCommitment.read("$basedir/../instances/$name.json.gz")
end

function read(path::AbstractString)::UnitCommitmentInstance
    if endswith(path, ".gz")
        return read(GZip.gzopen(path))
    else
        return read(open(path))
    end
end

function read(file::IO)::UnitCommitmentInstance
    return from_json(JSON.parse(file, dicttype=()->DefaultOrderedDict(nothing)))
end

function from_json(json; fix=true)
    units = Unit[]
    buses = Bus[]
    contingencies = Contingency[]
    lines = TransmissionLine[]
    loads = PriceSensitiveLoad[]
    T = json["Parameters"]["Time (h)"]

    name_to_bus = Dict{String, Bus}()
    name_to_line = Dict{String, TransmissionLine}()
    name_to_unit = Dict{String, Unit}()

    function timeseries(x; default=nothing)
        x !== nothing || return default
        x isa Array || return [x for t in 1:T]
        return x
    end

    function scalar(x; default=nothing)
        x !== nothing || return default
        x
    end

    # Read parameters
    power_balance_penalty = timeseries(json["Parameters"]["Power balance penalty (\$/MW)"],
                                       default=[1000.0 for t in 1:T])
    shortfall_penalty = timeseries(json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
                                   default=[0. for t in 1:T])

    # Read buses
    for (bus_name, dict) in json["Buses"]
        bus = Bus(bus_name,
                  length(buses),
                  timeseries(dict["Load (MW)"]),
                  Unit[],
                  PriceSensitiveLoad[])
        name_to_bus[bus_name] = bus
        push!(buses, bus)
    end

    # Read units
    for (unit_name, dict) in json["Generators"]
        bus = name_to_bus[dict["Bus"]]

        # Read production cost curve
        K = length(dict["Production cost curve (MW)"])
        curve_mw   = hcat([timeseries(dict["Production cost curve (MW)"][k]) for k in 1:K]...)
        curve_cost = hcat([timeseries(dict["Production cost curve (\$)"][k]) for k in 1:K]...)
        min_power = curve_mw[:, 1]
        max_power = curve_mw[:, K]
        min_power_cost = curve_cost[:, 1]
        segments = CostSegment[]
        for k in 2:K
            amount = curve_mw[:, k] - curve_mw[:, k-1]
            cost = (curve_cost[:, k] - curve_cost[:, k-1]) ./ amount
            replace!(cost, NaN=>0.0)
            push!(segments, CostSegment(amount, cost))
        end

        # Read startup costs
        startup_delays = scalar(dict["Startup delays (h)"], default=[1])
        startup_costs  = scalar(dict["Startup costs (\$)"], default=[0.])
        startup_categories = StartupCategory[]
        for k in 1:length(startup_delays)
            push!(startup_categories, StartupCategory(startup_delays[k],
                                                      startup_costs[k]))
        end

        # Read and validate initial conditions
        initial_power = scalar(dict["Initial power (MW)"], default=nothing)
        initial_status = scalar(dict["Initial status (h)"], default=nothing)
        if initial_power === nothing
            initial_status === nothing || error("unit $unit_name has initial status but no initial power")
        else
            initial_status !== nothing || error("unit $unit_name has initial power but no initial status")
            initial_status != 0 || error("unit $unit_name has invalid initial status")
            if initial_status < 0 && initial_power > 1e-3
                error("unit $unit_name has invalid initial power")
            end
        end

        unit = Unit(unit_name,
                    bus,
                    max_power,
                    min_power,
                    timeseries(dict["Must run?"], default=[false for t in 1:T]),
                    min_power_cost,
                    segments,
                    scalar(dict["Minimum uptime (h)"], default=1),
                    scalar(dict["Minimum downtime (h)"], default=1),
                    scalar(dict["Ramp up limit (MW)"], default=1e6),
                    scalar(dict["Ramp down limit (MW)"], default=1e6),
                    scalar(dict["Startup limit (MW)"], default=1e6),
                    scalar(dict["Shutdown limit (MW)"], default=1e6),
                    initial_status,
                    initial_power,
                    timeseries(dict["Provides spinning reserves?"],
                               default=[true for t in 1:T]),
                    startup_categories)
        push!(bus.units, unit)
        name_to_unit[unit_name] = unit
        push!(units, unit)
    end

    # Read reserves
    reserves = Reserves(zeros(T))
    if "Reserves" in keys(json)
        reserves.spinning = timeseries(json["Reserves"]["Spinning (MW)"],
                                       default=zeros(T))
    end

    # Read transmission lines
    if "Transmission lines" in keys(json)
        for (line_name, dict) in json["Transmission lines"]
            line = TransmissionLine(line_name,
                                    length(lines) + 1,
                                    name_to_bus[dict["Source bus"]],
                                    name_to_bus[dict["Target bus"]],
                                    scalar(dict["Reactance (ohms)"]),
                                    scalar(dict["Susceptance (S)"]),
                                    timeseries(dict["Normal flow limit (MW)"],
                                               default=[1e8 for t in 1:T]),
                                    timeseries(dict["Emergency flow limit (MW)"],
                                               default=[1e8 for t in 1:T]),
                                    timeseries(dict["Flow limit penalty (\$/MW)"],
                                               default=[5000.0 for t in 1:T]))
            name_to_line[line_name] = line
            push!(lines, line)
        end
    end

    # Read contingencies
    if "Contingencies" in keys(json)
        for (cont_name, dict) in json["Contingencies"]
            affected_units = Unit[]
            affected_lines = TransmissionLine[]
            if "Affected lines" in keys(dict)
                affected_lines = [name_to_line[l] for l in dict["Affected lines"]]
            end
            if "Affected units" in keys(dict)
                affected_units = [name_to_unit[u] for u in dict["Affected units"]]
            end
            cont = Contingency(cont_name, affected_lines, affected_units)
            push!(contingencies, cont)
        end
    end

    # Read price-sensitive loads
    if "Price-sensitive loads" in keys(json)
        for (load_name, dict) in json["Price-sensitive loads"]
            bus = name_to_bus[dict["Bus"]]
            load = PriceSensitiveLoad(load_name,
                                      bus,
                                      timeseries(dict["Demand (MW)"]),
                                      timeseries(dict["Revenue (\$/MW)"]),
                                     )
            push!(bus.price_sensitive_loads, load)
            push!(loads, load)
        end
    end

    instance = UnitCommitmentInstance(T,
                                      power_balance_penalty,
                                      shortfall_penalty,
                                      units,
                                      buses,
                                      lines,
                                      reserves,
                                      contingencies,
                                      loads)
    if fix
        UnitCommitment.fix!(instance)
    end
    return instance
end

"""
    slice(instance, range)

Creates a new instance, with only a subset of the time periods.
This function does not modify the provided instance. The initial
conditions are also not modified.

Example
-------

    # Build a 2-hour UC instance
    instance = UnitCommitment.read_benchmark("test/case14")
    modified = UnitCommitment.slice(instance, 1:2)

"""
function slice(instance::UnitCommitmentInstance, range::UnitRange{Int})::UnitCommitmentInstance
    modified = deepcopy(instance)
    modified.time = length(range)
    modified.power_balance_penalty = modified.power_balance_penalty[range]
    modified.reserves.spinning = modified.reserves.spinning[range]
    for u in modified.units
        u.max_power = u.max_power[range]
        u.min_power = u.min_power[range]
        u.must_run = u.must_run[range]
        u.min_power_cost = u.min_power_cost[range]
        u.provides_spinning_reserves = u.provides_spinning_reserves[range]
        for s in u.cost_segments
            s.mw = s.mw[range]
            s.cost = s.cost[range]
        end
    end
    for b in modified.buses
        b.load = b.load[range]
    end
    for l in modified.lines
        l.normal_flow_limit = l.normal_flow_limit[range]
        l.emergency_flow_limit = l.emergency_flow_limit[range]
        l.flow_limit_penalty = l.flow_limit_penalty[range]
    end
    for ps in modified.price_sensitive_loads
        ps.demand = ps.demand[range]
        ps.revenue = ps.revenue[range]
    end
    return modified
end

export UnitCommitmentInstance


text:
function scn(a::Number; digits::Int=1, exp_digits::Int=1)
	typeof(a) <: Float64 ? nothing : return nothing
end

function scn(a::Float64; digits::Int=1, exp_digits::Int=1)
	@assert digits >= 0
    # a = m x 10^e
    if a == 0
        e = 0
        m = 0.0
    else
        e = Int(floor(log(abs(a))/log(10)))
        m = a*exp(-e*log(10))
    end

    m = round(m, digits=digits)
	if m == 10.0
		m = 1.0
		e += 1
	end
    if digits == 0
        m = Int(floor(m))
		strm = string(m)
	else
		strm = string(m)
		is_neg = m < 0.
		strm = strm*"0"^max(0, 2+digits+is_neg-length(strm))
    end
    sgn = a >= 0 ? " " : ""
    sgne = e >= 0 ? "+" : "-"

	stre = string(abs(e))
	stre = "0"^max(0, exp_digits - length(stre)) * stre
    return "$sgn$(strm)e$sgne$(stre)"
end

function delta!(Δx::SizedArray{Tuple{nx},T,1,1}, x::SizedArray{Tuple{nx},T,1,1},
    x_ref::SizedArray{Tuple{nx},T,1,1}) where {nx,T}
    Δx .= x
    Δx .-= x_ref
    return nothing
end

function delta!(Δx::SizedArray{Tuple{nx},T,1,1}, x::AbstractArray{T},
    x_ref::AbstractArray{T}) where {nx,T}
    Δx .= x
    Δx .-= x_ref
    return nothing
end

function set!(s::SubArray, x::SizedArray{Tuple{nx},T,1,1}) where {nx,T}
    s .= x
    return nothing
end

function setminus!(s::SubArray, x::SizedArray{Tuple{nx},T,1,1}) where {nx,T}
    s .= -1.0 .* x
    return nothing
end

function convert_video_to_gif(video_file_path::String, output_path::String="output.gif";
    framerate::Int=30, start_time=0., duration=1e3, overwrite=false, width::Int=1080, height::Int=-2, hq_colors::Bool=false)
    output_path = abspath(output_path)

    if !isfile(video_file_path)
        error("Could not find the input file $video_file_path")
    end
    if isfile(output_path) && !overwrite
        error("The output path $output_path already exists. To overwrite that file, you can pass `overwrite=true` to this function")
    end

    mktempdir() do tmpdir
        # run(MeshCat.unpack_cmd(video_file_path, tmpdir, ".mp4", nothing)) # unpack the .tar file
        # cmd = ["-r", string(framerate), "-i", "%07d.png", "-vcodec", "libx264", "-preset", "slow", "-crf", "18"]
        color_map = hq_colors ?
            "[0:v] fps=$framerate, scale=$width:$height,split [a][b];[a] palettegen=stats_mode=single [p];[b][p] paletteuse=new=1" :
            "[0:v] fps=$framerate, scale=$width:$height,split [a][b];[a] palettegen [p];[b][p] paletteuse"
        cmd = ["-ss", string(start_time), "-t", string(duration), "-i", video_file_path, "-filter_complex", color_map]
        if overwrite
            push!(cmd, "-y")
        end
        push!(cmd, output_path)

        cd(tmpdir) do
            FFMPEG.exe(cmd...)
        end
    end
    @info("Saved output as $output_path")
    return output_path
end

function convert_frames_to_video_and_gif(filename, overwrite::Bool=true)
    MeshCat.convert_frames_to_video(
        homedir() * "/Downloads/$filename.tar",
        homedir() * "/Documents/video/$filename.mp4", overwrite=overwrite)

    convert_video_to_gif(
        homedir() * "/Documents/video/$filename.mp4",
        homedir() * "/Documents/video/$filename.gif", overwrite=overwrite)
    return nothing
end

function save_markdown(path::String, content::String; overwrite::Bool=true)
	mode = overwrite ? "w" : "a"
	io = open(path, mode)
	write(io, content)
	close(io)
	return nothing
end

function horizontal_line(n::Int)
	@assert n >= 1
	out = "|" * " --- |"^n
	return out
end

function content_line(c::AbstractVector)
	n = length(c)
	@assert n >= 1
	out = "|"
	for i = 1:n
		out *= string(c[i]) * "|"
	end
	return out
end

function save_expressions(expr::Dict{Symbol,Expr},
	path::AbstractString="expr.jld2"; overwrite::Bool=false)
	path = abspath(path)
    if isfile(path)
        if overwrite
            rm(path)
        end
        if !overwrite
            @warn "file exists -- not overwriting"
        end
    end
	@save path expr
	@info("Saved output as $path")
	return nothing
end

function load_expressions(path::AbstractString="expr.jld2")
	path = abspath(path)
	if !isfile(path)
		error("Could not find the input file $path")
	end
	@load path expr
	return expr
end

function module_dir()
	return joinpath(@__DIR__, "..")
end

function axes_pair_to_quaternion(n1, n2)
	if norm(n1 + n2, Inf) < 1e-5
		n2 = n2 + 1e-5ones(3)
	end

	reg(x) = 1e-20 * (x == 0) + x
	# provides the quaternion that rotates n1 into n2, assuming n1 and n2 are normalized
	n1 ./= reg(norm(n1))
	n2 ./= reg(norm(n2))
	n3 = skew(n1)*n2
	cθ = n1' * n2 # cosine
	sθ = norm(n3) # sine
	axis = n3 ./ reg(sθ)
	tanθhalf = sθ / reg(1 + cθ)
	q = [1; tanθhalf * axis]
	q /= norm(q)
	return Quaternion(q...)
end


text:
require("ArgParse")
require("src/Spark")
using ArgParse
using Spark

function parse_cli()
  s = ArgParseSettings("start_worker.jl")

  @add_arg_table s begin
    "port"
      help = "Listening port on worker for RPC calls"
      arg_type = Int64
      default = 6666
  end

  return parse_args(s)
end

function main()
  args = parse_cli()
  port = args["port"]
  worker = Worker("127.0.0.1", port)
  start(worker)
  println("done")
end

main()