text:
"""
Unpack a rotation-parameterized BVN covariance matrix.

Args:
 - ab: The ratio of the minor to the major axis.
 - angle: Rotation angle (in radians)
 - scale: The major axis.

Returns:
  The 2x2 covariance matrix parameterized by the inputs.
"""
function get_bvn_cov{NumType <: Number}(
    ab::NumType, angle::NumType, scale::NumType)

   if NumType <: AbstractFloat
       @assert 0 < scale
       @assert 0 < ab <= 1.
   end

   cp = cos(angle)
   sp = sin(angle)
   ab_term = (ab ^ 2 - 1)
   scale_squared = scale ^ 2
   off_diag_term = -scale_squared * cp * sp * ab_term
   NumType[ scale_squared * (1 + ab_term * (sp ^ 2))    off_diag_term;
            off_diag_term    scale_squared * (1 + ab_term * (cp ^ 2))]
end

"""
Pre-allocated memory for quantities related to derivatives of bivariate
normals.
"""
type BivariateNormalDerivatives{NumType <: Number}

  # Pre-allocated memory for py1, py2, and f when evaluating BVNs
  py1::Array{NumType, 1}
  py2::Array{NumType, 1}
  f_pre::Array{NumType, 1}

  # Derivatives of a bvn with respect to (x, sig).
  bvn_x_d::Array{NumType, 1}
  bvn_sig_d::Array{NumType, 1}
  bvn_xx_h::Array{NumType, 2}
  bvn_xsig_h::Array{NumType, 2}
  bvn_sigsig_h::Array{NumType, 2}

  # intermediate values used in d bvn / d(x, sig)
  dpy1_dsig::Array{NumType, 1}
  dpy2_dsig::Array{NumType, 1}

  # TODO: delete this, it is now in BvnComponent
  dsiginv_dsig::Array{NumType, 2}

  # Derivatives of a bvn with respect to (u, shape)
  bvn_u_d::Array{NumType, 1}
  bvn_uu_h::Array{NumType, 2}
  bvn_s_d::Array{NumType, 1}
  bvn_ss_h::Array{NumType, 2}
  bvn_us_h::Array{NumType, 2}

  function BivariateNormalDerivatives(ThisNumType::DataType)
    py1 = zeros(ThisNumType, 1)
    py2 = zeros(ThisNumType, 1)
    f_pre = zeros(ThisNumType, 1)

    bvn_x_d = zeros(ThisNumType, 2)
    bvn_sig_d = zeros(ThisNumType, 3)
    bvn_xx_h = zeros(ThisNumType, 2, 2)
    bvn_xsig_h = zeros(ThisNumType, 2, 3)
    bvn_sigsig_h = zeros(ThisNumType, 3, 3)

    dpy1_dsig = zeros(ThisNumType, 3)
    dpy2_dsig = zeros(ThisNumType, 3)
    dsiginv_dsig = zeros(ThisNumType, 3, 3)

    # Derivatives wrt u.
    bvn_u_d = zeros(ThisNumType, 2)
    bvn_uu_h = zeros(ThisNumType, 2, 2)

    # Shape deriviatives.  Here, s stands for "shape".
    bvn_s_d = zeros(ThisNumType, length(gal_shape_ids))

    # The hessians.
    bvn_ss_h = zeros(ThisNumType, length(gal_shape_ids), length(gal_shape_ids))
    bvn_us_h = zeros(ThisNumType, 2, length(gal_shape_ids))

    new(py1, py2, f_pre,
        bvn_x_d, bvn_sig_d, bvn_xx_h, bvn_xsig_h, bvn_sigsig_h,
        dpy1_dsig, dpy2_dsig,
        dsiginv_dsig,
        bvn_u_d, bvn_uu_h, bvn_s_d, bvn_ss_h, bvn_us_h)
  end
end

"""
Relevant parameters of a bivariate normal distribution.

Args:
  the_mean: The mean as a 2x1 column vector
  the_cov: The covaraiance as a 2x2 matrix
  weight: A scalar weight

Attributes:
   the_mean: The mean argument
   precision: The inverse of the_cov
   z: The weight times the normalizing constant.
"""
immutable BvnComponent{NumType <: Number}
    the_mean::Vector{NumType}
    precision::Matrix{NumType}
    z::NumType
    dsiginv_dsig::Matrix{NumType}

    function BvnComponent{T1 <: Number, T2 <: Number, T3 <: Number}(
        the_mean::Vector{T1}, the_cov::Matrix{T2}, weight::T3;
        calculate_siginv_deriv::Bool=true)

      ThisNumType = promote_type(T1, T2, T3)
      the_det = the_cov[1,1] * the_cov[2,2] - the_cov[1,2] * the_cov[2,1]
      c = 1 ./ (the_det^.5 * 2pi)

      if calculate_siginv_deriv
        # Derivatives of Sigma^{-1} with respect to sigma.  These are the second
        # derivatives of log|Sigma| with respect to sigma.
        # dsiginv_dsig[a, b] is the derivative of sig^{-1}[a] / d sig[b]
        dsiginv_dsig = Array(ThisNumType, 3, 3)

        precision = the_cov^-1

        dsiginv_dsig[1, 1] = -precision[1, 1] ^ 2
        dsiginv_dsig[1, 2] = -2 * precision[1, 1] * precision[1, 2]
        dsiginv_dsig[1, 3] = -precision[1, 2] ^ 2

        dsiginv_dsig[2, 1] = -precision[1, 1] * precision[2, 1]
        dsiginv_dsig[2, 2] =
          -(precision[1, 1] * precision[2, 2] + precision[1, 2] ^ 2)
        dsiginv_dsig[2, 3] = -precision[2, 2] * precision[1, 2]

        dsiginv_dsig[3, 1] = -precision[1, 2] ^ 2
        dsiginv_dsig[3, 2] = - 2 * precision[2, 2] * precision[2, 1]
        dsiginv_dsig[3, 3] = -precision[2, 2] ^ 2
        new{ThisNumType}(the_mean, precision, c * weight, dsiginv_dsig)
      else
        new{ThisNumType}(the_mean, the_cov^-1, c * weight, zeros(ThisNumType, 0, 0))
      end
    end
end

"""
Return quantities related to the pdf of an offset bivariate normal.

Args:
  - bvn_derivs: BivariateNormalDerivatives to be updated in place
  - bmc: A bivariate normal component
  - x: A 2x1 vector containing a mean offset to be applied to bmc

Returns:
  In bvn_derivs, sets these values in place:
  - py1: The first row of the precision times (x - the_mean)
  - py2: The second row of the precision times (x - the_mean)
  - The density of the bivariate normal times the weight.
"""
function eval_bvn_pdf!{NumType <: Number}(
    bvn_derivs::BivariateNormalDerivatives{NumType},
    bmc::BvnComponent{NumType}, x::Vector{Float64})

  bvn_derivs.py1[1] =
    bmc.precision[1,1] * (x[1] - bmc.the_mean[1]) +
    bmc.precision[1,2] * (x[2] - bmc.the_mean[2])
  bvn_derivs.py2[1] =
    bmc.precision[2,1] * (x[1] - bmc.the_mean[1]) +
    bmc.precision[2,2] * (x[2] - bmc.the_mean[2])
  bvn_derivs.f_pre[1] =
    bmc.z * exp(-0.5 * ((x[1] - bmc.the_mean[1]) * bvn_derivs.py1[1] +
                        (x[2] - bmc.the_mean[2]) * bvn_derivs.py2[1]))
end

##################
# Derivatives

"""
Calculate the value, gradient, and hessian of
  -0.5 * x' sigma^-1 x - 0.5 * log|sigma|
with respect to x and sigma.  This assumes that
  bvn_derivs.py1, .py2, and .f_pre have already been populated
  witha call to eval_bvn_pdf!.

Args:
  - elbo_vars: A data structure with pre-allocated intermediate variables.
  - bvn: A bivariate normal component to get derivatives for.
  - calculate_x_hess: Whether to calcualte x Hessian terms.
  - calculate_sigma_hess: Whether to calcualte sigma Hessian terms.
"""
function get_bvn_derivs!{NumType <: Number}(
    bvn_derivs::BivariateNormalDerivatives{NumType},
    bvn::BvnComponent{NumType}, calculate_x_hess::Bool,
    calculate_sigma_hessian::Bool)

  # Gradient with respect to x.
  bvn_x_d = bvn_derivs.bvn_x_d
  bvn_x_d[1] = -bvn_derivs.py1[1]
  bvn_x_d[2] = -bvn_derivs.py2[1]

  if calculate_x_hess
    bvn_xx_h = bvn_derivs.bvn_xx_h

    # Hessian terms involving only x
    bvn_xx_h[1, 1] = -bvn.precision[1, 1]
    bvn_xx_h[2, 2] = -bvn.precision[2, 2]
    bvn_xx_h[1, 2] = bvn_xx_h[2, 1] = -bvn.precision[1 ,2]
  end

  # The first term is the derivative of -0.5 * x' Sigma^{-1} x
  # The second term is the derivative of -0.5 * log|Sigma|
  bvn_sig_d = bvn_derivs.bvn_sig_d
  bvn_sig_d[1] =
    0.5 * bvn_derivs.py1[1] * bvn_derivs.py1[1] - 0.5 * bvn.precision[1, 1]
  bvn_sig_d[2] =
    bvn_derivs.py1[1] * bvn_derivs.py2[1]             - bvn.precision[1, 2]
  bvn_sig_d[3] =
    0.5 * bvn_derivs.py2[1] * bvn_derivs.py2[1] - 0.5 * bvn.precision[2, 2]

  if calculate_sigma_hessian

    # Hessian calculation for terms containing sigma.

    # Derivatives of py1 and py2 with respect to s11, s12, s22 in that order.
    # These are used for the hessian calculations.
    dpy1_dsig = bvn_derivs.dpy1_dsig
    dpy1_dsig[1] = -bvn_derivs.py1[1] * bvn.precision[1,1]
    dpy1_dsig[2] = -bvn_derivs.py2[1] * bvn.precision[1,1] -
                    bvn_derivs.py1[1] * bvn.precision[1,2]
    dpy1_dsig[3] = -bvn_derivs.py2[1] * bvn.precision[1,2]

    dpy2_dsig = bvn_derivs.dpy2_dsig
    dpy2_dsig[1] = -bvn_derivs.py1[1] * bvn.precision[1,2]
    dpy2_dsig[2] = -bvn_derivs.py1[1] * bvn.precision[2,2] -
                    bvn_derivs.py2[1] * bvn.precision[1,2]
    dpy2_dsig[3] = -bvn_derivs.py2[1] * bvn.precision[2,2]

    # Hessian terms involving only sigma
    bvn_sigsig_h = bvn_derivs.bvn_sigsig_h
    for s_ind=1:3
      # Differentiate with respect to s_ind second.
      bvn_sigsig_h[1, s_ind] = #bvn_sigsig_h[s_ind, 1] =
        bvn_derivs.py1[1] * dpy1_dsig[s_ind] - 0.5 * bvn.dsiginv_dsig[1, s_ind]

      # d log|sigma| / dsigma12 is twice lambda12.
      bvn_sigsig_h[2, s_ind] =
        bvn_derivs.py1[1] * dpy2_dsig[s_ind] + bvn_derivs.py2[1] * dpy1_dsig[s_ind] -
        bvn.dsiginv_dsig[2, s_ind]

      bvn_sigsig_h[3, s_ind] = #bvn_sigsig_h[s_ind, 3] =
        bvn_derivs.py2[1] * dpy2_dsig[s_ind] - 0.5 * bvn.dsiginv_dsig[3, s_ind]
    end

    # Hessian terms involving both x and sigma.
    # Note that dpyA / dxB = bvn.precision[A, B]
    bvn_xsig_h = bvn_derivs.bvn_xsig_h
    for x_ind=1:2
      bvn_xsig_h[x_ind, 1] = bvn_derivs.py1[1] * bvn.precision[1, x_ind]
      bvn_xsig_h[x_ind, 2] =
        bvn_derivs.py1[1] * bvn.precision[2, x_ind] +
        bvn_derivs.py2[1] * bvn.precision[1, x_ind]
      bvn_xsig_h[x_ind, 3] = bvn_derivs.py2[1] * bvn.precision[2, x_ind]
    end
  end
end

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

"""
The derivatives of sigma with respect to the galaxy shape parameters.  In
each case, sigma is [Sigma11, Sigma12, Sigma22], and the galaxy shape
parameters are indexed by GalaxyShapeParams.
 - j: A Jacobian matrix Sigma x GalaxyShapeParams of
      d Sigma / d GalaxyShapeParams
 - t: A Sigma x GalaxyShapeParams x GalaxyShapeParams tensor of second
      derivatives d2 Sigma / d GalaxyShapeParams d GalaxyShapeParams.
"""
type GalaxySigmaDerivs{NumType <: Number}
  j::Matrix{NumType}
  t::Array{NumType, 3}
end

"""
Args:
  - e_angle: Phi in the notes
  - e_axis: Rho in the notes
  - e_scale: Lower case sigma in the notes
  - XiXi: The value of sigma.

Note that nubar is not included.
"""
function GalaxySigmaDerivs{NumType <: Number}(
    e_angle::NumType, e_axis::NumType, e_scale::NumType,
    XiXi::Matrix{NumType}; calculate_tensor::Bool=true)

  cos_sin = cos(e_angle)sin(e_angle)
  sin_sq = sin(e_angle)^2
  cos_sq = cos(e_angle)^2

  j = Array(NumType, 3, length(gal_shape_ids))
  for i = 1:3
    j[i, gal_shape_ids.e_axis] =
      2 * e_axis * e_scale^2 * [sin_sq, -cos_sin, cos_sq][i]
    j[i, gal_shape_ids.e_angle] =
      e_scale^2 * (e_axis^2 - 1) * [2cos_sin, sin_sq - cos_sq, -2cos_sin][i]
    j[i, gal_shape_ids.e_scale] = (2XiXi ./ e_scale)[[1, 2, 4][i]]
  end

  t = Array(NumType, 3, length(gal_shape_ids), length(gal_shape_ids))
  if calculate_tensor
    # Second derivatives.

    for i = 1:3
      # Second derivatives involving e_scale
      t[i, gal_shape_ids.e_scale, gal_shape_ids.e_scale] =
        (2 * XiXi ./ (e_scale ^ 2))[[1, 2, 4][i]]
      t[i, gal_shape_ids.e_scale, gal_shape_ids.e_axis] =
        (2 * j[i, gal_shape_ids.e_axis] ./ e_scale)
      t[i, gal_shape_ids.e_scale, gal_shape_ids.e_angle] =
        (2 * j[i, gal_shape_ids.e_angle] ./ e_scale)

      t[i, gal_shape_ids.e_axis, gal_shape_ids.e_scale] =
        t[i, gal_shape_ids.e_scale, gal_shape_ids.e_axis]
      t[i, gal_shape_ids.e_angle, gal_shape_ids.e_scale] =
        t[i, gal_shape_ids.e_scale, gal_shape_ids.e_angle]

      # Remaining second derivatives involving e_angle
      t[i, gal_shape_ids.e_angle, gal_shape_ids.e_angle] =
        2 * e_scale^2 * (e_axis^2 - 1) *
        [cos_sq - sin_sq, 2cos_sin, sin_sq - cos_sq][i]
      t[i, gal_shape_ids.e_angle, gal_shape_ids.e_axis] =
        2 * e_scale^2 * e_axis * [2cos_sin, sin_sq - cos_sq, -2cos_sin][i]
      t[i, gal_shape_ids.e_axis, gal_shape_ids.e_angle] =
        t[i, gal_shape_ids.e_angle, gal_shape_ids.e_axis]

      # The second derivative involving only e_axis.
      t[i, gal_shape_ids.e_axis, gal_shape_ids.e_axis] =
        2 * e_scale^2 * [sin_sq, -cos_sin, cos_sq][i]
    end
  else
    fill!(t, 0.0)
  end

  GalaxySigmaDerivs(j, t)
end

"""
The convolution of a one galaxy component with one PSF component.
It also contains the derivatives of sigma with respect to the shape parameters.
It does not contain the derivatives with respect to other parameters
(u and e_dev) because they have easy expressions in terms of other known
quantities.

Args:
 - e_dev_dir: "Theta direction": this is 1 or -1, depending on whether
     increasing e_dev increases the weight of this GalaxyCacheComponent
     (1) or decreases it (-1).
 - e_dev_i: The weight given to this type of galaxy for this celestial object.
     This is either e_dev or (1 - e_dev).
 - gc: The galaxy component to be convolved
 - pc: The psf component to be convolved
 - u: The location of the celestial object in pixel coordinates as a 2x1 vector
 - e_axis: The ratio of the galaxy minor axis to major axis (0 < e_axis <= 1)
 - e_scale: The scale of the galaxy major axis

Attributes:
 - e_dev_dir: Same as input
 - e_dev_i: Same as input
 - bmc: A BvnComponent with the convolution.
 - dSigma: A 3x3 matrix containing the derivates of
     [Sigma11, Sigma12, Sigma22] (in the rows) with respect to
     [e_axis, e_angle, e_scale] (in the columns)
"""
immutable GalaxyCacheComponent{NumType <: Number}
    e_dev_dir::Float64
    e_dev_i::NumType
    bmc::BvnComponent{NumType}
    sig_sf::GalaxySigmaDerivs{NumType}
    # [Sigma11, Sigma12, Sigma22] x [e_axis, e_angle, e_scale]
end

function GalaxyCacheComponent{NumType <: Number}(
    e_dev_dir::Float64, e_dev_i::NumType,
    gc::GalaxyComponent, pc::PsfComponent, u::Vector{NumType},
    e_axis::NumType, e_angle::NumType, e_scale::NumType,
    calculate_derivs::Bool, calculate_hessian::Bool)

  # Declare in advance to save memory allocation.
  const empty_sig_sf =
    GalaxySigmaDerivs(Array(NumType, 0, 0), Array(NumType, 0, 0, 0))

  XiXi = get_bvn_cov(e_axis, e_angle, e_scale)
  mean_s = NumType[pc.xiBar[1] + u[1], pc.xiBar[2] + u[2]]
  var_s = pc.tauBar + gc.nuBar * XiXi
  weight = pc.alphaBar * gc.etaBar  # excludes e_dev

  # d siginv / dsigma is only necessary for the Hessian.
  bmc = BvnComponent{NumType}(
    mean_s, var_s, weight,
    calculate_siginv_deriv=calculate_derivs && calculate_hessian)

  if calculate_derivs
    sig_sf = GalaxySigmaDerivs(
      e_angle, e_axis, e_scale, XiXi, calculate_tensor=calculate_hessian)
    sig_sf.j .*= gc.nuBar
    if calculate_hessian
      # The tensor is only needed for the Hessian.
      sig_sf.t .*= gc.nuBar
    end
  else
    sig_sf = empty_sig_sf
  end

  GalaxyCacheComponent(e_dev_dir, e_dev_i, bmc, sig_sf)
end

###################################################
# Transform derivatives into the model parameterization.

"""
Transform the bvn derivatives and hessians from (x) to the
galaxy parameters (u). Updates bvn_u_d and bvn_uu_h in place.

bvn_x_d and bvn_xx_h should already have been set using get_bvn_derivs!()
"""
function transform_bvn_ux_derivs!{NumType <: Number}(
    bvn_derivs::BivariateNormalDerivatives{NumType},
    wcs_jacobian::Array{Float64, 2}, calculate_hessian::Bool)

  # Gradient calculations.

  # Note that dxA_duB = -wcs_jacobian[A, B].  (It is minus the jacobian
  # because the object position affects the bvn.the_mean term, which is
  # subtracted from the pixel location as defined in bvn_sf.d.)
  bvn_u_d = bvn_derivs.bvn_u_d
  bvn_x_d = bvn_derivs.bvn_x_d
  bvn_u_d[1] =
    -(bvn_x_d[1] * wcs_jacobian[1, 1] + bvn_x_d[2] * wcs_jacobian[2, 1])
  bvn_u_d[2] =
    -(bvn_x_d[1] * wcs_jacobian[1, 2] + bvn_x_d[2] * wcs_jacobian[2, 2])

  if calculate_hessian
    # Hessian calculations.

    bvn_uu_h = bvn_derivs.bvn_uu_h
    bvn_xx_h = bvn_derivs.bvn_xx_h
    fill!(bvn_uu_h, 0.0)
    # Second derivatives involving only u.
    # As above, dxA_duB = -wcs_jacobian[A, B] and d2x / du2 = 0.
    # TODO: time consuming **************
    @inbounds for x_id2 in 1:2, x_id1 in 1:2, u_id2 in 1:2
      inner_term = bvn_xx_h[x_id1, x_id2] * wcs_jacobian[x_id2, u_id2]
      @inbounds for u_id1 in 1:u_id2
        bvn_uu_h[u_id1, u_id2] += inner_term * wcs_jacobian[x_id1, u_id1]
      end
    end
    @inbounds bvn_uu_h[2, 1] = bvn_uu_h[1, 2]
  end
end

"""
Transform all the bvn derivatives from x and sigma to u and the model
parameters.

These values should already have been set using get_bvn_derivs!():
bvn_x_d
bvn_xx_h
bvn_sig_d
bvn_sigsig_h
bvn_xsig_h
"""
function transform_bvn_derivs!{NumType <: Number}(
    bvn_derivs::BivariateNormalDerivatives{NumType},
    sig_sf::GalaxySigmaDerivs{NumType},
    wcs_jacobian::Array{Float64, 2}, calculate_hessian::Bool)

  # Transform the u derivates first.
  # bvn_x_d and bvn_xx_h should already have been set using get_bvn_derivs!()
  transform_bvn_ux_derivs!(bvn_derivs, wcs_jacobian, calculate_hessian)

  # Gradient calculations.

  # Use the chain rule for the shape derviatives.
  # TODO: time consuming **************
  bvn_s_d = bvn_derivs.bvn_s_d
  bvn_sig_d = bvn_derivs.bvn_sig_d

  fill!(bvn_s_d, 0.0)
  @inbounds for shape_id in 1:length(gal_shape_ids), sig_id in 1:3
    bvn_s_d[shape_id] += bvn_sig_d[sig_id] * sig_sf.j[sig_id, shape_id]
  end

  if calculate_hessian
    # Hessian calculations.

    bvn_ss_h = bvn_derivs.bvn_ss_h
    bvn_us_h = bvn_derivs.bvn_us_h

    fill!(bvn_ss_h, 0.0)
    fill!(bvn_us_h, 0.0)

    # Second derviatives involving only shape parameters.
    # TODO: time consuming **************
    @inbounds for shape_id2 in 1:length(gal_shape_ids), shape_id1 in 1:shape_id2
      @inbounds for sig_id1 in 1:3
        bvn_ss_h[shape_id1, shape_id2] +=
          bvn_sig_d[sig_id1] * sig_sf.t[sig_id1, shape_id1, shape_id2]
      end
    end

    bvn_sigsig_h = bvn_derivs.bvn_sigsig_h
    @inbounds for sig_id1 in 1:3, sig_id2 in 1:3,
                  shape_id2 in 1:length(gal_shape_ids)
      inner_term =
        bvn_sigsig_h[sig_id1, sig_id2] * sig_sf.j[sig_id2, shape_id2]
      @inbounds for shape_id1 in 1:shape_id2
        bvn_ss_h[shape_id1, shape_id2] +=
          inner_term * sig_sf.j[sig_id1, shape_id1]
      end
    end

    @inbounds for shape_id2 in 1:length(gal_shape_ids), shape_id1 in 1:shape_id2
      bvn_ss_h[shape_id2, shape_id1] = bvn_ss_h[shape_id1, shape_id2]
    end

    # Second derivates involving both a shape term and a u term.
    # TODO: time consuming **************
    bvn_xsig_h = bvn_derivs.bvn_xsig_h
    @inbounds for shape_id in 1:length(gal_shape_ids),
                  u_id in 1:2, sig_id in 1:3, x_id in 1:2
      bvn_us_h[u_id, shape_id] +=
        bvn_xsig_h[x_id, sig_id] * sig_sf.j[sig_id, shape_id] *
        (-wcs_jacobian[x_id, u_id])
    end
  end
end


text:
using Test
using ScottishTaxBenefitModel
using .Randoms
using .ExampleHelpers

@testset "Randoms" begin
    rm  = "123456789"
    @test randchunk( rm, 1 ) ≈ 0.12345678 # default 8 digits
    @test randchunk( rm, 1, 9 ) ≈ 0.123456789
    @test testp( rm, 0.2346, 2 ) # random from string is less than thresh
    @test ! testp( rm, 0.2344, 2 ) # random from string is > thresh
    r1 = strtobi(mybigrand())
    print( "r1 = $r1")
    c = randchunk( r1, 1 )
    print( "c=$c")
    b = testp( r1, 0.2, 1 )
    print( "b=$b" )
    if c > 0.2
        @test ! b
    else
        @test b
    end

    t = 0
    n = 5_000_000
    for i in 1:n
        r = strtobi(mybigrand())
        if testp( r, 0.2, 10 )
            t += 1
        end
    end
    pr = t/n
    @test isapprox( pr, 0.2, atol=0.001)

end

ext: jl

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
begin
    using Luxor, Random, Colors

    Drawing(512, 512, "hello-world.png")
    background("black")
    sethue("white")
    translate(50, 50)
    radius = 30
    grid = GridHex(O, radius, 2)
    t = Table(10, 1)
    for (pt, n) in t
        randomhue()
        setgray(rescale(n, 1, length(t), 0, 1))
        p = nextgridpoint(grid)
        ngon(p, radius - 5, 6, π / 2, :fillstroke)
        sethue("white")
        fontsize(25)
        text(string(n), p, halign = :center)
    end

    translate(100, 0)
    radius = 30
    grid = GridHex(O, radius, 2)
    t = Table(10, 1)
    list = ["A" "B" "C" "D" "E" "F" "G" "H" "I" "J"]
    for (pt, n) in t
        randomhue()
        setgray(rescale(n, 1, length(t), 0, 1))
        p = nextgridpoint(grid)
        ngon(p, radius - 5, 6, π / 2, :fillstroke)
        sethue("white")
        fontsize(25)
        text(string(list[n]), p, halign = :center)
    end
    origin()
    translate(40, -150)

    circle(Point(0, 0), 90, :stroke)
    setdash("dot")
    circle(Point(0, 0), 80, :stroke)
    circle(Point(0, 0), 50, :fill)
    circlepath(Point(0, 0), 70, :path)
    newsubpath()
    circlepath(Point(0, 0), 60, :fill, reversepath = true)

    origin()
    translate(0, 70)
    setline(0.5)
    apolygon = star(O, 50, 5, 0.5, 0, vertices = true)
    for n = 1:500
        apoint = randompoint(Point(-50, -50), Point(50, 50))
        isinside(apoint, apolygon) ? circle(apoint, 3, :fill) : circle(apoint, 0.5, :stroke)
    end

    origin()
    sethue("white")
    translate(100, 70)
    box(Point(0, 0), 70, 70, :fill)
    box(Point(0, 0), 90, 90, :stroke)

    origin()
    translate(-50, -30)
    sethue("white")
    fontsize(30)
    text("HELLOWORLD!!")
    translate(0, 40)
    text("1 2 3 4 5 6 7 8 9 0")
    setdash("solid")
    origin()
    translate(150, -150)
    for i = 1:10
        setline(i)
        x = LinRange(1, 100, 10)
        rule(Point(x[i], 0), π / 2, boundingbox = BoundingBox() * 0.3)
    end

    origin()
    translate(200, 375)
    fib = [1, 1, 2, 3, 5, 8, 13, 21, 34]
    tiles = boxmap(fib, BoundingBox()[1], 200, 125)

    for (n, t) in enumerate(tiles)
        randomhue()
        bb = BoundingBox(t)
        sethue("black")
        box(bb - 3, :stroke)
        randomhue()
        box(bb - 5, :fill)

        sethue("white")
       fontsize(boxwidth(bb) > boxheight(bb) ? (boxheight(bb)/4)+20 : (boxwidth(bb)/4)+20)
		text(
            string(sort(fib, rev = true)[n]),
            midpoint(bb[1], bb[2]),
            halign = :center,
            valign = :middle,
        )
    end

    origin()
    translate(150, 100)
    scale(0.3)
    sethue("white")
    box(Point(160, 100), 350, 300, :fill)
    julialogo()

    finish()
    preview()

end

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
"""
Package functions for interacting with Travis.

$(EXPORTS)
"""
module Travis

using Compat, DocStringExtensions
import Compat.Pkg

export genkeys

import Compat.LibGit2.GITHUB_REGEX

"""
$(SIGNATURES)

Generate ssh keys for package `package` to automatically deploy docs from Travis to GitHub
pages. `package` can be either the name of a package or a path. Providing a path allows keys
to be generated for non-packages or packages that are not found in the Julia `LOAD_PATH`.
Use the `remote` keyword to specify the user and repository values.

This function requires the following command lines programs to be installed:

- `which`
- `git`
- `travis`
- `ssh-keygen`

# Examples

```jlcon
julia> using Documenter

julia> Travis.genkeys("MyPackageName")
[ ... output ... ]

julia> Travis.genkeys("MyPackageName", remote="organization")
[ ... output ... ]

julia> Travis.genkeys("/path/to/target/directory")
[ ... output ... ]
```
"""
function genkeys(package; remote="origin")
    # Error checking. Do the required programs exist?
    success(`which which`)      || error("'which' not found.")
    success(`which git`)        || error("'git' not found.")
    success(`which ssh-keygen`) || error("'ssh-keygen' not found.")

    directory = "docs"
    filename  = ".documenter"

    path = isdir(package) ? package : Pkg.dir(package, directory)
    isdir(path) || error("`$path` not found. Provide a package name or directory.")

    cd(path) do
        # Check for old '$filename.enc' and terminate.
        isfile("$filename.enc") &&
            error("$package already has an ssh key. Remove it and try again.")

        # Are we in a git repo?
        success(`git status`) || error("'Travis.genkey' only works with git repositories.")

        # Find the GitHub repo org and name.
        user, repo =
            let r = readchomp(`git config --get remote.$remote.url`)
                m = match(GITHUB_REGEX, r)
                m === nothing && error("no remote repo named '$remote' found.")
                m[2], m[3]
            end

        # Generate the ssh key pair.
        success(`ssh-keygen -N "" -f $filename`) || error("failed to generated ssh key pair.")

        # Prompt user to add public key to github then remove the public key.
        let url = "https://github.com/$user/$repo/settings/keys"
            Compat.@info("add the public key below to $url with read/write access:")
            println("\n", read("$filename.pub", String))
            rm("$filename.pub")
        end

        # Base64 encode the private key and prompt user to add it to travis. The key is
        # *not* encoded for the sake of security, but instead to make it easier to
        # copy/paste it over to travis without having to worry about whitespace.
        let url = "https://travis-ci.org/$user/$repo/settings"
            Compat.@info("add a secure environment variable named 'DOCUMENTER_KEY' to $url with value:")
            println("\n", base64encode(read(".documenter", String)), "\n")
            rm(filename)
        end
    end
end

end # module


text:
"""
    macro stack(into, structs...)

Stack the fields of multiple structs and create a new one.
The first argument is the name of the new struct followed by the
ones to be stacked. Parametric types are not supported and
the fieldnames needs to be unique.

Example:

    @stack Baz Foo Bar

Creates `Baz` with the concatenated fields of `Foo` and `Bar`

"""
macro stack(into, structs...)
    fields = []
    for _struct in structs
        names = fieldnames(getfield(__module__, _struct))
        types = [fieldtype(getfield(__module__, _struct), field) for field in names]
        for (n, t) in zip(names, types)
            push!(fields, :($n::$t))
        end
    end

    esc(
        quote
            struct $into
                $(fields...)
            end
        end
    )
end

"""
    unpack(x::CompressionHeader)

Return the following information:
- Name of compression algorithm
- Level of the compression
- compressedbytes and uncompressedbytes according to [uproot3](https://github.com/scikit-hep/uproot3/blob/54f5151fb7c686c3a161fbe44b9f299e482f346b/uproot3/source/compressed.py#L132)
"""
function unpack(x::CompressionHeader)
    algname = String(x.algo)
    ver = Int(x.method)
    # shift without casting to `Int` will give you 0x00 because we're shifting 0 bits into UInt8
    compressedbytes = x.c1 + (Int(x.c2) << 8) + (Int(x.c3) << 16)
    uncompressedbytes = x.u1 + (Int(x.u2) << 8) + (Int(x.u3) << 16)

    return algname, ver, compressedbytes, uncompressedbytes
end


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

abstract type AbstractWorkerPool end

# An AbstractWorkerPool should implement
#
# `push!` - add a new worker to the overall pool (available + busy)
# `put!` - put back a worker to the available pool
# `take!` - take a worker from the available pool (to be used for remote function execution)
# `length` - number of workers available in the overall pool
# `isready` - return false if a `take!` on the pool would block, else true
#
# The default implementations of the above (on a AbstractWorkerPool) require fields
#    channel::Channel{Int}
#    workers::Set{Int}
#

mutable struct WorkerPool <: AbstractWorkerPool
    channel::Channel{Int}
    workers::Set{Int}
    ref::RemoteChannel

    WorkerPool(c::Channel, ref::RemoteChannel) = new(c, Set{Int}(), ref)
end

function WorkerPool()
    wp = WorkerPool(Channel{Int}(typemax(Int)), RemoteChannel())
    put!(wp.ref, WeakRef(wp))
    wp
end

"""
    WorkerPool(workers::Vector{Int})

Create a WorkerPool from a vector of worker ids.
"""
function WorkerPool(workers::Vector{Int})
    pool = WorkerPool()
    foreach(w->push!(pool, w), workers)
    return pool
end

# On workers where this pool has been serialized to, instantiate with a dummy local channel.
WorkerPool(ref::RemoteChannel) = WorkerPool(Channel{Int}(1), ref)

function serialize(S::AbstractSerializer, pool::WorkerPool)
    # Allow accessing a worker pool from other processors. When serialized,
    # initialize the `ref` to point to self and only send the ref.
    # Other workers will forward all put!, take!, calls to the process owning
    # the ref (and hence the pool).
    Serialization.serialize_type(S, typeof(pool))
    serialize(S, pool.ref)
end

deserialize(S::AbstractSerializer, t::Type{T}) where {T<:WorkerPool} = T(deserialize(S))

wp_local_push!(pool::AbstractWorkerPool, w::Int) = (push!(pool.workers, w); put!(pool.channel, w); pool)
wp_local_length(pool::AbstractWorkerPool) = length(pool.workers)
wp_local_isready(pool::AbstractWorkerPool) = isready(pool.channel)

function wp_local_put!(pool::AbstractWorkerPool, w::Int)
    # In case of default_worker_pool, the master is implictly considered a worker, i.e.,
    # it is not present in pool.workers.
    # Confirm the that the worker is part of a pool before making it available.
    w in pool.workers && put!(pool.channel, w)
    w
end

function wp_local_workers(pool::AbstractWorkerPool)
    if length(pool) == 0 && pool === default_worker_pool()
        return [1]
    else
        return collect(pool.workers)
    end
end

function wp_local_nworkers(pool::AbstractWorkerPool)
    if length(pool) == 0 && pool === default_worker_pool()
        return 1
    else
        return length(pool.workers)
    end
end

function wp_local_take!(pool::AbstractWorkerPool)
    # Find an active worker
    worker = 0
    while true
        if length(pool) == 0
            if pool === default_worker_pool()
                # No workers, the master process is used as a worker
                worker = 1
                break
            else
                throw(ErrorException("No active worker available in pool"))
            end
        end

        worker = take!(pool.channel)
        if id_in_procs(worker)
            break
        else
            delete!(pool.workers, worker) # Remove invalid worker from pool
        end
    end
    return worker
end

function remotecall_pool(rc_f, f, pool::AbstractWorkerPool, args...; kwargs...)
    worker = take!(pool)
    try
        rc_f(f, worker, args...; kwargs...)
    finally
        put!(pool, worker)
    end
end

# Check if pool is local or remote and forward calls if required.
# NOTE: remotecall_fetch does it automatically, but this will be more efficient as
# it avoids the overhead associated with a local remotecall.

for func = (:length, :isready, :workers, :nworkers, :take!)
    func_local = Symbol(string("wp_local_", func))
    @eval begin
        function ($func)(pool::WorkerPool)
            if pool.ref.where != myid()
                return remotecall_fetch(ref->($func_local)(fetch(ref).value), pool.ref.where, pool.ref)
            else
                return ($func_local)(pool)
            end
        end

        # default impl
        ($func)(pool::AbstractWorkerPool) = ($func_local)(pool)
    end
end

for func = (:push!, :put!)
    func_local = Symbol(string("wp_local_", func))
    @eval begin
        function ($func)(pool::WorkerPool, w::Int)
            if pool.ref.where != myid()
                return remotecall_fetch((ref, w)->($func_local)(fetch(ref).value, w), pool.ref.where, pool.ref, w)
            else
                return ($func_local)(pool, w)
            end
        end

        # default impl
        ($func)(pool::AbstractWorkerPool, w::Int) = ($func_local)(pool, w)
    end
end

"""
    remotecall(f, pool::AbstractWorkerPool, args...; kwargs...) -> Future

`WorkerPool` variant of `remotecall(f, pid, ....)`. Waits for and takes a free worker from `pool` and performs a `remotecall` on it.
"""
remotecall(f, pool::AbstractWorkerPool, args...; kwargs...) = remotecall_pool(remotecall, f, pool, args...; kwargs...)

"""
    remotecall_wait(f, pool::AbstractWorkerPool, args...; kwargs...) -> Future

`WorkerPool` variant of `remotecall_wait(f, pid, ....)`. Waits for and takes a free worker from `pool` and
performs a `remotecall_wait` on it.
"""
remotecall_wait(f, pool::AbstractWorkerPool, args...; kwargs...) = remotecall_pool(remotecall_wait, f, pool, args...; kwargs...)

"""
    remotecall_fetch(f, pool::AbstractWorkerPool, args...; kwargs...) -> result

`WorkerPool` variant of `remotecall_fetch(f, pid, ....)`. Waits for and takes a free worker from `pool` and
performs a `remotecall_fetch` on it.
"""
remotecall_fetch(f, pool::AbstractWorkerPool, args...; kwargs...) = remotecall_pool(remotecall_fetch, f, pool, args...; kwargs...)

"""
    remote_do(f, pool::AbstractWorkerPool, args...; kwargs...) -> nothing

`WorkerPool` variant of `remote_do(f, pid, ....)`. Waits for and takes a free worker from `pool` and
performs a `remote_do` on it.
"""
remote_do(f, pool::AbstractWorkerPool, args...; kwargs...) = remotecall_pool(remote_do, f, pool, args...; kwargs...)

const _default_worker_pool = Ref{Union{WorkerPool, Nothing}}(nothing)

"""
    default_worker_pool()

`WorkerPool` containing idle `workers` - used by `remote(f)` and [`pmap`](@ref) (by default).
"""
function default_worker_pool()
    # On workers retrieve the default worker pool from the master when accessed
    # for the first time
    if _default_worker_pool[] === nothing
        if myid() == 1
            _default_worker_pool[] = WorkerPool()
        else
            _default_worker_pool[] = remotecall_fetch(()->default_worker_pool(), 1)
        end
    end
    return _default_worker_pool[]
end

"""
    remote([::AbstractWorkerPool], f) -> Function

Return an anonymous function that executes function `f` on an available worker
using [`remotecall_fetch`](@ref).
"""
remote(f) = (args...; kwargs...)->remotecall_fetch(f, default_worker_pool(), args...; kwargs...)
remote(p::AbstractWorkerPool, f) = (args...; kwargs...)->remotecall_fetch(f, p, args...; kwargs...)

mutable struct CachingPool <: AbstractWorkerPool
    channel::Channel{Int}
    workers::Set{Int}

    # Mapping between a tuple (worker_id, f) and a remote_ref
    map_obj2ref::Dict{Tuple{Int, Function}, RemoteChannel}

    function CachingPool()
        wp = new(Channel{Int}(typemax(Int)), Set{Int}(), Dict{Int, Function}())
        finalizer(clear!, wp)
        wp
    end
end

serialize(s::AbstractSerializer, cp::CachingPool) = throw(ErrorException("CachingPool objects are not serializable."))

"""
    CachingPool(workers::Vector{Int})

An implementation of an `AbstractWorkerPool`.
[`remote`](@ref), [`remotecall_fetch`](@ref),
[`pmap`](@ref) (and other remote calls which execute functions remotely)
benefit from caching the serialized/deserialized functions on the worker nodes,
especially closures (which may capture large amounts of data).

The remote cache is maintained for the lifetime of the returned `CachingPool` object.
To clear the cache earlier, use `clear!(pool)`.

For global variables, only the bindings are captured in a closure, not the data.
`let` blocks can be used to capture global data.

# Examples
```julia
const foo = rand(10^8);
wp = CachingPool(workers())
let foo = foo
    pmap(wp, i -> sum(foo) + i, 1:100);
end
```

The above would transfer `foo` only once to each worker.

"""
function CachingPool(workers::Vector{Int})
    pool = CachingPool()
    for w in workers
        push!(pool, w)
    end
    return pool
end

"""
    clear!(pool::CachingPool) -> pool

Removes all cached functions from all participating workers.
"""
function clear!(pool::CachingPool)
    for (_,rr) in pool.map_obj2ref
        finalize(rr)
    end
    empty!(pool.map_obj2ref)
    pool
end

exec_from_cache(rr::RemoteChannel, args...; kwargs...) = fetch(rr)(args...; kwargs...)
function exec_from_cache(f_ref::Tuple{Function, RemoteChannel}, args...; kwargs...)
    put!(f_ref[2], f_ref[1])        # Cache locally
    f_ref[1](args...; kwargs...)
end

function remotecall_pool(rc_f, f, pool::CachingPool, args...; kwargs...)
    worker = take!(pool)
    f_ref = get(pool.map_obj2ref, (worker, f), (f, RemoteChannel(worker)))
    isa(f_ref, Tuple) && (pool.map_obj2ref[(worker, f)] = f_ref[2])   # Add to tracker

    try
        rc_f(exec_from_cache, worker, f_ref, args...; kwargs...)
    finally
        put!(pool, worker)
    end
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
using sem, Test, ModelingToolkit, LinearAlgebra, SparseArrays

function impfun(x)

    S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
        0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
        0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
        0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
        0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
        0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
        0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
        0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
        0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
        0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
        0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
        0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
        0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
        0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

    F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 1 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 1 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 1 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 1 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 1 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 1 0 0 0]

    A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
        0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
        0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
        0  0  0  0  0  0  0  0  0  0  0     0     1     0
        0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
        0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
        0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
        0  0  0  0  0  0  0  0  0  0  0     0     0     1
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
        0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
        0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

    impcov = F*inv(I-A)*S*permutedims(inv(I-A))*permutedims(F)

    return impcov
end

@variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

start_val = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3)
    )

new_val = fill(1.0, 31)

new_val_2 = rand(31)

##Symbolic

implysym_test = ImplySymbolic(A, S, F, x, start_val)

implysym_test(start_val)

@test implysym_test.imp_cov == impfun(start_val)

implysym_test(new_val)

@test implysym_test.imp_cov == impfun(new_val)

implysym_test(new_val_2)

@test implysym_test.imp_cov ≈ impfun(new_val_2)

invia = I + A
next_term = A^2

while nnz(next_term) != 0
    global invia += next_term
    global next_term *= A
end

ModelingToolkit.simplify.(invia)

invia2 = I + A + A^2 + A^3

ModelingToolkit.simplify.(invia2)

imp_cov_sym = F*invia2*S*permutedims(invia2)*permutedims(F)

imp_cov_sym = Array(imp_cov_sym)
imp_cov_sym .= ModelingToolkit.simplify.(imp_cov_sym)

imp_fun =
    eval(ModelingToolkit.build_function(
        imp_cov_sym,
        x
    )[2])

mat = rand(11,11)

imp_fun(mat, start_val)

imp_fun(mat, new_val_2)

mat ≈ impfun(new_val_2)

##Sparse
implyspa_test = ImplySparse(A, S, F, x, start_val)

implyspa_test(start_val)

@test implyspa_test.imp_cov == impfun(start_val)

implyspa_test(new_val)

@test implyspa_test.imp_cov .== impfun(new_val)

implyspa_test(new_val_2)

@test implyspa_test.imp_cov .== impfun(new_val_2)


text:
using Documenter
using DataStructures

DocMeta.setdocmeta!(DataStructures, :DocTestSetup, :(using DataStructures); recursive=true)

makedocs(
    sitename = "DataStructures.jl",
    pages = [
        "index.md",
        "deque.md",
        "circ_buffer.md",
        "circ_deque.md",
        "stack_and_queue.md",
        "priority-queue.md",
        "fenwick.md",
        "accumulators.md",
        "disjoint_set.md",
        "heaps.md",
        "ordered_containers.md",
        "default_dict.md",
        "robin_dict.md",
        "swiss_dict.md",
        "trie.md",
        "linked_list.md",
        "mutable_linked_list.md",
        "intset.md",
        "sorted_containers.md",
        "dibit_vector.md",
        "red_black_tree.md",
        "avl_tree.md",
        "splay_tree.md",
    ],
    modules = [DataStructures],
    format = Documenter.HTML()
)

deploydocs(
    repo = "github.com/JuliaCollections/DataStructures.jl.git",
)


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id="first", value = 1),
    dcc_input(id="second"),
    dcc_input(id="third")
end

callback!(
    ClientsideFunction("clientside", "add1_break_at_11"),
    app, callid"first.value=>second.value"
    )

callback!(
    ClientsideFunction("clientside", "add1_break_at_11"),
    app, callid"second.value=>third.value"
    )

run_server(app)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

mutable struct IterateChannel{A <: AbstractArray} <: AbstractChannel{A}
    decisions::Dict{Int,A}
    cond_take::Condition
    IterateChannel(decisions::Dict{Int,A}) where A <: AbstractArray = new{A}(decisions, Condition())
end

function put!(channel::IterateChannel, t, x)
    channel.decisions[t] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::IterateChannel, t)
    x = fetch(channel, t)
    delete!(channel.decisions, t)
    return x
end

isready(channel::IterateChannel) = length(channel.decisions) > 1
isready(channel::IterateChannel, t) = haskey(channel.decisions, t)

function fetch(channel::IterateChannel, t)
    wait(channel, t)
    return channel.decisions[t]
end

function wait(channel::IterateChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

RemoteIterates{A} = RemoteChannel{IterateChannel{A}}

const MetaData = Dict{Tuple{Int,Symbol},Any}

has_metadata(data::MetaData, type::Symbol) = haskey(data, (0, type))
has_metadata(data::MetaData, idx::Integer, type::Symbol) = haskey(data, (idx, type))
get_metadata(data::MetaData, type::Symbol) = data[(0, type)]
get_metadata(data::MetaData, idx::Integer, type::Symbol) = data[(idx, type)]
function set_metadata!(data::MetaData, type::Symbol, value)
    data[(0, type)] = value
    return nothing
end
function set_metadata!(data::MetaData, idx::Integer, type::Symbol, value)
    data[(idx, type)] = value
    return nothing
end

mutable struct MetaChannel <: AbstractChannel{Any}
    metadata::MetaData
    cond_take::Condition
    MetaChannel() = new(Dict{Tuple{Int,Symbol},Any}(), Condition())
end

function put!(channel::MetaChannel, t, key, x)
    channel.metadata[(t,key)] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::MetaChannel, t, key)
    x = fetch(channel, t, key)
    delete!(channel.metadata, (t,k))
    return x
end

isready(channel::MetaChannel) = length(channel.metadata) > 1
isready(channel::MetaChannel, t, key) = haskey(channel.metadata, (t,key))

function fetch(channel::MetaChannel, t, key)
    wait(channel, t, key)
    return channel.metadata[(t,key)]
end

function wait(channel::MetaChannel, t, key)
    while !isready(channel, t, key)
        wait(channel.cond_take)
    end
end

MetaDataChannel = RemoteChannel{MetaChannel}

has_metadata(data::MetaDataChannel, type::Symbol) = isready(data, 0, type)
has_metadata(data::MetaDataChannel, idx::Integer, type::Symbol) = isready(data, idx, type)
get_metadata(data::MetaDataChannel, type::Symbol) = fetch(data, 0, type)
get_metadata(data::MetaDataChannel, idx::Integer, type::Symbol) = fetch(data, idx, type)
function set_metadata!(data::MetaDataChannel, type::Symbol, value)
    put!(data, 0, type, value)
    return nothing
end
function set_metadata!(data::MetaDataChannel, idx::Integer, type::Symbol, value)
    put!(data, idx, type, value)
    return nothing
end


text:
module TestAddKwonly
try
    using Test
catch
    using Base.Test
end
using Reconstructables: @add_kwonly, add_kwonly
import Reconstructables
include("utils.jl")

@add_kwonly function f(a, b; c=3, d=4)
  (a, b, c, d)
end
@test f(1, 2) == (1, 2, 3, 4)
@test f(a=1, b=2) == (1, 2, 3, 4)
@test_throws Reconstructables.UndefKeywordError f()

@add_kwonly g(a, b; c=3, d=4) = (a, b, c, d)
@test g(1, 2) == (1, 2, 3, 4)
@test g(a=1, b=2) == (1, 2, 3, 4)

@add_kwonly h(; c=3, d=4) = (c, d)
@test h() == (3, 4)

@add_kwonly with_type(a::A; b::B=2) where {A, B} = (a, A, b, B)
@test with_type(a=10, b=10.0) == (10, Int, 10.0, Float64)

@add_kwonly typed_kwonly(; a::A=1, b::B=2) where {A, B} = (a, A, b, B)
@test typed_kwonly(a=10, b=10.0) == (10, Int, 10.0, Float64)

@add_kwonly with_kwargs(a; b=2, kwargs...) =
    (a, b, Any[(k, v) for (k, v) in kwargs])
@test with_kwargs(a=10, x=20) == (10, 2, Any[(:x, 20)])

if VERSION >= v"0.7-"
    @eval @add_kwonly required_kwargs(a; b, c=3) = (a, b, c)
    @test required_kwargs(1, b=2) == (1, 2, 3)
    @test required_kwargs(a=1, b=2) == (1, 2, 3)
end

@test_error begin
    @eval @add_kwonly i(c=3, d=4) = (c, d)
end (err) -> begin
    @test err isa ErrorException
    @test contains(err.msg,
                   "At least one positional mandatory argument is required")
end

@test_error begin
    @eval @add_kwonly if false
        g(a, b; c=3, d=4) = (a, b, c, d)
    end
end (err) -> begin
    @test err isa ErrorException
    @test contains(err.msg,
                   "add_only does not work with expression if")
end

@test_error begin
    add_kwonly(:(f(x, f(x)) = x))
end (err) -> begin
    @test err isa ErrorException
    @test contains(err.msg,
                   "Not expecting to see: f(x)")
end

let io = IOBuffer()
    showerror(io, Reconstructables.UndefKeywordError(:X))
    msg = String(take!(copy(io)))
    @test contains(msg, "UndefKeywordError: keyword argument X not assigned")
end

end


text:
# Write your package code here.

using simulationQP934

data = Data();
box = initial_point(data);
nc, fatm, natm = linkedlist(box,data);
force_vec = fvec(data)

println(utotal(box, data, fatm, natm, nc))

force!(box, data, force_vec, fatm, natm, nc)

println(force_vec)

#using BenchmarkTools, Test

#@btime utotal($box, $data, $fatm, $natm, $nc)
#@btime utotal_parallel($box, $data, $fatm, $natm, $nc)
#@btime utotal_parallel2($box,$data,$fatm,$natm,$nc)


text:
function Base.summary(io::IO, A::AbstractDimArray{T,N}) where {T,N}
    if N == 0
        print(io, "0-dimensional ")
    elseif N == 1
        print(io, size(A, 1), "-element ")
    else
        print(io, join(size(A), "×"), " ")
    end
    print(io, string(nameof(typeof(A)), "{$T,$N}"))
end

function Base.show(io::IO, mime::MIME"text/plain", A::AbstractDimArray{T,N}) where {T,N}
    lines = 0
    summary(io, A)
    _printname(io, name(A))
    lines += _printdims(io, mime, dims(A))
    !(isempty(dims(A)) || isempty(refdims(A))) && println(io)
    lines += _printrefdims(io, mime, refdims(A))
    println(io)
    ds = displaysize(io)
    ioctx = IOContext(io, :displaysize => (ds[1] - lines, ds[2]))

    # Printing the array data is optional, subtypes can
    # show other things here instead.
    show_after(ioctx, mime, A)

    return nothing
end

function Base.show(io::IO, mime::MIME"text/plain", stack::AbstractDimStack)
    print(io, nameof(typeof(stack)))
    _printdims(io, mime, dims(stack))
    nlayers = length(keys(stack))
    layers_str = nlayers == 1 ? "layer" : "layers"
    printstyled(io, "\nand "; color=:light_black)
    print(io, "$nlayers $layers_str:\n")
    maxlen = reduce(max, map(length ∘ string, keys(stack)))
    for key in keys(stack)
        layer = stack[key]
        pkey = rpad(key, maxlen)
        printstyled(io, "  :$pkey", color=:yellow)
        print(io, string(" ", eltype(layer)))
        field_dims = DD.dims(layer)
        n_dims = length(field_dims)
        printstyled(io, " dims: "; color=:light_black)
        if n_dims > 0
            for (d, dim) in enumerate(field_dims)
                _show_dimname(io, dim)
                d != length(field_dims) && print(io, ", ")
            end
            print(io, " (")
            for (d, dim) in enumerate(field_dims)
                print(io, "$(length(dim))")
                d != length(field_dims) && print(io, '×')
            end
            print(io, ')')
        end
        print(io, '\n')
    end

    md = metadata(stack)
    if !(md isa NoMetadata)
        n_metadata = length(md)
        if n_metadata > 0
            printstyled(io, "\nwith "; color=:light_black)
            show(io, mime, md)
        end
    end

    # Show anything else subtypes want to append
    show_after(io, mime, stack)

    return nothing
end

function Base.show(io::IO, mime::MIME"text/plain", dim::Dimension)
    get(io, :compact, false) && return _show_compact(io::IO, dim)

    printstyled(io, nameof(typeof(dim)); color=_dimcolor(io))
    if dim isa Dim
        printstyled(io, "{"; color=:red)
        printstyled(io, string(":", name(dim)); color=:yellow)
        printstyled(io, "}"; color=:red)
    else
        if nameof(typeof(dim)) != name(dim)
            _printname(io, name(dim))
        end
    end
    print(io, ": ")
    printstyled(io, "\n  val: "; color=:light_black)
    _printdimindex(io, val(dim))

    if !(mode(dim) isa AutoMode)
        printstyled(io, "\n  mode: "; color=:light_black)
        show(io, mime, mode(dim))
    end
    if !(metadata(dim) isa NoMetadata)
        printstyled(io, "\n  metadata: "; color=:light_black)
        show(io, mime, metadata(dim))
    end
    println(io)
    show(io, mime, typeof(dim))
end

function Base.show(io::IO, mime::MIME"text/plain", metadata::Metadata{N}) where N
    print(io, "Metadata")
    if N !== Nothing
        print(io, "{")
        show(io, N)
        print(io, "}")
    end
    printstyled(io, " of "; color=:light_black)
    show(io, mime, val(metadata))
end
Base.show(io::IO, mime::MIME"text/plain", mode::IndexMode) = _printmode(io, mode)

function Base.show(io::IO, mime::MIME"text/plain", mode::AbstractSampled)
    _printmode(io, mode)
    print(io, ": ")
    _printorder(io, mode)
    print(io, " ")
    _printspan(io, mode)
    print(io, " ")
    _printsampling(io, mode)
end

function Base.show(io::IO, mime::MIME"text/plain", mode::AbstractCategorical)
    _printmode(io, mode)
    print(io, ": ")
    _printorder(io, mode)
end

# Semi-interface methods for adding addional `show` text
# for AbstractDimArray/AbstractDimStack subtypes
# TODO actually document in the interface
show_after(io, mime, stack::DimStack) = nothing
show_after(io::IO, mime, A::AbstractDimArray) = show_array(io, mime, parent(A))

# Showing the array is optional for AbstractDimArray
# `show_array` must be called from `show_after`.
function show_array(io::IO, mime, A::AbstractArray{T,0}) where T
    print(_ioctx(io, T), "\n", A[])
end
function show_array(io::IO, mime, A::AbstractArray{T,1}) where T
    Base.print_matrix(_ioctx(io, T), A)
end
function show_array(io::IO, mime, A::AbstractArray{T,2}) where T
    Base.print_matrix(_ioctx(io, T), A)
end
function show_array(io::IO, mime, A::AbstractArray{T,N}) where {T,N}
    o = ones(Int, N-2)
    frame = A[:, :, o...]
    onestring = join(o, ", ")
    println(io, "[:, :, $(onestring)]")
    Base.print_matrix(_ioctx(io, T), frame)
    nremaining = prod(size(A, d) for d=3:N) - 1
    nremaining > 0 && print(io, "\n[and ", nremaining," more slices...]")
end

function _show_dimname(io, dim::Dim)
    color = DD._dimcolor(io)
    printstyled(io, "Dim{"; color=color)
    printstyled(io, string(":", name(dim)); color=:yellow)
    printstyled(io, "}"; color=color)
end
function _show_dimname(io, dim::Dimension)
    printstyled(io, DD.dim2key(dim); color = DD._dimcolor(io))
end

# short printing version for dimensions
function _show_compact(io::IO, dim::Dimension)
    nm = nameof(typeof(dim))
    nchars = length(string(nm))
    printstyled(io, nm; color=_dimcolor(io))
    if name(dim) != nm
        print(io, " (")
        print(io, name(dim))
        print(io, ")")
        nchars += length(string(name(dim))) + 3
    end
    _printdimproperties(io, dim, nchars)
end
function _show_compact(io::IO, dim::Dim)
    color = _dimcolor(io)
    printstyled(io, "Dim{"; color=color)
    printstyled(io, string(":", name(dim)); color=:yellow)
    printstyled(io, "}"; color=color)
    nchars = 5 + length(string(name(dim)))
    _printdimproperties(io, dim, nchars)
end

_dimcolor(io) = get(io, :is_ref_dim, false) ? :magenta : :red

function _printname(io::IO, name)
    if !(name == Symbol("") || name isa NoName)
        printstyled(io, string(" :", name); color=:yellow)
    end
end

function _printdims(io::IO, mime, dims::Tuple)
    if isempty(dims)
        print(io, ": ")
        return 0
    end
    printstyled(io, " with dimensions: "; color=:light_black)
    return _layout_dims(io, mime, dims)
end

function _printrefdims(io::IO, mime, refdims::Tuple)
    if isempty(refdims)
        return 0
    end
    printstyled(io, "and reference dimensions: "; color=:light_black)
    ctx = IOContext(io, :is_ref_dim=>true, :show_dim_val=>true)
    lines = _layout_dims(ctx, mime, refdims)
    return lines
end

function _layout_dims(io, mime, dims::Tuple)
    length(dims) > 0 || return 0
    ctx = IOContext(io, :compact=>true)
    if all(m -> m isa NoIndex, mode(dims))
        for d in dims[1:end-1]
            show(ctx, mime, d)
            print(io, ", ")
        end
        show(ctx, mime, dims[end])
        return 0
    else # Dims get a line each
        lines = 3
        print(io, "\n  ")
        for d in dims[1:end-1]
            show(ctx, mime, d)
            print(io, ",")
            lines += 2 # Often they wrap
            print(io, "\n  ")
        end
        show(ctx, mime, dims[end])
        return lines
    end
end

function _printdimproperties(io, dim::Dimension, nchars)
    mode(dim) isa NoIndex && return nothing
    iobuf = IOBuffer()
    print(io, ": ")
    nchars += 1
    _printdimindex(iobuf, val(dim))
    dimvalstr = String(take!(iobuf))
    printstyled(io, dimvalstr; color=:cyan)
    nchars += length(dimvalstr)
    if !(mode(dim) isa AutoMode)
        show(iobuf, MIME"text/plain"(), mode(dim))
        termwidth = displaysize(stdout)[2]
        modestr = String(take!(iobuf))
        if nchars + length(modestr) + 6 > termwidth
            print(io, "\n   ")
        end
        print(io, string(" ", modestr))
    end
    return nothing
end

_printdimindex(io, A) = printstyled(io, A; color=:cyan)
_printdimindex(io, A::AbstractRange) = printstyled(io, A; color=:cyan)
function _printdimindex(io, v::AbstractVector)
    s = string(eltype(v)) * "["
    svals =
    if length(v) > 2 && eltype(v) <: Dates.TimeType
        "$(v[1]), …, $(v[end])"
    elseif length(v) > 5
        "$(v[1]), $(v[2]), …, $(v[end-1]), $(v[end])"
    else
        join((string(va) for va in v), ", ")
    end
    printstyled(io, s * svals * "]"; color=:cyan)
end

_printmode(io, mode) = print(io, nameof(typeof(mode)))
_printorder(io, mode) = print(io, nameof(typeof(order(mode))))
_printspan(io, mode) = print(io, nameof(typeof(span(mode))))
_printsampling(io, mode) = print(io, nameof(typeof(sampling(mode))))

function _ioctx(io, T)
    IOContext(io, :compact=>true, :limit=>true, :typeinfo=>T)
end

function Base.show(io::IO, mime::MIME"text/plain", dims::Tuple{<:Dimension,Vararg{<:Dimension}})
    _layout_dims(io, mime, dims)
end


text:
# override Documenter.jl/src/Writers/HTMLWriter.jl

import Documenter: Anchors, Builder, Documents, Expanders, Documenter, Utilities, Writers
import Documenter.Utilities.DOM: DOM, Tag, @tags
import Documenter.Writers.HTMLWriter: pagetitle, domify, open_output
import Documenter.Writers.HTMLWriter: navhref, relhref
import Documenter.Writers.HTMLWriter: mdconvert, mdflatten
import Documenter.Writers.HTMLWriter: getpage, get_url
import Documenter.Writers.HTMLWriter: render_head, render_topbar, render_navmenu, render_article
import Documenter.Writers.HTMLWriter: normalize_css, google_fonts, fontawesome_css, highlightjs_css
import Documenter.Writers.HTMLWriter: analytics_script, requirejs_cdn, asset_links
import Documenter.Writers.HTMLWriter: canonical_link_element

function Documenter.Writers.HTMLWriter.render_head(ctx, navnode)
    @tags head meta link script title
    src = get_url(ctx, navnode)

    page_title = "$(mdflatten(pagetitle(ctx, navnode))) · $(ctx.doc.user.sitename)"
    css_links = [
        normalize_css,
        google_fonts,
        fontawesome_css,
        highlightjs_css,
    ]
    head(
        meta[:charset=>"UTF-8"],
        meta[:name => "viewport", :content => "width=device-width, initial-scale=1.0"],
        title(page_title),

        analytics_script(ctx.settings.analytics),

        canonical_link_element(ctx.settings.canonical, src),

        # Stylesheets.
        map(css_links) do each
            link[:href => each, :rel => "stylesheet", :type => "text/css"]
        end,

        script("documenterBaseURL=\"$(relhref(src, "."))\""),
        script[
            :src => requirejs_cdn,
            Symbol("data-main") => relhref(src, ctx.documenter_js)
        ],

        script[:src => relhref(src, "siteinfo.js")],
        script[:src => relhref(src, "../versions.js")],

        # Custom user-provided assets.
        asset_links(src, ctx.local_assets),

        # juliakorea custom css
        link[:href => relhref(src, "assets/custom.css"), :rel => "stylesheet", :type => "text/css"],

        # juliakorea Korean word break
        script[:src => "/js/jquery-1.8.3.min.js"],
        script[:src => "/js/jquery.word-break-keep-all.min.js"],
        script("\$(document).ready(function() { \$('p').wordBreakKeepAll(); });")
    )
end

function Documenter.Writers.HTMLWriter.render_page(ctx, navnode)
    @tags html body

    page = getpage(ctx, navnode)

    head = render_head(ctx, navnode)
    navmenu = render_navmenu(ctx, navnode)
    article = render_article(ctx, navnode)

    htmldoc = DOM.HTMLDocument(
        html[:lang=>"ko"](
            head,
            body(navmenu, article)
        )
    )

    open_output(ctx, navnode) do io
        print(io, htmldoc)
    end
end

const t_Previous = "이전글"
const t_Next     = "다음글"
t_Edit_on(host)  = " Edit on $host"

function Documenter.Writers.HTMLWriter.render_article(ctx, navnode)
    @tags article header footer nav ul li hr span a

    header_links = map(Documents.navpath(navnode)) do nn
        title = mdconvert(pagetitle(ctx, nn); droplinks=true)
        nn.page === nothing ? li(title) : li(a[:href => navhref(ctx, nn, navnode)](title))
    end

    topnav = nav(ul(header_links))

    # Set the logo and name for the "Edit on.." button.
    host_type = Utilities.repo_host_from_url(ctx.doc.user.repo)
    if host_type == Utilities.RepoGitlab
        host = "GitLab"
        logo = "\uf296"
    elseif host_type == Utilities.RepoGithub
        host = "GitHub"
        logo = "\uf09b"
    elseif host_type == Utilities.RepoBitbucket
        host = "BitBucket"
        logo = "\uf171"
    else
        host = ""
        logo = "\uf15c"
    end
    hoststring = isempty(host) ? " source" : " on $(host)"

    if !ctx.settings.disable_git
        pageurl = get(getpage(ctx, navnode).globals.meta, :EditURL, getpage(ctx, navnode).source)
        if Utilities.isabsurl(pageurl)
            url = pageurl
        else
            if !(pageurl == getpage(ctx, navnode).source)
                # need to set users path relative the page itself
                pageurl = joinpath(first(splitdir(getpage(ctx, navnode).source)), pageurl)
            end
            url = Utilities.url(ctx.doc.user.repo, pageurl, commit=ctx.settings.edit_branch)
        end
        if url !== nothing
            edit_verb = (ctx.settings.edit_branch === nothing) ? "View" : "Edit"
            push!(topnav.nodes, a[".edit-page", :href => url](span[".fa"](logo), " $(edit_verb)$hoststring"))
        end
    end
    art_header = header(topnav, hr(), render_topbar(ctx, navnode))

    # build the footer with nav links
    art_footer = footer(hr())
    if navnode.prev !== nothing
        direction = span[".direction"](t_Previous)
        title = span[".title"](mdconvert(pagetitle(ctx, navnode.prev); droplinks=true))
        link = a[".previous", :href => navhref(ctx, navnode.prev, navnode)](direction, title)
        push!(art_footer.nodes, link)
    end

    if navnode.next !== nothing
        direction = span[".direction"](t_Next)
        title = span[".title"](mdconvert(pagetitle(ctx, navnode.next); droplinks=true))
        link = a[".next", :href => navhref(ctx, navnode.next, navnode)](direction, title)
        push!(art_footer.nodes, link)
    end

    pagenodes = domify(ctx, navnode)
    article["#docs"](art_header, pagenodes, art_footer)
end


text:
module SocpSolver

export Problem, LinearProgram, QuadraticProgram, Settings, solve

import Base: size, copy, iterate, length, getindex, ∈, *, \, /
using Printf
using Parameters
using SparseArrays
using LinearAlgebra
using AMD
using MathOptInterface

include("utils.jl")
include("qdldl.jl")
include("cone.jl")
include("problem.jl")
include("types.jl")
include("scaling.jl")
include("settings.jl")
include("solution.jl")
include("equil.jl")
include("kkt.jl")
include("abs_kkt.jl")
include("term.jl")
include("ipm.jl")
include("printing.jl")
include("moi_wrapper.jl")

end


text:
using Cliffords, Formatting

function pauli_rand_clifford_circuit(pseqs, ofile)
    # Pauli ops
    ps = [Id,X,Y,Z]
    # Pauli ops as Clifford group ops
    c_ps = map(localclifford,[1,3,6,9])

    for n in 1:length(pseqs)
        prop = Id # Pauli correction propagated through sequence
        rp   = Id
        for i in 1:length(pseqs[n])-1
            c = pseqs[n][i] == 0 ? localclifford(1) : localclifford(pseqs[n][i])
            # pick a new random pauli
            rp = rand(ps)
            # combine new Pauli with propagated Pauli to get new propagated Pauli
            # compute new "compiled" pulse
            prop, comp = propagate_and_compile(c,rp,prop)
            @printf ofile "%d," comp
        end
        c = pseqs[n][end] == 0 ? localclifford(1) : localclifford(pseqs[n][end])
        rp = rand(ps)
        prop, comp = propagate_and_compile(c,rp,prop)
        @printf ofile "%d\n" Cliffords.localcliffordindex(Clifford(rand([Id,Z])*prop)*(c*Clifford(rp)))
        #println(Cliffords.localcliffordindex(Clifford(rand([Id,Z])*prop)*(c*Clifford(rp))))
    end
end

function pauli_rand_clifford_circuit2(pseqs, ofile)
    # Pauli ops
    paulis = [Id,X,Y,Z]
    # Pauli ops as Clifford group ops
    cliffs = [1,3,6,9]
    c_ps = map(localclifford,[1,3,6,9])

    for n in 1:length(pseqs)
        prop = Id # Pauli correction propagated through sequence
        rp   = Id
        for i in 1:length(pseqs[n])-1
            c = pseqs[n][i] == 0 ? localclifford(1) : localclifford(pseqs[n][i])
            # pick a new random pauli
            rp = rand(1:4)
            # combine new Pauli with propagated Pauli to get new propagated Pauli
            # compute new "compiled" pulse
            prop, _ = propagate_and_compile(c,paulis[rp],prop)
            @printf ofile "%d,%d," cliffs[rp] pseqs[n][i]
        end
        c = pseqs[n][end] == 0 ? localclifford(1) : localclifford(pseqs[n][end])
        rp = rand(1:4)
        prop, _ = propagate_and_compile(c,paulis[rp],prop)
        @printf ofile "%d,%d,%d\n" cliffs[rp] pseqs[n][end] Cliffords.localcliffordindex(Clifford(rand([Id,Z])*prop))
        # println(cliffs[rp],",",pseqs[n][end],",",Cliffords.localcliffordindex(Clifford(rand([Id,Z])*prop)))
    end
end

# As is, the clifford multiplication is still expensive, so we need to memoize it (cache it)
# in order to speed up the computation.
let
    global propagate_and_compile

    table = Dict{Tuple{Clifford,Pauli,Pauli},Tuple{Pauli,Int}}()

    function propagate_and_compile(cliff,randp,propp)
        if !((cliff,randp,propp) in keys(table))
            prop = cliff*(randp*propp)
            table[(cliff,randp,propp)] = (prop,Cliffords.localcliffordindex(cliff*Clifford(randp)))
        end
        return table[(cliff,randp,propp)]
    end
end

begin
    # set random seed (useful for testing)
    srand(0xcafe)

    infile = open(ARGS[1])
    seqs = readlines(infile)
    close(infile)

    count = parse(Int,ARGS[2])

    pseqs = [map(x->parse(Int,x), split(s,',')) for s in seqs];

    L = ceil(Int,log(10,count)+1)
    for i in 1:count
        label = format("{1:0$(L)d}",i)
        println(replace(ARGS[1],r".csv$","-r$(label).csv"))
        rfile = open(replace(ARGS[1],r".csv$","-r$(label).csv"),"w")
        pauli_rand_clifford_circuit(pseqs, rfile)
        close(rfile)
    end
end

max_stars_repo_name: marcusps/RGST-analysis

max_issues_repo_name: marcusps/RGST-analysis

max_forks_repo_name: marcusps/RGST-analysis


text: println("This char is illegal:", "ç")

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
"""
	parratio(eigvec::AbstractMatrix{<:Number})

Compute participation ratio of an array of eigenvectors, supplied as a Matrix
"""
function parratio(eigvec::AbstractMatrix{<:Number})
	denom = sum(abs2, eigvec, dims=2).^2
	numer = sum(x->abs2(x)^2, eigvec,dims=2)

	return numer./denom
end

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
using MonteCarloMeasurements, Distributions
particles(m)
# (x = 0.00257 ± 1.0, y = Particles{Float64,1000}[1.0 ± 1.7, 1.0 ± 1.7, 1.0 ± 1.7, 0.998 ± 1.7, 1.0 ± 1.8, 1.0 ± 1.7, 0.998 ± 1.7, 0.997 ± 1.7, 0.999 ± 1.7, 0.996 ± 1.8])
# The way I get there... it's not so pretty. Let's start with μ and σ:

μ = Particles(1000,Normal(0,1))
# Part1000(0.000235 ± 1.0)

σ = Particles(1000,Normal(0,1))^2
# Part1000(1.008 ± 1.5)
# Surprisingly, this works just fine:

@unsafe d = Normal(μ,σ)
# Normal{Particles{Float64,1000}}(
# μ: 0.000235 ± 1.0
# σ: 1.01 ± 1.5
# )
# Now trying the obvious thing,

@unsafe x = rand(d)
# Part1000(-0.03248 ± 1.0)
# It thinks it works, but it's really horribly broken:

(x-μ)/σ
# Part1000(-0.03246 ± 2.6e-8)
# So I got it to work with a little helper function:

N = 1000
parts(x::Normal{P} where {P <: AbstractParticles}) = Particles(length(x.μ), Normal()) * x.σ + x.μ
# parts(x::Sampleable{F,S}) where {F,S} = Particles(N,x)
# parts(x::Integer) = parts(float(x))
# parts(x::Real) = parts(repeat([x],N))
# parts(x::AbstractArray) = Particles(x)
# parts(p::Particles) = p

@unsafe x = parts(d)
# Part1000(0.05171 ± 2.37)

(x-μ)/σ
# Part1000(0.0007729 ± 1.0)
# Much better.

# This is fine for one distribution, but most don't compose as nicely as a normal.

# Many distributions break entirely, seemingly because of argument checking in Distributions.jl:

@unsafe Bernoulli(Particles(1000,Beta(2,3)))
# ERROR: TypeError: non-boolean (Particles{Float64,1000}) used in boolean context

import MonteCarloMeasurements.particletype
struct ParticleDistribution{D}
    d::D
end

function ParticleDistribution(d::Type{<:Distribution}, p...)
    @unsafe ParticleDistribution(d(p...))
end
particletype(pd::ParticleDistribution) = particletype(getfield(pd.d,1))

function Base.rand(d::ParticleDistribution{D}) where D
    T,N = particletype(d)
    i = rand(1:N)
    d2 = MonteCarloMeasurements.replace_particles(d.d,P->P isa AbstractParticles, P->P[i])
    rand(d2)
end

pd = ParticleDistribution(Bernoulli, Particles(1000,Beta(2,3)))
@btime rand($pd)
@btime rand(Bernoulli(0.3))
@code_warntype rand(pd)
@code_warntype MonteCarloMeasurements.replace_particles(pd.d,P->P isa AbstractParticles, P->P[1])

using Distributions, MonteCarloMeasurements, LinearAlgebra, TransformVariables

z = Particles(1000,MvNormal(zeros(2,2) + I))
t = as((y=asℝ₊,p=as𝕀))

t([0.0,0.0])
# (y = 1.0, p = 0.5)

t([0.0,0.0]) |> inverse(t)

foreach(register_primitive, [<=, >=, <, >])

tz = t(z)
it = inverse(t)
@bymap it(tz)

using MonteCarloMeasurements, Test
p = 0 ± 1

@test MonteCarloMeasurements.Ngetter([p,p],p) == 500
@test MonteCarloMeasurements.Ngetter([p,p]) == 500

f(x) = 2x
f(x,y) = 2x + y

@test f(p) ≈ @bymap f(p)
@test f(p,p) ≈ @bymap f(p,p)

@bymap f(p,p)

##
using Zygote, MonteCarloMeasurements
p = 1 ± 0.1
f(x) = x^2
f'(p)


text:
"""
    optimize!(fdf, M::Wrapper, x0)

Find an optimizer for `fdf`, starting with the initial approximation `x0`.
`fdf(x, g)` must return a tuple (f(x), ∇f(x)) and, if `g` is mutable, overwrite
it with the gradient.
"""
function optimize!(fdf, M::Wrapper, x0; reset=true, tracking=nothing, verbosity=0)
    return __optim_impl!(fdf, M, x0, tracking; reset, verbosity)
end

function __optim_impl!(fdf, M::Wrapper, x0, logger::AbstractLogger; reset=true)
    optfn! = OptFunc(fdf, M)

    with_logger(logger) do
        n = length(x0)
        @logmsg OptLogLevel """

        ==OPTIMIZATION START==
        # First $n value$(n == 1 ? "" : "s") - argument vector
        # next value - function value
        # last $n value$(n == 1 ? "" : "s") - gradient vector
        """
        @logmsg OptLogLevel "==SOLVER INITIALIZATION=="
        init!(optfn!, M, x0; reset)
        @logmsg OptLogLevel "==SOLVER INITIALIZED=="
        for niter in Iterators.countfrom(1)
            @logmsg OptLogLevel "# Iteration $niter"
            step!(optfn!, M)
            stopcond(M) && break
        end
        return convstat(M)
    end
end

function __optim_impl!(
    fdf, M::Wrapper, x0, tracking::IO;
    reset=true, verbosity, kw...
)
    logger = ConsoleLogger(
        tracking,
        LogLevel(-10 * verbosity);
        meta_formatter=metafmt_noprefix_nosuffix, show_limited=false,
    )
    __optim_impl!(fdf, M, x0, logger; reset)
end

function __optim_impl!(
    fdf, M::Wrapper, x0, tracking::AbstractString;
    reset=true, verbosity::Integer
)
    open(tracking, "w") do log_io
        logger = ConsoleLogger(
            log_io,
            LogLevel(-10 * verbosity);
            meta_formatter=metafmt_noprefix_nosuffix, show_limited=false,
        )
        return __optim_impl!(fdf, M, x0, logger; reset)
    end
end

function __optim_impl!(fdf, M::Wrapper, x0, ::Nothing; reset=true, kw...)
    __optim_impl!(fdf, M, x0, NullLogger(); reset)
end

"""
    optimize!(
        fdf, M::OptBuffer, x₀;
        gtol=1e-6,
        convcond=nothing,
        maxiter=100,
        maxcalls=nothing,
        reset=true,
        constrain_step=nothing,
        tracking=nothing,
        verbosity=0
    )

Find an optimizer for `fdf`, starting with the initial approximation `x₀`.

# Arguments
- `M::OptBuffer`: the core method to use for optimization
- `fdf(x, g)::Function`: function to optimize. It must return a tuple (f(x), ∇f(x)) and,
    if `g` is mutable, overwrite
    it with the gradient.
- `x0`: initial approximation

# Keywords

## Convergence criteria

There are two options to specify convergence criterion.
The default is by `gtol` and the second by custom stop `convcond`.

- `gtol::Real`: (default stop criterion) stop optimization when the gradient's 2-norm is less
- `convcond=(x, xpre, y, ypre, g)->Bool`: function, custom stop criterion based on argument
    values, function values and `g`radient. If `nothing` (default), corresponds to `gtol`,
    and when specified, the `gtol`-criterion is ignored.

Example (default criterion): `convcond = (x, xpre, y, ypre, g) -> norm(g, 2) ≤ gtol`.

## Limitting optimization

(Un)Limit optimization process by specifing either `maxiter` and/or `maxcalls`.

- `maxiter::Integer`: force stop optimization after this number of iterations
    (use `nothing` or a negative value to not constrain iteration number)
- `maxcalls::Integer`: force stop optimization after this number of function calls
    (use `nothing` or a negative value to not constrain call number)

## Optimization constrains

The inequality constrains of optimization is handled by `constrain_step`.

- `constrain_step(x0, d)`: a function to constrain step from `x0` in the direction `d`.
    It must return a real-numbered value `α` such that `x0 + αd` is the maximum allowed step

## Initializing

- `reset=true`: a value to pass as a keyword argument to the optimizer `init!` method

## Optimization path

- `tracking::Union{Nothing,IO,AbstractString}`: IO stream or a file name to log the
    optimization process or `nothing` to disable logging (default: `nothing`)
- `verbosity::Integer`: verbosity of logging. `0` (default) disables tracking. `1` logs all
    points of objective function evaluation with corresponding values and gradients.
    `2` shows additional statistics regarding the line search. Option ignored if
    `tracking == nothing`.
"""
function optimize!(
    fdf, M::OptBuffer, x0;
    gtol=convert(float(eltype(x0)), 1e-6),
    convcond=nothing,
    maxiter::Integer=100,
    maxcalls=nothing,
    reset=true,
    constrain_step=nothing,
    tracking=nothing,
    verbosity::Integer=0,
)
    grad_tol = (isnothing(gtol) || gtol < 0) ? zero(eltype(x0)) : gtol
    convcond = isnothing(convcond) ? stopbygradient(grad_tol) : convcond
    iter_limit = (isnothing(maxiter) || maxiter < 0) ? typemax(Int) : convert(Int, maxiter)
    call_limit = (isnothing(maxcalls) || maxcalls < 0) ?
        typemax(Int) :
        convert(Int, maxcalls)
    M = BasicConvergenceStats(M; convcond, iter_limit, call_limit)
    if !isnothing(constrain_step)
        M = ConstrainStepSize(constrain_step, M)
    end
    optimize!(fdf, M, x0; reset, tracking, verbosity)
end

function optimize!(fdf, M::Type{<:OptBuffer}, x0; kw...)
    return optimize!(fdf, M(x0), x0; kw...)
end

"""
    optimize(
        fdf, x₀;
        method,
        kw...
    )

Find an optimizer for `fdf`, starting with the initial approximation `x₀`.
    `method` keyword chooses a specific optimization method. See [`optimize!`](@ref) for
    the description of other keywords.
"""
function optimize(fdf, x0; method, kw...)
    return optimize!(fdf, method, x0; kw...)
end

"""
    Downhill.solver(
        M::OptBuffer;
        gtol = 1e-6,
        maxiter = 100,
        maxcalls = nothing,
        constrain_step=nothing,
    )

Return the wrapper object for a chosen method to solve an optimization problem with given
    parameters. For the description of keywords, see [`optimize!`](@ref)
"""
function solver(
    M::OptBuffer;
    gtol=convert(eltype(argumentvec(M)), 1e-6),
    convcond=nothing,
    maxiter=100,
    maxcalls=nothing,
    constrain_step=nothing,
)
    grad_tol = (isnothing(gtol) || gtol < 0) ? zero(eltype(x0)) : gtol
    convcond = isnothing(convcond) ? stopbygradient(grad_tol) : convcond
    iter_limit = (isnothing(maxiter) || maxiter < 0) ? typemax(Int) : convert(Int, maxiter)
    call_limit = (isnothing(maxcalls) || maxcalls < 0) ?
        typemax(Int) :
        convert(Int, maxcalls)
    M = BasicConvergenceStats(M; convcond, iter_limit, call_limit)
    if !isnothing(constrain_step)
        M = ConstrainStepSize(constrain_step, M)
    end
    return M
end

"""
    Downhill.solver(
        M::DataType, x;
        gtol=1e-6, maxiter = 100, maxcalls = nothing, constrain_step)

Return the wrapper object for a chosen method to solve an optimization problem with given
    parameters compatible with the dimensions of `x`.
"""
function solver(M::Type{<:OptBuffer}, x::AbstractVector; kw...)
    return solver(M(x); kw...)
end


text:
using PyPlot
PyPlot.ioff()
push!(LOAD_PATH, "/home/berceanu/Development/topo-photon/anc/scripts")
import HH
import BP
using Base.Test

# calculates zero-point energy, with and without nonabelian correction
ηzpenew(q::Int, κ::Float64) = (α=1/q; 4π*α/κ * endiffnonab(q, κ))
ηzpeold(q::Int, κ::Float64) = (α=1/q; 4π*α/κ * endiff(q, κ))
ηzpeterm(q::Int, κ::Float64) = (α=1/q; 4π*α/κ * endiffterm(q, κ))
# calculates energy difference between numerical (exact) and
#    theoretical energies, with the nonab corr.
endiffnonab(q::Int, κ::Float64) = er(q,κ) - (et(q,κ) + δE(q,κ))
# calculates energy difference between numerical (exact) and
#    theoretical energies (without nonab corr)
endiff(q::Int, κ::Float64) =  er(q,κ) - et(q,κ)
# calculates energy difference between numerical (exact) and
#    theoretical energies, with first term of the nonabelian correction
endiffterm(q::Int, κ::Float64) = er(q,κ) - (et(q,κ) + δEterm(q, κ))
et(q::Int,κ::Float64) =  et(q,1,κ)
# calculates theoretical energy, without non-abelian correction
function et(q::Int,p::Int,κ::Float64)
    α = p/q
    gs = Array(Float64, 25, 25)
    HH.hhgrstate!(gs, p, q)
    e1 = mean(gs)
    e1 + 1/2*κ/(2π*α)
end
er(q::Int,κ::Float64) = er(q,1,15,κ)
# calculates numerical (exact) energy
function er(q::Int,p::Int,N::Int,κ::Float64)
    M = spzeros(Complex{Float64}, N^2,N^2)
    α = p/q
    BP.buildham_exact!(M, N,α,κ)
    real(eigs(M, nev=1, which=:SR, ritzvec=false)[1][1])
end
# calculates average over MBZ of the non-abelian correction to the
#    1(st) band for p=1 for certain trap
δE(q::Int,κ::Float64) = δE(1,q,1, collect(linspace(-π/q, π/q, 20)),
   collect(linspace(-π, π, 20)), κ)
# calculates the average (over specified points) non-abelian energy
#    correction to n(th) band
δE(n::Int,q::Int,p::Int, px::Array{Float64, 1},py::Array{Float64,
   1},κ::Float64) = mean([δE(n,q,p, x,y,κ) for y in py, x in px])
# calculates non-abelian energy correction to 1(st) band considering
#    only the first term of the sum
δEterm(q::Int, κ::Float64) = mean([δEterm(q, x,y, κ) for y in
   linspace(-π, π, 20), x in linspace(-π/q, π/q, 20)])
# calculates non-abelian energy correction to n(th) band at position
#    (kₓ⁰,ky) in the MBZ
δE(n::Int,q::Int,p::Int, k₀x::Float64,ky::Float64,κ::Float64) =
   κ/2*sum([n′ != n ? norm(A(n,n′,q,p, k₀x,ky))^2 : 0.0 for n′ in 1:q])
# calculates non-abelian energy correction to 1(st) band at position
#    (kₓ⁰,ky) in the MBZ, considering only the first term of the sum
δEterm(q::Int, k₀x::Float64,ky::Float64,κ::Float64) = κ/2 *
   norm(A(1,2,q,1, k₀x,ky))^2
# calculates Berry connection A (vector quantity) at position (kₓ⁰,ky)
#    in the MBZ for bands n and n′
function A(n::Int,n′::Int,q::Int,p::Int, k₀x::Float64,ky::Float64)
    # initializing hamiltonian matrices
    for M = (:H, :∇Hx, :∇Hy)
        @eval ($M) = zeros(Complex{Float64}, ($q,$q))
    end
    # top right corner
    H[1,q] = -exp(-im*ky)
    ∇Hy[1,q] = im*exp(-im*ky) # ∂ky
    #bottom left corner
    H[q,1] = -exp(im*ky)
    ∇Hy[q,1] = -im*exp(im*ky) # ∂ky
    #upper subdiagonal
    ius = diagind(H,1)
    #lower subdiagonal
    ils = diagind(H,-1)
    H[ius] = -exp(im*ky)
    H[ils] = -exp(-im*ky)
    ∇Hy[ius] = -im*exp(im*ky)
    ∇Hy[ils] = im*exp(-im*ky)
    #main diagonal
    α = p/q
    for j in 1:q
        H[j,j] =  -2*cos(k₀x + 2*π*α*j)
      ∇Hx[j,j] =   2*sin(k₀x + 2*π*α*j)
    end
    # diagonalize HH Hamiltonian -> E's and u's (eigs and eigvs)
    F = eigfact(H)
    U = F[:vectors]
    E = F[:values]
    # calculate denominator Float64
    denominator = E[n′] - E[n]
    # calculate expectation value along kₓ (xnumerator)
    xnumerator = dot(U[:,n], ∇Hx*U[:,n′])
    # calculate expectation value along ky (ynumerator)
    ynumerator = dot(U[:,n], ∇Hy*U[:,n′])
    return [im * xnumerator / denominator, im * ynumerator / denominator]
end

#plotting
qs = 4:20
y1 = [ηzpeold(q, 0.02)::Float64 for q in qs]
y2 = [ηzpenew(q, 0.02)::Float64 for q in qs]
# plot also energy correction using just first term of the sum
y4 = [ηzpeterm(q, 0.02)::Float64 for q in qs]
y3 = HH.ηzpe(qs,0.02)
@test_approx_eq_eps y1 y3 1e-6

# matrix holding level spacing errors as columns
ηL = Array(Float64, length(qs),4)
# populating matrix
for i = 1:4
    ηL[:,i] = HH.ηlev(qs,0.02, i-1,i)
end

# matplotlib parameters
matplotlib["rcParams"][:update](Dict("axes.labelsize" => 22,
                                     "axes.titlesize" => 20,
                                     "font.size" => 18,
                                     "legend.fontsize" => 14,
                                     "axes.linewidth" => 1.5,
                                     "font.family" => "serif",
                                     "font.serif" => "Computer Modern Roman",
                                     "xtick.labelsize" => 20,
                                     "xtick.major.size" => 5.5,
                                     "xtick.major.width" => 1.5,
                                     "ytick.labelsize" => 20,
                                     "ytick.major.size" => 5.5,
                                     "ytick.major.width" => 1.5,
                                     "text.usetex" => true,
                                     "figure.autolayout" => true))

## checking δE(kₓ,k_y) flat for q = 5
# system parameters
q = 5
r = 11 # points in MBZ

# N should be an odd multiple of q
N = r*q # zero-padded system size
l = div(N-1,2)
x = -l:l
δk = 2π/N #resolution in mom space
k = x * δk

#generate all p values inside MBZ
xmbz = -div(r-1,2):div(r-1,2)
kxmbz = xmbz * δk

data = [δE(1,5,1, x,y,0.02)::Float64 for y in k, x in kxmbz]
a, b = extrema(data)

fig, ax = plt[:subplots](figsize=(5, 5))

img = ax[:imshow](data, origin="upper", ColorMap("gist_heat_r"),
                 interpolation="none",
                 extent=[-π/q, π/q, -π, π],
                 aspect=1/q)
ax[:set_ylabel](L"$p_y$",labelpad=-1)
ax[:set_xlabel](L"$p_x^0$", labelpad=-1)
ax[:xaxis][:set_ticks]([-π/q,0,π/q])
ax[:xaxis][:set_ticklabels]([L"$-\pi/5$", L"$0$", L"$\pi/5$"])
ax[:yaxis][:set_ticks]([-π,0,π])
ax[:yaxis][:set_ticklabels]([L"$-\pi$", L"$0$", L"$\pi$"])
cbaxes = fig[:add_axes]([0.25, 0.04, 0.65, 0.015])
cbar = fig[:colorbar](img, cax=cbaxes, orientation="horizontal")
cbar[:set_ticks]([a,b])
cbar[:set_ticklabels]([L"$5.378 \times 10^{-3}$", L"$11.680 \times 10^{-3}$"])
cbar[:set_label](L"$\delta E(5,0.02)$", rotation=0, labelpad=-20, y=0.5)
cbar[:solids][:set_edgecolor]("face")
fig[:savefig]("../../figures/correction_mbz.pdf", transparent=true,
              pad_inches=0.0, bbox_inches="tight")
plt[:close](fig)

# various line types
lines = ["-","--","-.",":"]

fig, axes = plt[:subplots](2, figsize=(8, 5))
for (i, ax) in enumerate(axes)
    if i == 1 # first panel, with zero-point-energy error
        # marker = "o"
        ax[:plot](qs, y1, "black", ls="dashed", label = "old") # $E_{ex} - E_{th}$
        ax[:plot](qs, y2, "black", label="new") # $E_{ex} - (E_{th} + \delta E)$
        ax[:set_xticklabels]([])
        ax[:yaxis][:set_ticks]([-3.6, -2.4, -1.2, 0., 1.2])
        ax[:set_yticklabels]([L"$-3.6$", L"$-2.4$", L"$-1.2$", L"$0$",
           L"$1.2$"])
        ax[:set_ylabel](L"$\eta_{zpe}$")
    else # second panel, with level spacing error
        for i = 1:4
            ax[:plot](qs, ηL[:,i], "black", ls=lines[i]) # $\kappa=0.02$
        end
        ax[:set_ylabel](L"$\eta_{lev}$")
        ax[:yaxis][:set_ticks]([-1, 0, 1, 2, 3])
        ax[:set_xlabel](L"$q$")
    end
    ax[:set_xlim](qs[1], qs[end])
end
fig[:savefig]("../../figures/nonabcorr.pdf", transparent=true,
   pad_inches=0.0, bbox_inches="tight")
plt[:close](fig)

max_stars_repo_name: berceanu/topo-photon

max_issues_repo_name: berceanu/topo-photon

max_forks_repo_name: berceanu/topo-photon


text:
#Define function to convert independent variational parameters
#to network parameters a, b and W with translation symmetry
function WToPar(RBMPar::RBM_par, WRBM::Array{Complex{Float64},1})

  C = Array{Complex{Float64}}(undef,0)
  W_pre = Array{Array{Int64}}(undef,alpha)

  for h=1:alpha
    W_pre[h] = WvDerivativeTi.cod(h-1,L_x,L_y)
  end

  for i=1:alpha
    append!(C,WRBM[i].*ones(nspins))
  end

  W_s = WRBM[alpha+1:end]
  RBMPar.b[:] = C

  for i=1:nspins
    for h=1:alpha
       RBMPar.W[i,1+nspins*(h-1):h*nspins]=W_s[W_pre[h][:,i]]
    end
  end

end

#Define function to convert network parameters a,b and W
#to independent variational parameters with translation symmetry
function ParToW(RBMPar::RBM_par,W_RBM::Array{Complex{Float64},1})

    for h=1:alpha
        W_RBM[h] = RBMPar.b[1 + (h-1)*nspins]
    end

    W_RBM[alpha+1:end] = RBMPar.W[1,:]

end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
import Base.push!

function push!(B :: Union{UniformScaling{T},Matrix{T}},
               sk:: Vector{T},
               yk:: Vector{T}) where T
#               scaling :: Bool = true) where T

    scaling = true
    denom = yk'*sk
    if (denom > 1.0e-20)
        #self-scaled version aka Luenberger

        #divide vector yk by denom
        By = B*(yk/denom)

        ytBy = (yk'*By) # inspired by Luenberger
        #ytBy = (yk'*yk)  # similar to L-BFGS

        γ = 1.0
        if scaling
            γ = T(1) / (ytBy)
        end

        # refactor to yield O(n²) complexity... Much more efficient!
        Byst = By*sk'

        #divide vector sk by denom
        sst = (sk/denom)*sk'

        #B = γ * (B - (Byst' + Byst)/denom + (ytBy/denom)*sst/denom)  + sst/denom
        # avoids dividing matrices by denom
        #
        B = γ * (B  - (Byst' + Byst) + ytBy*sst)  + sst
        B = T(0.5)*(B + B')  # make sure B is symmetric

        # Equivalent math formulation, for testing purposes.
        # math formula, involves matrix multiplications, thus O(n³) complexity
        #M = I - sk*yk'/denom
        #B1 = γ * M*B1*M' + sk*sk'/denom
        #B1 = T(0.5)*(B1 + B1')  # make sure B is symmetric
        #println( " In push! différence entre B et B1 : " , norm(B - B1))
    else
        @warn "No update,  denom = ", denom
        #B = Matrix(eltype(yk).(I(length(yk))))
    end

    return B
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
using HttpCommon
using Compat
using Test

# headers
h = HttpCommon.headers()
@test isa(h, Headers)
@test isa(HttpCommon.headers(), Headers)
# Request
@test sprint(show, Request()) == "Request(:/, 0 headers, 0 bytes in body)"
@test sprint(show, Request("GET", "/get", HttpCommon.headers(), "")) ==
        "Request(:/, 4 headers, 0 bytes in body)"

# Cookie
@test sprint(show, Cookie("test","it")) ==
        "Cookie(test, it, 0 attributes)"

# Response
h_abstract = Dict{AbstractString, AbstractString}(h)
@test sprint(show, Response(200, h, "HttpCommon")) ==
        "Response(200 OK, 4 headers, 10 bytes in body)"
@test sprint(show, Response(200, h_abstract, "HttpCommon")) ==
        "Response(200 OK, 4 headers, 10 bytes in body)"
@test sprint(show, Response(200, h, UInt8[1,2,3])) ==
        "Response(200 OK, 4 headers, 3 bytes in body)"
@test sprint(show, Response(200, h_abstract, UInt8[1,2,3])) ==
        "Response(200 OK, 4 headers, 3 bytes in body)"
@test sprint(show, Response(200, h)) ==
        "Response(200 OK, 4 headers, 0 bytes in body)"
@test sprint(show, Response(200, h_abstract)) ==
        "Response(200 OK, 4 headers, 0 bytes in body)"
@test sprint(show, Response(200, "HttpCommon")) ==
        "Response(200 OK, 4 headers, 10 bytes in body)"
@test sprint(show, Response(200, UInt8[1,2,3])) ==
        "Response(200 OK, 4 headers, 3 bytes in body)"
@test sprint(show, Response("HttpCommon", h)) ==
        "Response(200 OK, 4 headers, 10 bytes in body)"
@test sprint(show, Response("HttpCommon", h_abstract)) ==
        "Response(200 OK, 4 headers, 10 bytes in body)"
@test sprint(show, Response(UInt8[1,2,3], h)) ==
        "Response(200 OK, 4 headers, 3 bytes in body)"
@test sprint(show, Response(UInt8[1,2,3], h_abstract)) ==
        "Response(200 OK, 4 headers, 3 bytes in body)"
@test sprint(show, Response("HttpCommon")) ==
        "Response(200 OK, 4 headers, 10 bytes in body)"
@test sprint(show, Response(UInt8[1,2,3])) ==
        "Response(200 OK, 4 headers, 3 bytes in body)"
@test sprint(show, Response(200)) ==
        "Response(200 OK, 4 headers, 0 bytes in body)"
@test sprint(show, Response()) ==
        "Response(200 OK, 4 headers, 0 bytes in body)"
@test sprint(show, Response(999)) ==
        "Response(999 Unknown Code, 4 headers, 0 bytes in body)"

# Escape HTML
@test escapeHTML(" foo bar") ==
        "<script type='text/javascript'>alert('sucker');</script> foo bar"

# Parse URL query strings
@test parsequerystring("foo=%3Ca%20href%3D%27foo%27%3Ebar%3C%2Fa%3Erun%26%2B%2B&bar=123") ==
        Dict("foo" => "barrun&++", "bar" => "123")

begin
  substrings = split("a%20=1&b=%202,b,c", ",")
  @test parsequerystring(substrings[1]) == Dict("a " => "1", "b" => " 2")
end
@test parsequerystring("") == Dict()
@test_throws ArgumentError parsequerystring("looknopairs")
@test_throws ArgumentError parsequerystring("good=pair&badpair")


text:
export addblock!!,
       setindex,
       setindex!!,

# blocksparse/blocksparsetensor.jl
       randomBlockSparseTensor,

# dense.jl
       randomTensor,
       randomDenseTensor,

# empty.jl
       Empty,
       EmptyTensor,
       EmptyBlockSparseTensor


text:
@enum GitMode mode_dir=0o040000 mode_normal=0o100644 mode_executable=0o100755 mode_symlink=0o120000 mode_submodule=0o160000
Base.string(mode::GitMode) = string(UInt32(mode); base=8)
Base.print(io::IO, mode::GitMode) = print(io, string(mode))

struct GitTreeEntry
    mode::GitMode
    name::String
    hash::SHA1
end

struct GitTree <: GitObject
    entries::Vector{GitTreeEntry}
end

Base.sizeof(entry::GitTreeEntry) =
    ndigits(UInt32(entry.mode); base=8) + 1 + sizeof(entry.name) + 1 + 20
Base.sizeof(tree::GitTree) = sum(sizeof, tree.entries)

function Base.read(io::IO, ::Type{GitTreeEntry})
    mode = GitMode(parse(UInt32,readuntil(io, ' '),base=8))
    name = readuntil(io,'\0')
    hash = read(io, SHA1)
    return GitTreeEntry(mode, name, hash)
end

function GitTree(data::Vector{UInt8})
    entries = GitTreeEntry[]
    io = IOBuffer(data)
    while bytesavailable(io) > 0
        entry = read(io, GitTreeEntry)
        push!(entries, entry)
    end
    GitTree(entries)
end

function oid(tree::GitTree)
    ctx = SHA.SHA1_CTX()
    SHA.update!(ctx, Vector{UInt8}("tree $(sizeof(tree))\0"))
    for entry in tree.entries
        SHA.update!(ctx, Vector{UInt8}("$(entry.mode) $(entry.name)\0"))
        SHA.update!(ctx, entry.hash.bytes)
    end
    return SHA1(SHA.digest!(ctx))
end

entryorder(entry::GitTreeEntry) = entry.mode == mode_dir ? entry.name*"/" : entry.name

"""
    GitTreeEntry(fpath[, name])

Construct a `GitTreeEntry` object based on `path`. If no `name` is provided, it is
assumed to be the `basename(path)`.
"""
function GitTreeEntry(path::AbstractString, name::AbstractString=basename(filepath))
    mode = isdir(path) ? mode_dir : mode_normal
    hash = oid(path)
    GitTreeEntry(mode, name, hash)
end

function GitTree(path::AbstractString, ignore::Glob.FilenameMatch=fn"")
    names = filter(name -> !occursin(ignore,name), readdir(path))
    entries = map(names) do name
        GitTreeEntry(joinpath(path, name), name)
    end
    sort!(entries, by=entryorder)
    GitTree(entries)
end

function oid(path::AbstractString, ignore::Glob.FilenameMatch=fn".git")
    if isdir(path)
        oid(GitTree(path, ignore))
    else
        oid(GitBlob(read(path)))
    end
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text: iseven(i), isodd(i)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
macro chk1(expr)
    quote
        x,i = $(esc(expr))
        if isnull(x)
            @goto err
        else
            get(x),i
        end
    end
end

function parse(str)
    nd, i = tryparsenext_datetime(str,1)
    if isnull(nd)
        throw(ParseError("Invalid DateTime"))
    end
    get(nd)
end

function parse{S<:AbstractString}(ads::Vector{S})
    l = length(ads)
    a = Array(DateTime, l)
    for i in 1:l
        a[i] = parse(ads[i])
    end
    a
end

function tryparsenext_datetime(str,i)
    R = Nullable{DateTime}
    dy, i = @chk1 tryparsenext_base10(str,i,4)
    c,  i = @chk1 tryparsenext_char(str,i,'-')
    dm, i = @chk1 tryparsenext_base10(str,i,2)
    c,  i = @chk1 tryparsenext_char(str,i,'-')
    dd, i = @chk1 tryparsenext_base10(str,i,2)
    c,  i = @chk1 tryparsenext_char(str,i,('T',' '))
    th, i = @chk1 tryparsenext_base10(str,i,2)
    c,  i = @chk1 tryparsenext_char(str,i,':')
    tm, i = @chk1 tryparsenext_base10(str,i,2)
    c,  i = @chk1 tryparsenext_char(str,i,':')
    ts, i = @chk1 tryparsenext_base10(str,i,2)

    nc, i = tryparsenext_char(str,i,'.')
    if isnull(nc)
        d = DateTime(dy,dm,dd,th,tm,ts)
    else
        tms,i = @chk1 tryparsenext_base10_frac(str,i,3)
        d = DateTime(dy,dm,dd,th,tm,ts,tms)
    end
    return R(d), i

    @label err
    return R(), i
end

@inline function tryparsenext_base10_digit(str,i)
    R = Nullable{Int}
    done(str,i) && @goto err
    c,ii = next(str,i)
    '0' <= c <= '9' || @goto err
    return R(c-'0'), ii

    @label err
    return R(), i
end

@inline function tryparsenext_base10(str,i,n)
    R = Nullable{Int}
    r = 0
    for j = 1:n
        d,i = @chk1 tryparsenext_base10_digit(str,i)
        r = r*10 + d
    end
    return R(r), i

    @label err
    return R(), i
end

@inline function tryparsenext_base10_frac(str,i,maxdig)
    R = Nullable{Int}
    r,i = @chk1 tryparsenext_base10_digit(str,i)
    for j = 2:maxdig
        nd, i = tryparsenext_base10_digit(str,i)
        if isnull(nd)
            for k = j:maxdig
                r *= 10
            end
            return R(r), i
        end
        d = get(nd)
        r = 10*r + d
    end
    return R(r),i

    @label err
    return R(), i
end

@inline function tryparsenext_char(str,i,cc::Char)
    R = Nullable{Char}
    done(str,i) && @goto err
    c,ii = next(str,i)
    c == cc || @goto err
    return R(c), ii

    @label err
    return R(), i
end
@inline function tryparsenext_char(str,i,CC::Tuple{Char,Char})
    R = Nullable{Char}
    done(str,i) && @goto err
    c,ii = next(str,i)
    c == CC[1] || c == CC[2] || @goto err
    return R(c), ii

    @label err
    return R(), i
end


text:
push!(LOAD_PATH,"../src/")

using Documenter
using Sylvia

makedocs(
    sitename = "Sylvia",
    format = :html,
    modules = [Sylvia],
    pages = [
        "Home" => "index.md"
        # "Basics" => "basics.md"
        "Symbols" => "symbols.md"
        "Expressions" => "expressions.md"
        "Contexts" => "contexts.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/dalum/Sylvia.jl.git"
)