text:
# This file contains functions for setting optimization constraints on the FFD
# mapping.

@doc """
### numLinearPlaneConstraints

Returns the number of plane constraints (see setLinearPlaneConstraints! for
further details).

**Arguments**

* `map` : PumiMapping object
* `di` : the index the defines the constrained plane (e.g. `di`=2 -> j planes)

"""->
function numLinearPlaneConstraints(map::PumiMapping{Tffd}, di::Int) where Tffd
  it1 = mod(di,3)+1
  it2 = mod(di+1,3)+1
  num_per_plane = (map.nctl[it1]*map.nctl[it2] - 3)*3
  return map.nctl[di]*num_per_plane
end

@doc """
### countVarsLinearPlaneConstraints!

Updates array `cntEq` with the number of variables involved in each linear-plane
constraint.

**Arguments**

* `map` : PumiMapping object
* `di` : the index the defines the constrained plane (e.g. `di`=2 -> j planes)
* `cntEq` : array with number of variables involved in each equality constraint
* `fncidx` : index to the next equality constraint to add to cntEq

**Returns** `fncidx` + number of linear-plane constraints + 1

"""->
function countVarsLinearPlaneConstraints!(map::PumiMapping{Tffd}, di::Int,
                                          cntEq::AbstractArray{Int,1},
                                          fncidx::Int) where Tffd
  num_cnstr = numLinearPlaneConstraints(map, di)
  cntEq[fncidx:fncidx+num_cnstr-1] = 4
  return fncidx+num_cnstr
end

@doc """
### setLinearPlaneConstraints!

Set constraints that relate CPs in given planes to CPs at the corners of the
plane.  Basically, a slice of CPs in the FFD are constrained to remain in a
plane, and shear and translate according to the CPs at the corners of the plane.

**Arguments**

* `map`  : PumiMapping object
* `di`   : the index the defines the constrained plane (e.g. `di`=2 -> j planes)
* `iLfun`, `jLvar` : the k^th element in the gradient matrix has row index
*                    `iLfun(k)` and column index `jLvar(k)`
* `LinG` : `LinG(k) = A(iLfun(k),jLvar(k))`
* `Flow`, `Fupp` : the upper and lower bounds on all constraints
* `fncidx` : function counter; next function to use in iLfun
* `ptr` : index pointer to next free space in iLfun and jLvar

**Returns** `fncidx` + number of linear-plane constraints + 1, and `ptr` + total
  number of non-zeros in the sparse contraints + 1

"""->
function setLinearPlaneConstraints!(map::PumiMapping{Tffd}, di::Int, iLfun::AbstractArray{Int,1},
    jLvar::AbstractArray{Int,1}, LinG::AbstractArray{Tffd,1},
    Flow::AbstractArray{Tffd,1}, Fupp::AbstractArray{Tffd,1},
    fncidx::Int, ptr::Int) where Tffd

  it1 = mod(di,3)+1
  it2 = mod(di+1,3)+1

  jkm = zeros(Int,3)
  cp1_jkm = zeros(jkm); cp2_jkm = zeros(jkm); cp3_jkm = zeros(jkm)
  cp1_idx = zeros(jkm); cp2_idx = zeros(jkm); cp3_idx = zeros(jkm)
  cp1_xyz = zeros(3); cp2_xyz = zeros(3); cp3_xyz = zeros(3)
  u1 = zeros(3); u2 = zeros(3); u3 = zeros(3); e1 = zeros(3); e2 = zeros(3)
  x = zeros(3)

  # loop over the CP planes in the direction di
  for jdi = 1:map.nctl[di]
    jkm[di] = jdi

    # set the (j,k,m) indices for the CPs at the corners of the plane
    cp1_jkm[di] = jdi
    cp1_jkm[it1] = 1
    cp1_jkm[it2] = 1

    cp2_jkm[di] = jdi
    cp2_jkm[it1] = 1
    cp2_jkm[it2] = map.nctl[it2]

    cp3_jkm[di] = jdi
    cp3_jkm[it1] = map.nctl[it1]
    cp3_jkm[it2] = 1

    # get the coordinates of the CPs at the corners of the plane, and some unit
    # vectors
    cp1_xyz = map.cp_xyz[:,cp1_jkm[1],cp1_jkm[2],cp1_jkm[3]]
    cp2_xyz = map.cp_xyz[:,cp2_jkm[1],cp2_jkm[2],cp2_jkm[3]]
    cp3_xyz = map.cp_xyz[:,cp3_jkm[1],cp3_jkm[2],cp3_jkm[3]]

    u1[:] = cp2_xyz - cp1_xyz
    fac1 = 1./norm(u1)
    u1[:] *= fac1

    u2[:] = cp3_xyz - cp1_xyz
    fac2 = 1./norm(u2)
    u2[:] *= fac2

    u3[:] = cross(u1,u2)
    u3[:] *= 1./norm(u3)
    e1[:] = cross(u2,u3)
    e2[:] = cross(u3,u1)

    # get the indices of the CPs at the corner of the plane
    cp1_idx = map.cp_idx[:,cp1_jkm[1],cp1_jkm[2],cp1_jkm[3]]
    cp2_idx = map.cp_idx[:,cp2_jkm[1],cp2_jkm[2],cp2_jkm[3]]
    cp3_idx = map.cp_idx[:,cp3_jkm[1],cp3_jkm[2],cp3_jkm[3]]

    # loop over the CPs that lie in the jdi-th plane
    for jit1 = 1:map.nctl[it1]
      jkm[it1] = jit1
      for jit2 = 1:map.nctl[it2]
        jkm[it2] = jit2
        if jkm == cp1_jkm || jkm == cp2_jkm || jkm == cp3_jkm
          # the corner CPs are independent, so skip them
          continue
        end
        # get the dependency coefficients
        x[:] = map.cp_xyz[:,jkm[1],jkm[2],jkm[3]] - cp1_xyz
        u1_coeff, u2_coeff = calcObliqueProjection(x, u1, u2, e1, e2)
        coeff2 = u1_coeff*fac1
        coeff3 = u2_coeff*fac2
        coeff1 = -coeff2 - coeff3

        # insert into sparse storage
        for bdi = 1:3
          # set data for dependent CP coord
          iLfun[ptr] = fncidx
          jLvar[ptr] = map.cp_idx[bdi,jkm[1],jkm[2],jkm[3]]
          LinG[ptr] = 1.0
          ptr += 1
          # set data for independent CP1
          iLfun[ptr] = fncidx
          jLvar[ptr] = cp1_idx[bdi]
          LinG[ptr] = -coeff1
          ptr += 1
          # set data for independent CP2
          iLfun[ptr] = fncidx
          jLvar[ptr] = cp2_idx[bdi]
          LinG[ptr] = -coeff2
          ptr += 1
          # set data for independent CP3
          iLfun[ptr] = fncidx
          jLvar[ptr] = cp3_idx[bdi]
          LinG[ptr] = -coeff3
          ptr += 1
          Flow[fncidx] = 0.0
          Fupp[fncidx] = 0.0
          fncidx += 1
        end

      end # jit2 loop
    end # jit1 loop
  end # jdi loop
  return fncidx, ptr
end

@doc """
### calcObliqueProjection

Basis coefficients for u1 and u2 based on an oblique projection along e1, e2.

**Arguments**

* `x` : vector whose representation in (`u1`,`u2`) is sought
* `u1`,`u2` : unit basis vectors, not necessarily orthogonal
* `e1`,`e2` : covariant basis vectors satisfying dot(e1,u2) = 0, dot(e2,u1) = 0

"""->
function calcObliqueProjection(x::AbstractArray{T,1}, u1::AbstractArray{T,1},
                               u2::AbstractArray{T,1}, e1::AbstractArray{T,1},
                               e2::AbstractArray{T,1}) where T
  return dot(e1,x)/dot(e1,u1), dot(e2,x)/dot(e2,u2)
end

@doc """
### numLinearCornerConstraints

Returns the number of corner constraints (see setLinearCornerConstraints! for
further details).

**Arguments**

* `map` : PumiMapping object
* `di` : the index the defines the constrained plane (e.g. `di`=2 -> j planes)

"""->
function numLinearCornerConstraints(map::PumiMapping{Tffd}, di::Int) where Tffd
  return map.nctl[di]*2
end

@doc """
### countVarsLinearCornerConstraints!

Updates array `cntEq` with the number of variables involved in each linear-corner
constraint.

**Arguments**

* `map` : PumiMapping object
* `di` : the index the defines the constrained plane (e.g. `di`=2 -> j planes)
* `cntEq` : array with number of variables involved in each equality constraint
* `fncidx` : index to the next equality constraint to add to cntEq

**Returns** `fncidx` + number of linear-plane constraints + 1

"""->
function countVarsLinearCornerConstraints!(map::PumiMapping{Tffd}, di::Int,
                                           cntEq::AbstractArray{Int,1},
                                           fncidx::Int) where Tffd
  num_cnstr = numLinearCornerConstraints(map, di)
  cntEq[fncidx:fncidx+num_cnstr-1] = 4
  return fncidx+num_cnstr
end

@doc """
### setLinearCornerConstraints!

This constraint is used to keep each plane in the `di` direction square; this
should be used in conjuction with setLinearPlaneConstraints!.  It also assumes
that the di coordinates of the CPs are the same.

**Arguments**

* `map`  : PumiMapping object
* `di`   : the index the defines the constrained plane (e.g. `di`=2 -> j planes)
* `iLfun`, `jLvar` : the k^th element in the gradient matrix has row index
*                    `iLfun(k)` and column index `jLvar(k)`
* `LinG` : `LinG(k) = A(iLfun(k),jLvar(k))`
* `Flow`, `Fupp` : the upper and lower bounds on all constraints
* `fncidx` : function counter; next function to use in iLfun
* `ptr` : index pointer to next free space in iLfun and jLvar

**Returns** `fncidx` + number of linear-plane constraints + 1, and `ptr` + total
  number of non-zeros in the sparse contraints + 1

"""->
function setLinearCornerConstraints!(map::PumiMapping{Tffd}, di::Int, iLfun::AbstractArray{Int,1},
    jLvar::AbstractArray{Int,1}, LinG::AbstractArray{Tffd,1},
    Flow::AbstractArray{Tffd,1}, Fupp::AbstractArray{Tffd,1},
    fncidx::Int, ptr::Int) where Tffd

  it1 = mod(di,3)+1
  it2 = mod(di+1,3)+1

  jkm = zeros(Int,3)
  cp1_jkm = zeros(jkm); cp2_jkm = zeros(jkm); cp3_jkm = zeros(jkm)
  cp1_idx = zeros(jkm); cp2_idx = zeros(jkm); cp3_idx = zeros(jkm)
  cp1_xyz = zeros(3); cp2_xyz = zeros(3); cp3_xyz = zeros(3)
  u1 = zeros(3); tvec = zeros(3)

  # loop over the CP planes in the direction di
  for jdi = 1:map.nctl[di]
    jkm[di] = jdi

    # set the (j,k,m) indices for the CPs at the corners of the plane
    cp1_jkm[di] = jdi
    cp1_jkm[it1] = 1
    cp1_jkm[it2] = 1

    cp2_jkm[di] = jdi
    cp2_jkm[it1] = 1
    cp2_jkm[it2] = map.nctl[it2]

    cp3_jkm[di] = jdi
    cp3_jkm[it1] = map.nctl[it1]
    cp3_jkm[it2] = 1

    # get the coordinates of the CPs at the corners of the plane, and some unit
    # vectors
    cp1_xyz = map.cp_xyz[:,cp1_jkm[1],cp1_jkm[2],cp1_jkm[3]]
    cp2_xyz = map.cp_xyz[:,cp2_jkm[1],cp2_jkm[2],cp2_jkm[3]]
    cp3_xyz = map.cp_xyz[:,cp3_jkm[1],cp3_jkm[2],cp3_jkm[3]]

    u1[:] = cp2_xyz - cp1_xyz
    fac1 = 1./norm(u1)
    u1[:] *= fac1

    tvec[it1] = -u1[it2]
    tvec[it2] = u1[it1]
    tvec[di] = 0.0
    alpha = dot(tvec, cp3_xyz - cp1_xyz)

    # u1 and tvec should be orthogonal
    @assert( dot(u1,tvec) < 1e-14 )

    # First equation, for coordinate it1
    # x3[it1] = -alpha*fac1*(x2[it2] - x1[it2]) + x1[it1]

    # set data for dependent cp3
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp3_idx[it1]
    LinG[ptr] = 1.0
    ptr += 1

    # set data for independent cp1
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp1_idx[it1]
    LinG[ptr] = -1.0
    ptr += 1
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp1_idx[it2]
    LinG[ptr] = -alpha*fac1
    ptr += 1

    # set data for independent cp2
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp2_idx[it2]
    LinG[ptr] = alpha*fac1
    ptr += 1

    Flow[fncidx] = 0.0
    Fupp[fncidx] = 0.0
    fncidx += 1

    # Second equation, for coordinate it2
    # x3[it2] = alpha*fac1*(x2[it1] - x1[it1]) + x1[it2]

    # set data for dependent cp3
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp3_idx[it2]
    LinG[ptr] = 1.0
    ptr += 1

    # set data for independent cp1
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp1_idx[it2]
    LinG[ptr] = -1.0
    ptr += 1
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp1_idx[it1]
    LinG[ptr] = alpha*fac1
    ptr += 1

    # set data for independent cp2
    iLfun[ptr] = fncidx
    jLvar[ptr] = cp2_idx[it1]
    LinG[ptr] = -alpha*fac1
    ptr += 1

    Flow[fncidx] = 0.0
    Fupp[fncidx] = 0.0
    fncidx += 1
  end

  return fncidx, ptr
end

@doc """
### numLinearStretchConstraints

Returns the number of stretch constraints (see setLinearStretchConstraints! for
further details).

**Arguments**

* `map` : PumiMapping object

"""->
function numLinearStretchConstraints(map::PumiMapping{Tffd}) where Tffd
  return map.nctl[1]*map.nctl[2]*map.nctl[3] - 1
end

@doc """
### countVarsLinearStretchConstraints!

Updates array `cntEq` with the number of variables involved in each
linear-stretch constraint.

**Arguments**

* `map` : PumiMapping object
* `cntEq` : array with number of variables involved in each equality constraint
* `fncidx` : index to the next equality constraint to add to cntEq

**Returns** `fncidx` + number of linear-plane constraints + 1

"""->
function countVarsLinearStretchConstraints!(map::PumiMapping{Tffd},
                                            cntEq::AbstractArray{Int,1},
                                            fncidx::Int) where Tffd
  num_cnstr = numLinearStretchConstraints(map)
  cntEq[fncidx:fncidx+num_cnstr-1] = 2
  return fncidx+num_cnstr
end

@doc """
### setLinearStretchConstraints!

This constraint is used to scale all the control points' `di` coordinate based
on a single master control point.  Useful when span is permitted to change and
you want each CP in a plane to stretch in proportion.

**Arguments**

* `map`  : PumiMapping object
* `di`   : the index the defines the constrained plane (e.g. `di`=2 -> j planes)
* `iLfun`, `jLvar` : the k^th element in the gradient matrix has row index
*                    `iLfun(k)` and column index `jLvar(k)`
* `LinG` : `LinG(k) = A(iLfun(k),jLvar(k))`
* `Flow`, `Fupp` : the upper and lower bounds on all constraints
* `fncidx` : function counter; next function to use in iLfun
* `ptr` : index pointer to next free space in iLfun and jLvar

**Returns** `fncidx` + number of linear-plane constraints + 1, and `ptr` + total
  number of non-zeros in the sparse contraints + 1

"""->
function setLinearStretchConstraints!(map::PumiMapping{Tffd}, di::Int, iLfun::AbstractArray{Int,1},
    jLvar::AbstractArray{Int,1}, LinG::AbstractArray{Tffd,1},
    Flow::AbstractArray{Tffd,1}, Fupp::AbstractArray{Tffd,1},
    fncidx::Int, ptr::Int) where Tffd

  # get master CP index and xyz
  master_jkm = ones(Int,3)
  master_jkm[di] = map.nctl[di]
  master_idx = map.cp_idx[di,master_jkm[1],master_jkm[2],master_jkm[3]]
  master_xyz = map.cp_xyz[di,master_jkm[1],master_jkm[2],master_jkm[3]]

  for k = 1:map.nctl[3]
    for j = 1:map.nctl[2]
      for i = 1:map.nctl[1]
        if i == master_jkm[1] && j == master_jkm[2] && k == master_jkm[3]
          # if this is the master node, skip it
          continue
        end

        # set data for dependent CP
        iLfun[ptr] = fncidx
        jLvar[ptr] = map.cp_idx[di,i,j,k]
        LinG[ptr] = 1.0
        ptr += 1

        # set data for master CP
        iLfun[ptr] = fncidx
        jLvar[ptr] = master_idx
        LinG[ptr] = map.cp_xyz[di,i,j,k]/master_xyz
        ptr += 1

        Flow[fncidx] = 0.0
        Fupp[fncidx] = 0.0
        fncidx += 1
      end
    end
  end
  return fncidx, ptr
end

@doc """
### numLinearRootConstraints

Returns the number of root constraints (see setLinearRootConstraints! for
further details).

**Arguments**

* `map` : PumiMapping object

"""->
function numLinearRootConstraints(map::PumiMapping{Tffd}) where Tffd
  return 5
end

@doc """
### countVarsLinearRootConstraints!

Updates array `cntEq` with the number of variables involved in each
linear-root constraint.

**Arguments**

* `map` : PumiMapping object
* `cntEq` : array with number of variables involved in each equality constraint
* `fncidx` : index to the next equality constraint to add to cntEq

**Returns** `fncidx` + number of linear-plane constraints + 1

"""->
function countVarsLinearRootConstraints!(map::PumiMapping{Tffd},
                                         cntEq::AbstractArray{Int,1},
                                         fncidx::Int) where Tffd
  for bdi = 1:3
    cntEq[fncidx] = 1
    fncidx += 1
  end
  cntEq[fncidx] = 2
  fncidx += 1
  cntEq[fncidx] = 2
  fncidx += 1
  return fncidx
end

@doc """
### setLinearRootConstraints!

This constraint sets some very specific linear constraints at the root of the
FFD, i.e. where the FFD touches the symmetric plane.  One constraint freezes one
CP completely, and another CP is keep aligned (in some direction) with the
frozen CP.  **Assumes the symmetry plane is at the low end of the CP indices**.

**Arguments**

* `map`  : PumiMapping object
* `di`   : defines the symmetry plane (e.g. `di`=2 -> y plane)
* `iLfun`, `jLvar` : the k^th element in the gradient matrix has row index
*                    `iLfun(k)` and column index `jLvar(k)`
* `LinG` : `LinG(k) = A(iLfun(k),jLvar(k))`
* `Flow`, `Fupp` : the upper and lower bounds on all constraints
* `fncidx` : function counter; next function to use in iLfun
* `ptr` : index pointer to next free space in iLfun and jLvar

**Returns** `fncidx` + number of linear-plane constraints + 1, and `ptr` + total
  number of non-zeros in the sparse contraints + 1

"""->
function setLinearRootConstraints!(map::PumiMapping{Tffd}, di::Int, iLfun::AbstractArray{Int,1},
    jLvar::AbstractArray{Int,1}, LinG::AbstractArray{Tffd,1},
    Flow::AbstractArray{Tffd,1}, Fupp::AbstractArray{Tffd,1},
    fncidx::Int, ptr::Int) where Tffd

  # frozen CP constraint
  for bdi = 1:3
    iLfun[ptr] = fncidx
    jLvar[ptr] = map.cp_idx[bdi,1,1,1]
    LinG[ptr] = 1.0
    ptr += 1

    Flow[fncidx] = map.cp_xyz[bdi,1,1,1]
    Fupp[fncidx] = Flow[fncidx]
    fncidx += 1
  end

  # Aligned CP: it is allowed to move in the direction it1 only
  it1 = mod(di,3)+1
  it2 = mod(di+1,3)+1
  jkm = ones(Int,3)
  jkm[it1] = map.nctl[it1]

  # constraint x2[di] = x1[di]
  iLfun[ptr] = fncidx
  jLvar[ptr] = map.cp_idx[di,jkm[1],jkm[2],jkm[3]]
  LinG[ptr] = 1.0
  ptr += 1
  iLfun[ptr] = fncidx
  jLvar[ptr] = map.cp_idx[di,1,1,1]
  LinG[ptr] = -1.0
  ptr += 1

  Flow[fncidx] = 0.0
  Fupp[fncidx] = 0.0
  fncidx += 1

  # constraint x2[it2] = x1[it2]
  iLfun[ptr] = fncidx
  jLvar[ptr] = map.cp_idx[it2,jkm[1],jkm[2],jkm[3]]
  LinG[ptr] = 1.0
  ptr += 1
  iLfun[ptr] = fncidx
  jLvar[ptr] = map.cp_idx[it2,1,1,1]
  LinG[ptr] = -1.0
  ptr += 1

  Flow[fncidx] = 0.0
  Fupp[fncidx] = 0.0
  fncidx += 1

  return fncidx, ptr
end


text:
relaxed_infeas_tols = (tol_infeas = 1e-9,)
insts = OrderedDict()
insts["minimal"] = [
    ((0.7, 4, 1e-3, true, true),),
    ((0.7, 4, 1e-3, false, true), :SOCExpPSD),
    ((2.0, 4, 1e-3, true, false), nothing, relaxed_infeas_tols),
    ((2.0, 4, 1e-3, false, false), :SOCExpPSD, relaxed_infeas_tols),
    ]
insts["fast"] = [
    ((1.0, 2, 1e-3, true, false),),
    ((1.0, 2, 1e-3, false, false), :SOCExpPSD),
    ((4.0, 2, 1e-3, true, false),),
    ((4.0, 2, 1e-3, false, false), :SOCExpPSD),
    ]
insts["various"] = vcat(insts["minimal"], insts["fast"])
return (ContractionJuMP, insts)

ext: jl


text:
x = 10;
if (x == 10)
    println("x0 == ",x);
elseif (x == 20)
    println("x1 == ",20);
else
    println("didnot match");
end;

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
# COP functions
function violationsSum(g::Vector, h::Vector; ε =0.0)
    sum_g = sum(max.(0.0, g))
    sum_h = 0.0

    for i = 1:length(h)
        if !isapprox(h[i], 0.0, atol=ε)
            sum_h += abs(h[i])
        end
    end

    return (sum_g/length(g)) + (sum_h / length(h))
end


text:
# Slave part of the Algorithm
using Random

include("..\\elements\\Elements.jl")
include("..\\file_input\\FileLoaders.jl")
include("..\\configurations\\MasterConfiguration.jl")
include("..\\instances\\SlaveInstance.jl")

mutable struct ODPSlaveAlgorithm
    # Slave instance
    instance::SlaveInstance
    # Heap of candidate points sent to the Master for the next iteration
    nextCand::Heap
    # Number of distances calculated in the current iteration
    totDistCount::Int64
    # Number of active points left after the current iteration
    currentActive::Int64

    ODPSlaveAlgorithm() = new()
end

# Factory method to create an ODPSlaveAlgorithm instance
function getSlaveAlgorithm(path::String)
    alg = ODPSlaveAlgorithm()
    alg.instance = createSlaveInstance(path)
    alg.nextCand = MinHeap()
    alg.totDistCount = 0
    alg.currentActive = 0
    return alg
end

# Some dataset info methods
function getDataInfo(this::ODPSlaveAlgorithm)
    info = DataInfo(this.instance.idHeader, getDatasetDimension(this.instance), getPointDimension(this.instance))
    return info
end

function getDataset(this::ODPSlaveAlgorithm)
    return this.instance.dataset
end

# Invoked by the Master to send configuration parameters to the Slave.
# Set ODPSlaveAlgorithm configuration parameters and load dataset.
function setParam(this::ODPSlaveAlgorithm, n::Int64, k::Int64, m::Int64, dist::String)
    this.instance.n = n
    this.instance.k = k
    this.instance.m = m
    this.instance.dist = DistanceFunction(dist)
    # Create a dataset from file in path
    loadDataset(this.instance)
    # Initializing Points NN with Max k elements
    for p in this.instance.dataset
        p.nn = NN(this.instance.k)
        # filling plotDataset
        # push!(plotDataset, p.coords)
    end

    # Retrieve active points
    this.currentActive = getDatasetDimension(this.instance)
end

# Invoked by the Master.
# Chose an initial set of candidate Points, chosen at random from dataset
function initCandidates(this::ODPSlaveAlgorithm)
    count = this.instance.m
    collection = this.instance.dataset
    res = Set{Point}()
    rng = MersenneTwister()
    i = 0

    while i < count
        randElement = rand(rng, collection, 1)[1]
        if !(randElement in res)
            push!(res, randElement)
            i+=1
        end
    end

    return res
end

function cmpCandCand(this::ODPSlaveAlgorithm, C::Set{Point})
    distCount = 0
    tmpC = Set{Point}()
    for el in C
        push!(tmpC, el)
        # remove point from dataset
        setdiff!(this.instance.dataset, [el])
        this.currentActive-=1
    end

    cut = NaN
    for p in C
        for q in tmpC
            dist = 0.
            if p.id == q.id
                cut = q
                distCount+=1
            else
                dist = euclidean(p.coords, q.coords)
                distCount+=1
            end
            updateMin(p.nn, dist)

            if p.id != q.id
                updateMin(q.nn, dist)
            end
        end
        setdiff!(tmpC, [cut])
    end
    return distCount
end

function cmpDatasetCand(this::ODPSlaveAlgorithm, C::Set{Point}, lb::Float64)
    distCount = 0
    d = .0

    for p in this.instance.dataset
        for q in C
            if max(p.nn.weight, q.nn.weight) >= lb
                d = euclidean(p.coords, q.coords)
                distCount+=1

                updateMin(p.nn, d)
                updateMin(q.nn, d)
            end
        end

        if p.active && (p.nn.weight < lb)
            p.active = false
            this.currentActive-=1
        end

        if p.nn.weight > lb
            updateMax(this.nextCand, p)
        end
    end
    return distCount
end

# Invoked by the Master.
# Execute an algorithm iteration. Input: Heap of Candidate Points, solvingSet lower bound
function computeIteration(this::ODPSlaveAlgorithm, C::Set{Point}, lb::Float64)
    distCount = 0

    # First cycle -> Compare Cnod with himself and compute weights
    distCount += cmpCandCand(this, C)

    # Preparing candidate heap for next iteration
    this.nextCand = MinHeap(this.instance.m)

    # Second cycle -> Compare Dataset elements with Cnod, filling nextCand (candidate points for next iteration)
    distCount += cmpDatasetCand(this, C, lb)

    this.totDistCount += distCount

    # Contains all the information sent by the Slave to the Master at the end of an iteration.
    iterationResult = IterationResult(this.nextCand, this.totDistCount, getDatasetDimension(this.instance), this.currentActive, 0)

    return iterationResult
end


text:
"""
SPADL
"""
Base.@kwdef mutable struct RegularSPADL
    action_id::Union{Nothing,Int} = nothing
    bodypart_id::Union{Nothing,Int} = nothing
    bodypart_name::Union{Nothing,String} = nothing
    end_x::Union{Nothing,Float64} = nothing
    end_y::Union{Nothing,Float64} = nothing
    game_id::Union{Nothing,Int} = nothing
    original_event_id::Union{Nothing,Int} = nothing
    period_id::Union{Nothing,String} = nothing
    result_id::Union{Nothing,Int} = nothing
    result_name::Union{Nothing,String} = nothing
    start_x::Union{Nothing,Float64} = nothing
    start_y::Union{Nothing,Float64} = nothing
    team_id::Union{Nothing,Int} = nothing
    time_seconds::Union{Nothing,Float64} = nothing
    type_id::Union{Nothing,Int} = nothing
    type_name::Union{Nothing,String} = nothing
end

"""
WyscoutEventTags
"""
Base.@kwdef mutable struct WyscoutEventTags
    goal::Bool = false
    own_goal::Bool = false
    assist::Bool = false
    key_pass::Bool = false
    counter_attack::Bool = false
    left_foot::Bool = false
    right_foot::Bool = false
    head_body::Bool = false
    direct::Bool = false
    indirect::Bool = false
    dangerous_ball_lost::Bool = false
    blocked::Bool = false
    high::Bool = false
    low::Bool = false
    interception::Bool = false
    clearance::Bool = false
    opportunity::Bool = false
    feint::Bool = false
    missed_ball::Bool = false
    free_space_right::Bool = false
    free_space_left::Bool = false
    take_on_left::Bool = false
    take_on_right::Bool = false
    sliding_tackle::Bool = false
    anticipated::Bool = false
    anticipation::Bool = false
    red_card::Bool = false
    yellow_card::Bool = false
    second_yellow_card::Bool = false
    position_goal_low_center::Bool = false
    position_goal_low_right::Bool = false
    position_goal_mid_center::Bool = false
    position_goal_mid_left::Bool = false
    position_goal_low_left::Bool = false
    position_goal_mid_right::Bool = false
    position_goal_high_center::Bool = false
    position_goal_high_left::Bool = false
    position_goal_high_right::Bool = false
    position_out_low_right::Bool = false
    position_out_mid_left::Bool = false
    position_out_low_left::Bool = false
    position_out_mid_right::Bool = false
    position_out_high_center::Bool = false
    position_out_high_left::Bool = false
    position_out_high_right::Bool = false
    position_post_low_right::Bool = false
    position_post_mid_left::Bool = false
    position_post_low_left::Bool = false
    position_post_mid_right::Bool = false
    position_post_high_center::Bool = false
    position_post_high_left::Bool = false
    position_post_high_right::Bool = false
    through::Bool = false
    fairplay::Bool = false
    lost::Bool = false
    neutral::Bool = false
    won::Bool = false
    accurate::Bool = false
    not_accurate::Bool = false
    # added offside
    offside::Bool = false
end

"""
WyscoutEventFixed
"""
Base.@kwdef mutable struct WyscoutEventFixed
    event_id::Int
    game_id::Int
    player_id::Int
    period_id::String
    start_x::Union{Nothing,Float64} = nothing
    start_y::Union{Nothing,Float64} = nothing
    end_x::Union{Nothing,Float64} = nothing
    end_y::Union{Nothing,Float64} = nothing
    team_id::Union{Nothing,Int} = nothing
    type_id::Union{Nothing,Int} = nothing
    subtype_id::Union{Nothing,Int} = nothing
    milliseconds::Float64
end

"""
    WyscoutData
"""
Base.@kwdef mutable struct WyscoutData
    event::WyscoutEvent
    tags::Union{Nothing,WyscoutEventTags} = nothing
    event_fixed::Union{Nothing,WyscoutEventFixed} = nothing
end

"""
    ActionsConfig
"""
Base.@kwdef mutable struct ActionsConfig
    field_length::Float64 = 105.0 # unit: meters
    field_width::Float64 = 68.0  # unit meters
    bodyparts::Vector{String} = ["foot", "head", "other", "head/other"]
    results::Vector{String} = ["fail","success","offside","owngoal","yellow_card","red_card"]
    actionstypes::Vector{String} = [
        "pass",
        "cross",
        "throw_in",
        "freekick_crossed",
        "freekick_short",
        "corner_crossed",
        "corner_short",
        "take_on",
        "foul",
        "tackle",
        "interception",
        "shot",
        "shot_penalty",
        "shot_freekick",
        "keeper_save",
        "keeper_claim",
        "keeper_punch",
        "keeper_pick_up",
        "clearance",
        "bad_touch",
        "non_action",
        "dribble",
        "goalkick"
    ]
end


text:
module inputs

using Pkg
using Pipe: @pipe

cd(dirname(@__DIR__)) # move pwd up one level from this file
# @info "The PWD for inputs module is: $(pwd())"
# @info "The __DIR__ inputs module is: $(@__DIR__)"

################
#    Day 1
################
ex_d1() = [1721, 979, 366, 299, 675, 1456]
input_d1() = readlines("input/day01.txt") .|> x->parse(Int, x)

################
#    Day 2
################
ex_d2() = @pipe "1-3 a: abcde
1-3 b: cdefg
2-9 c: ccccccccc" |> String.(split(_, "\n"))
input_d2() = @pipe read("input/day02.txt", String) |> String.(split(_, "\n"))

################
#    Day 3
################
ex_d3() = """..##.......
#...#...#..
.#....#..#.
..#.#...#.#
.#...##..#.
..#.##.....
.#.#.#....#
.#........#
#.##...#...
#...##....#
.#..#...#.#"""
input_d3() = read("input/day03.txt", String);

################
#    Day 4
################
raw_ex_d4 = """ecl:gry pid:860033327 eyr:2020 hcl:#fffffd
byr:1937 iyr:2017 cid:147 hgt:183cm

iyr:2013 ecl:amb cid:350 eyr:2023 pid:028048884
hcl:#cfa07d byr:1929

hcl:#ae17e1 iyr:2013
eyr:2024
ecl:brn pid:760753108 byr:1931
hgt:179cm

hcl:#cfa07d eyr:2025 pid:166559648
iyr:2011 ecl:brn hgt:59in"""
ex_d4() = @pipe String.(split(raw_ex_d4, "\n")) |> # equiv to output of readlines(input.txt)
            join(_, "\n") |>
            split(_, "\n\n") |>
            split.(_, r"\s+")

input_d4() = @pipe readlines("input/day04.txt") |>
            join(_, "\n") |>
            split(_, "\n\n") |>
            split.(_, r"\s+")
invalid_raw = @pipe """eyr:1972 cid:100
hcl:#18171d ecl:amb hgt:170 pid:186cm iyr:2018 byr:1926

iyr:2019
hcl:#602927 eyr:1967 hgt:170cm
ecl:grn pid:012533040 byr:1946

hcl:dab227 iyr:2012
ecl:brn hgt:182cm pid:021572410 eyr:2020 byr:1992 cid:277

hgt:59cm ecl:zzz
eyr:2038 hcl:74454a iyr:2023
pid:3556412378 byr:2007"""
invalid_d4() = @pipe String.(split(invalid_raw, "\n")) |> # equiv to output of readlines(input.txt)
            join(_, "\n") |>
            split(_, "\n\n") |>
            split.(_, r"\s+")

valid_raw = """pid:087499704 hgt:74in ecl:grn iyr:2012 eyr:2030 byr:1980
hcl:#623a2f

eyr:2029 ecl:blu cid:129 byr:1989
iyr:2014 pid:896056539 hcl:#a97842 hgt:165cm

hcl:#888785
hgt:164cm byr:2001 iyr:2015 cid:88
pid:545766238 ecl:hzl
eyr:2022

iyr:2010 hgt:158cm hcl:#b6652a ecl:blu byr:1944 eyr:2021 pid:093154719"""
valid_d4() = @pipe String.(split(valid_raw, "\n")) |> # equiv to output of readlines(input.txt)
            join(_, "\n") |>
            split(_, "\n\n") |>
            split.(_, r"\s+")

################
#    Day 5
################
input_d5() = readlines("input/day05.txt")

################
#    Day 6
################
ex1_d6() = @pipe """abcx
abcy
abcz""" |> String.(split(_, "\n")) |> join(_, "\n") |> split(_, "\n\n") |> split.(_, r"\s+")
ex2_d6() = @pipe """abc

a
b
c

ab
ac

a
a
a
a

b""" |> String.(split(_, "\n")) |> join(_, "\n") |> split(_, "\n\n") |> split.(_, r"\s+")
input_d6() = @pipe readlines("input/day06.txt") |> join(_,"\n") |> split(_, "\n\n") |> split.(_, r"\s+")

################
#    Day 7
################

ex1_d7() = @pipe """light red bags contain 1 bright white bag, 2 muted yellow bags.
dark orange bags contain 3 bright white bags, 4 muted yellow bags.
bright white bags contain 1 shiny gold bag.
muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.
shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.
dark olive bags contain 3 faded blue bags, 4 dotted black bags.
vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.
faded blue bags contain no other bags.
dotted black bags contain no other bags.""" |> String.(split(_, "\n"))
input_d7() = readlines("input/day07.txt")

################
#    Day 8
################
ex1_d8() = @pipe """nop +0
acc +1
jmp +4
acc +3
jmp -3
acc -99
acc +1
jmp -4
acc +6""" |> String.(split(_, "\n"))
input_d8() = readlines("input/day08.txt")

################
#    Day 9
################
ex1_d9() = @pipe """35
20
15
25
47
40
62
55
65
95
102
117
150
182
127
219
299
277
309
576""" |> String.(split(_, "\n")) |> parse.(Int, _)
input_d9() = parse.(Int, readlines("input/day09.txt"))

################
#    Day 10
################

ex0_d10() = @pipe """16
10
15
5
1
11
7
19
6
12
4""" |> String.(split(_, "\n")) |> parse.(Int, _)

ex1_d10() = @pipe """28
33
18
42
31
14
46
20
48
47
24
23
49
45
19
38
39
11
1
32
25
35
8
17
7
9
4
2
34
10
3""" |> String.(split(_, "\n")) |> parse.(Int, _)
input_d10() = parse.(Int, readlines("input/day10.txt"))

################
#    Day 11
################

ex_d11() = @pipe """L.LL.LL.LL
LLLLLLL.LL
L.L.L..L..
LLLL.LL.LL
L.LL.LL.LL
L.LLLLL.LL
..L.L.....
LLLLLLLLLL
L.LLLLLL.L
L.LLLLL.LL""" |> String.(split(_, "\n"))

ex1_d11() = @pipe """.......#.
...#.....
.#.......
.........
..#L....#
....#....
.........
#........
...#.....""" |> String.(split(_, "\n"))

ex2_d11() = @pipe """.##.##.
#.#.#.#
##...##
...L...
##...##
#.#.#.#
.##.##.""" |> String.(split(_, "\n"))

input_d11() = readlines("input/day11.txt")

################
#    Day 12
################

ex_d12() = @pipe """F10
N3
F7
R90
F11""" |> String.(split(_, "\n"))

input_d12() = readlines("input/day12.txt")

################
#    Day 13
################

ex_d13() = @pipe """939
7,13,x,x,59,x,31,19""" |> String.(split(_, "\n"))

input_d13() = @pipe """1000066
13,x,x,41,x,x,x,37,x,x,x,x,x,659,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,19,x,x,x,23,x,x,x,x,x,29,x,409,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,x,17""" |> String.(split(_, "\n"))

################
#    Day 14
################

ex1_d14() = """mask = XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X
mem[8] = 11
mem[7] = 101
mem[8] = 0"""

ex2_d14() = """mask = 000000000000000000000000000000X1001X
mem[42] = 100
mask = 00000000000000000000000000000000X0XX
mem[26] = 1"""

input_d14() = read("input/day14.txt", String)

################
#    Day 16
################

ex_d16() = """class: 1-3 or 5-7
row: 6-11 or 33-44
seat: 13-40 or 45-50

your ticket:
7,1,14

nearby tickets:
7,3,47
40,4,50
55,2,20
38,6,12"""

ex2_d16() = """class: 0-1 or 4-19
row: 0-5 or 8-19
seat: 0-13 or 16-19

your ticket:
11,12,13

nearby tickets:
3,9,18
15,1,5
5,14,9"""

input_d16() = read("input/day16.txt", String)

################
#    Day 17
################

ex_d17() = """.#.
..#
###"""

input_d17() = """##..#.#.
#####.##
#######.
#..#..#.
#.#...##
..#....#
....#..#
..##.#.."""

################
#    Day 18
################

ex1_d18() = "1 + 2 * 3 + 4 * 5 + 6"
ex2_d18() = "1 + (2 * 3) + (4 * (5 + 6))"
ex3_d18() = "2 * 3 + (4 * 5)"
ex4_d18() = "5 + (8 * 3 + 9 + 3 * 4 * 3)"
ex5_d18() = "5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))"
ex6_d18() = "((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2"
input_d18() = readlines("input/day18.txt")

################
#    Day 19
################

ex_d19() = """0: 4 1 5
1: 2 3 | 3 2
2: 4 4 | 5 5
3: 4 5 | 5 4
4: "a"
5: "b"

ababbb
bababa
abbbab
aaabbb
aaaabbb"""

input_d19() = read("input/day19.txt", String)

################
#    Day 20
################

ex_d20() = """Tile 2311:
..##.#..#.
##..#.....
#...##..#.
####.#...#
##.##.###.
##...#.###
.#.#.#..##
..#....#..
###...#.#.
..###..###

Tile 1951:
#.##...##.
#.####...#
.....#..##
#...######
.##.#....#
.###.#####
###.##.##.
.###....#.
..#.#..#.#
#...##.#..

Tile 1171:
####...##.
#..##.#..#
##.#..#.#.
.###.####.
..###.####
.##....##.
.#...####.
#.##.####.
####..#...
.....##...

Tile 1427:
###.##.#..
.#..#.##..
.#.##.#..#
#.#.#.##.#
....#...##
...##..##.
...#.#####
.#.####.#.
..#..###.#
..##.#..#.

Tile 1489:
##.#.#....
..##...#..
.##..##...
..#...#...
#####...#.
#..#.#.#.#
...#.#.#..
##.#...##.
..##.##.##
###.##.#..

Tile 2473:
#....####.
#..#.##...
#.##..#...
######.#.#
.#...#.#.#
.#########
.###.#..#.
########.#
##...##.#.
..###.#.#.

Tile 2971:
..#.#....#
#...###...
#.#.###...
##.##..#..
.#####..##
.#..####.#
#..#.#..#.
..####.###
..#.#.###.
...#.#.#.#

Tile 2729:
...#.#.#.#
####.#....
..#.#.....
....#..#.#
.##..##.#.
.#.####...
####.#.#..
##.####...
##..#.##..
#.##...##.

Tile 3079:
#.#.#####.
.#..######
..#.......
######....
####.#..#.
.#...#.##.
#.#####.##
..#.###...
..#.......
..#.###..."""

input_d20() = read("input/day20.txt", String)

################
#    Day 21
################

ex_d21() = """mxmxvkd kfcds sqjhc nhms (contains dairy, fish)
trh fvjkl sbzzf mxmxvkd (contains dairy)
sqjhc fvjkl (contains soy)
sqjhc mxmxvkd sbzzf (contains fish)"""

input_d21() = read("input/day21.txt", String)

################
#    Day 22
################

ex_d22() = """Player 1:
9
2
6
3
1

Player 2:
5
8
4
7
10"""

input_d22() = """Player 1:
8
19
46
11
36
10
35
9
24
22
50
1
34
7
18
28
3
38
43
2
6
42
23
12
20

Player 2:
39
27
44
29
5
48
30
32
15
31
14
21
49
17
45
47
16
26
33
25
13
41
4
40
37"""

################
#    Day 24
################

ex_d24() = """sesenwnenenewseeswwswswwnenewsewsw
neeenesenwnwwswnenewnwwsewnenwseswesw
seswneswswsenwwnwse
nwnwneseeswswnenewneswwnewseswneseene
swweswneswnenwsewnwneneseenw
eesenwseswswnenwswnwnwsewwnwsene
sewnenenenesenwsewnenwwwse
wenwwweseeeweswwwnwwe
wsweesenenewnwwnwsenewsenwwsesesenwne
neeswseenwwswnwswswnw
nenwswwsewswnenenewsenwsenwnesesenew
enewnwewneswsewnwswenweswnenwsenwsw
sweneswneswneneenwnewenewwneswswnese
swwesenesewenwneswnwwneseswwne
enesenwswwswneneswsenwnewswseenwsese
wnwnesenesenenwwnenwsewesewsesesew
nenewswnwewswnenesenwnesewesw
eneswnwswnwsenenwnwnwwseeswneewsenese
neswnwewnwnwseenwseesewsenwsweewe
wseweeenwnesenwwwswnew"""

input_d24() = read("input/day24.txt", String)

end


text:
module CouenneNL

export CouenneNLSolver

using Lazy: @forward
using AmplNLWriter
using AmplNLWriter: AmplNLMathProgModel, AmplNLNonlinearModel, AmplNLLinearQuadraticModel
using MathProgBase
using MathProgBase.SolverInterface

const depsjl = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsjl)
    error("CouenneNL not properly installed. Please run\nPkg.build(\"CouenneNL\")")
else
    include(depsjl)
end

struct CouenneNLSolver <: AbstractMathProgSolver
    solver::AmplNLSolver
end

function CouenneNLSolver(options::Vector{String}=String[], filename::String="")
    CouenneNLSolver(AmplNLSolver(couenne_executable, options, filename))
end

struct CouenneNLMathProgModel <: AbstractMathProgModel
    model::AmplNLMathProgModel
end

struct CouenneNLLinearQuadraticModel <: AbstractLinearQuadraticModel
    model::AmplNLLinearQuadraticModel
end

SolverInterface.LinearQuadraticModel(s::CouenneNLSolver) = CouenneNLLinearQuadraticModel(LinearQuadraticModel(s.solver))

struct CouenneNLNonlinearModel <: AbstractNonlinearModel
    model::AmplNLNonlinearModel
end

SolverInterface.NonlinearModel(s::CouenneNLSolver) = CouenneNLNonlinearModel(NonlinearModel(s.solver))

@forward CouenneNLSolver.solver AmplNLWriter.getsolvername
@forward CouenneNLNonlinearModel.model SolverInterface.loadproblem!
@forward CouenneNLLinearQuadraticModel.model SolverInterface.loadproblem!
@forward CouenneNLMathProgModel.model AmplNLWriter.load_A!
@forward CouenneNLMathProgModel.model AmplNLWriter.loadcommon!

@forward CouenneNLMathProgModel.model SolverInterface.getvartype
@forward CouenneNLMathProgModel.model SolverInterface.setvartype!
@forward CouenneNLMathProgModel.model SolverInterface.getsense
@forward CouenneNLMathProgModel.model SolverInterface.setsense!
@forward CouenneNLMathProgModel.model SolverInterface.setwarmstart!
@forward CouenneNLMathProgModel.model SolverInterface.optimize!
@forward CouenneNLMathProgModel.model SolverInterface.status
@forward CouenneNLMathProgModel.model SolverInterface.getsolution
@forward CouenneNLMathProgModel.model SolverInterface.getobjval
@forward CouenneNLMathProgModel.model SolverInterface.numvar
@forward CouenneNLMathProgModel.model SolverInterface.numconstr
@forward CouenneNLMathProgModel.model SolverInterface.getsolvetime
@forward CouenneNLMathProgModel.model AmplNLWriter.get_solve_result
@forward CouenneNLMathProgModel.model AmplNLWriter.get_solve_result_num
@forward CouenneNLMathProgModel.model AmplNLWriter.get_solve_message
@forward CouenneNLMathProgModel.model AmplNLWriter.get_solve_exitcode

@forward CouenneNLLinearQuadraticModel.model SolverInterface.getvartype
@forward CouenneNLLinearQuadraticModel.model SolverInterface.setvartype!
@forward CouenneNLLinearQuadraticModel.model SolverInterface.getsense
@forward CouenneNLLinearQuadraticModel.model SolverInterface.setsense!
@forward CouenneNLLinearQuadraticModel.model SolverInterface.setwarmstart!
@forward CouenneNLLinearQuadraticModel.model SolverInterface.optimize!
@forward CouenneNLLinearQuadraticModel.model SolverInterface.status
@forward CouenneNLLinearQuadraticModel.model SolverInterface.getsolution
@forward CouenneNLLinearQuadraticModel.model SolverInterface.getobjval
@forward CouenneNLLinearQuadraticModel.model SolverInterface.numvar
@forward CouenneNLLinearQuadraticModel.model SolverInterface.numconstr
@forward CouenneNLLinearQuadraticModel.model SolverInterface.getsolvetime
@forward CouenneNLLinearQuadraticModel.model AmplNLWriter.get_solve_result
@forward CouenneNLLinearQuadraticModel.model AmplNLWriter.get_solve_result_num
@forward CouenneNLLinearQuadraticModel.model AmplNLWriter.get_solve_message
@forward CouenneNLLinearQuadraticModel.model AmplNLWriter.get_solve_exitcode

@forward CouenneNLNonlinearModel.model SolverInterface.getvartype
@forward CouenneNLNonlinearModel.model SolverInterface.setvartype!
@forward CouenneNLNonlinearModel.model SolverInterface.getsense
@forward CouenneNLNonlinearModel.model SolverInterface.setsense!
@forward CouenneNLNonlinearModel.model SolverInterface.setwarmstart!
@forward CouenneNLNonlinearModel.model SolverInterface.optimize!
@forward CouenneNLNonlinearModel.model SolverInterface.status
@forward CouenneNLNonlinearModel.model SolverInterface.getsolution
@forward CouenneNLNonlinearModel.model SolverInterface.getobjval
@forward CouenneNLNonlinearModel.model SolverInterface.numvar
@forward CouenneNLNonlinearModel.model SolverInterface.numconstr
@forward CouenneNLNonlinearModel.model SolverInterface.getsolvetime
@forward CouenneNLNonlinearModel.model AmplNLWriter.get_solve_result
@forward CouenneNLNonlinearModel.model AmplNLWriter.get_solve_result_num
@forward CouenneNLNonlinearModel.model AmplNLWriter.get_solve_message
@forward CouenneNLNonlinearModel.model AmplNLWriter.get_solve_exitcode

end # module


text:
###############################################
#FEniCS tutorial demo program: Deflection of a membrane.
#  -Laplace(w) = p  in the unit circle
#            w = 0  on the boundary
#The load p is a Gaussian function centered at (0, 0.6).
###############################################

module ft02

using FenicsPy

# Create mesh and define function space
Ω = Circle(Point(0, 0), 1)
mesh = generate_mesh(Ω, 64)
V = FunctionSpace(mesh, "P", 2)

# Define boundary condition
w_D = Constant(0)

bc = DirichletBC(V, w_D, "on_boundary")

# Define load
β = 8
R0 = 0.6
p = Expression("4*exp(-pow(beta, 2)*(pow(x[0], 2) + pow(x[1] - R0, 2)))",
               degree=1, beta=β, R0=R0)

# Define variational problem
w = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(w), grad(v))*dx
L = p*v*dx

# Compute solution
w = FeFunction(V)
solve(a == L, w, bc)

p = interpolate(p, V)

# Save solution to file in VTK format
vtkfile_w = File("poisson_membrane/deflection.pvd")
vtkfile_w << w
vtkfile_p = File("poisson_membrane/load.pvd")
vtkfile_p << p

# Plot solution
#plot(w, title="Deflection")
#plot(p, title="Load")

# Curve plot along x = 0 comparing p and w
#import numpy as np
#import matplotlib.pyplot as plt
#tol = 0.001  # avoid hitting points outside the domain
#y = np.linspace(-1 + tol, 1 - tol, 101)
#points = [(0, y_) for y_ in y]  # 2D points
#w_line = np.array([w(point) for point in points])
#p_line = np.array([p(point) for point in points])
#plt.plot(y, 50*w_line, 'k', linewidth=2)  # magnify w
#plt.plot(y, p_line, 'b--', linewidth=2)
#plt.grid(True)
#plt.xlabel('$y$')
#plt.legend(['Deflection ($\\times 50$)', 'Load'], loc='upper left')
#plt.savefig('poisson_membrane/curves.pdf')
#plt.savefig('poisson_membrane/curves.png')

end # module ft02


text:
module Lighting
  using CG;
  export PointLight, Light, SceneLights;

  abstract Light;

  type PointLight <: Light
    position::Vec4f
  end

  type SceneLights
    lights::Vector{Light}
  end
end

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


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

using Test
import Plots
using PokerBots
using TexasHoldem
const PB = PokerBots
const TH = TexasHoldem

function plot_stats(stats, n_games)

    player_bank_rolls = stats[:player_bank_rolls]
    p = Plots.plot()
    for (player_bank_roll, pid, lf) in player_bank_rolls
        bot_name = nameof(typeof(lf))
        Plots.plot!(1:n_games, player_bank_roll;
            markershape=:auto,
            label="$bot_name prob=$(lf.prob)",
            xlabel="n games",
            ylabel="\$",
            legend = :topleft,
            dpi = 300
        )
    end
    # Plots.savefig("Stats.png")
end

@testset "Battle" begin
    n_tournaments = 50
    bots = (
        PB.BotRandom(0.4),
        PB.BotRandom(0.4),
        PB.BotRandom(0.6),
        PB.BotRandom(0.6),
    )
    @time stats, n_games = PB.battle!(n_tournaments, bots...)
    plot_stats(stats, n_games)
end


text:
# simulate MUSIC with a circular array
using DirectionFinding

begin
    fc = 1e9
    λ = 299.79e6/fc
    a1 = 0.95; a2 = 6;
    N₀ = 0.1
    ca = circulararray(11, 2λ, σ²=0.0001)
    s1 = SignalSource(a1, halfsine_train(1e6))
    s2 = SignalSource(a2, halfsine_train(1e6))
    sim = Simulation(ca, fc, N₀, s1, s2)
    mps = MUSIC(sim, 2.00001e6, 1000)
    ϕ = 0:0.01:2π-0.01
    peaks = findpeaks(mps, ϕ, n=2)
    println("Angles: $a1 -- $a2")
    println("Found: $(peaks[1]) -- $(peaks[2])")
    plot(ϕ, mps.(ϕ))
end


text:
module KernelsTests

using Test
using Gridap
using TensorValues
using Gridap.Kernels: CellSumKernel
using Gridap.Kernels: LinCombKernel
using Gridap.Kernels: VarinnerKernel
using Gridap.Kernels: PhysGradKernel
using Gridap.Kernels: IntegrateNumberKernel
using Gridap.Kernels: IntegrateArrayKernel
using Gridap.Kernels: IntegrateKernel
using Gridap.Kernels: CellNewAxisKernel

# NumberKernelFromFunction

k = NumberKernelFromFunction(-)

test_number_kernel(k,1,4,3)

test_number_kernel(k,-1.0,1.0)

k = NumberKernelFromFunction(sum)

test_number_kernel(k,7,[1,2,4])

# ArrayKernelFromBroadcastedFunction

k = ArrayKernelFromBroadcastedFunction(*)

test_array_kernel(k,[3,4,3],[1,2,3],[3,2,1])

test_array_kernel(k,[6,4,2],2,[3,2,1])

test_array_kernel(k,[6,4,2],[3,2,1],2)

test_array_kernel(k,[6,4,2],[3,2,1],2,1)

u = rand(Int,2,3,1)
v = rand(Int,1,3,4)
w = u .* v
test_array_kernel(k,w,u,v)
test_array_kernel(k,w,u,v,1)

k = ArrayKernelFromBroadcastedFunction(+)

v1 = VectorValue(2,3)
v2 = VectorValue(3,2)
v3 = VectorValue(1,2)
u = [v1,v2,v3]
v = v1
w = u .+ v
test_array_kernel(k,w,u,v)

w = broadcast(+,u,v,0)
test_array_kernel(k,w,u,v,0)

# CellSumKernel

D = 1
k = CellSumKernel{D}()
a = rand(2,3,4)
b = reshape(sum(a,dims=D),(3,4))
test_array_kernel(k,b,a)

D = 2
k = CellSumKernel{D}()
a = rand(2,3,4)
b = reshape(sum(a,dims=D),(2,4))
test_array_kernel(k,b,a)

D = 3
k = CellSumKernel{D}()
a = rand(2,3,4)
b = reshape(sum(a,dims=D),(2,3))
test_array_kernel(k,b,a)

# LinCombKernel

k = LinCombKernel()
a = rand(4,8)
b = rand(4)
r = reshape(sum(outer.(a,reshape(b,(4,1))),dims=1),(8,))
test_array_kernel(k,r,a,b)

k = LinCombKernel()
ai = VectorValue(2.0,3.1,4.3)
a = fill(ai,4,8)
b = rand(4)
r = reshape(sum(outer.(a,reshape(b,(4,1))),dims=1),(8,))
test_array_kernel(k,r,a,b)

k = LinCombKernel()
ai = VectorValue(2.0,3.1,4.3)
a = fill(ai,4,8)
b = fill(ai,4)
r = reshape(sum(outer.(a,reshape(b,(4,1))),dims=1),(8,))
test_array_kernel(k,r,a,b)

# VarinnerKernel

k = VarinnerKernel()
a = rand(4,8)
b = rand(8)
r = inner.(a,reshape(b,(1,8)))
test_array_kernel(k,r,a,b)

k = VarinnerKernel()
a = rand(4,8)
b = rand(5,8)
r = inner.(reshape(a,(4,1,8)),reshape(b,(1,5,8)))
test_array_kernel(k,r,a,b)

ai = VectorValue(2.0,3.1,4.3)
bi = VectorValue(2.1,1.1,4.1)

a = fill(ai,4,8)
b = fill(bi,8)
r = inner.(a,reshape(b,(1,8)))
test_array_kernel(k,r,a,b)

a = fill(ai,4,8)
b = fill(bi,5,8)
r = inner.(reshape(a,(4,1,8)),reshape(b,(1,5,8)))
test_array_kernel(k,r,a,b)

# PhysGradKernel

k = PhysGradKernel()

ai = TensorValue(1.0,2.0,0.0,2.0)
bi = 3.0
a = fill(ai,8)
b = fill(bi,5,8)
r = outer.(reshape(inv.(a),1,8),b)
test_array_kernel(k,r,a,b)

bi = VectorValue(3.4,2.5)
b = fill(bi,5,8)
r = reshape(inv.(a),1,8) .* b
test_array_kernel(k,r,a,b)

# IntegrateNumberKernel

k = IntegrateNumberKernel{Int}()

f = [1,3,4,3]
j = [1,2,4,1]
w = [2,4,5,9]
r = sum( f .* j .* w )
test_number_kernel(k,r,f,j,w)

# IntegrateArrayKernel

k = IntegrateArrayKernel()
f = rand(3,2,4)
j = rand(4)
w = rand(4)

r =  reshape( sum(f .* reshape(j,(1,1,4)) .* reshape(w,(1,1,4)), dims=3), (3,2))

test_array_kernel(k,r,f,j,w)

k = IntegrateKernel(1,Float64)
@test isa(k,IntegrateNumberKernel)

k = IntegrateKernel(2,Float64)
@test isa(k,IntegrateArrayKernel)

# CellNewAxisKernel

a = rand(Int,3,2)

k = CellNewAxisKernel{1}()
b = collect(reshape(a,(1,3,2)))
test_array_kernel(k,b,a)

k = CellNewAxisKernel{2}()
b = collect(reshape(a,(3,1,2)))
test_array_kernel(k,b,a)

k = CellNewAxisKernel{3}()
b = collect(reshape(a,(3,2,1)))
test_array_kernel(k,b,a)

end # module


text:
# StaticMaps/Masks.jl
#
# MIT License
#
# Copyright (c) 2019 Brandon Gomes
#
# 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.

"""
    Masks

...
"""
baremodule Masks

import Base: |, &

include("base.jl")

export StaticMask, apply, MaskNone, MaskAll, OR, AND

"""
    StaticMask

...
"""
abstract type StaticMask{F, Args} <: BroadcastedMap{F, Args} end

"""
    apply(::Type{<: StaticMask}, ::AbstractArray)

...
"""
apply(M::Type{<: StaticMask}, array::AbstractArray) = array[M(array)]

"""
    MaskNone

...
"""
abstract type MaskNone <: StaticMask{x -> true, ()} end

"""
    MaskAll

...
"""
abstract type MaskAll <: StaticMask{x -> false, ()} end

"""
    OR{A, B}

...
"""
abstract type OR{A <: StaticMask, B <: StaticMask} <: AbstractStaticMap end

"""
    AND{A, B}

...
"""
abstract type AND{A <: StaticMask, B <: StaticMask} <: AbstractStaticMap end

"""
    A | B

...
"""
Base.:|(A::Type{AA}, B::Type{BB}) where {AA <: StaticMask, BB <: StaticMask} = OR{A, B}

"""
    A & B

...
"""
Base.:&(A::Type{AA}, B::Type{BB}) where {AA <: StaticMask, BB <: StaticMask} = AND{A, B}

"""
    MaskNone | B === MaskNone

...
"""
Base.:|(A::Type{AA}, B::Type{BB}) where {AA <: MaskNone, BB <: StaticMask} = A

"""
    A | MaskNone === MaskNone

...
"""
Base.:|(A::Type{AA}, B::Type{BB}) where {AA <: StaticMask, BB <: MaskNone} = B

"""
    MaskAll | B === B

...
"""
Base.:|(A::Type{AA}, B::Type{BB}) where {AA <: MaskAll, BB <: StaticMask} = B

"""
    A | MaskAll === A

...
"""
Base.:|(A::Type{AA}, B::Type{BB}) where {AA <: StaticMask, BB <: MaskAll} = A

"""
    MaskNone & B === B

...
"""
Base.:&(A::Type{AA}, B::Type{BB}) where {AA <: MaskNone, BB <: StaticMask} = B

"""
    A & MaskNone === A

...
"""
Base.:&(A::Type{AA}, B::Type{BB}) where {AA <: StaticMask, BB <: MaskNone} = A

"""
    MaskAll & B === MaskAll

...
"""
Base.:&(A::Type{AA}, B::Type{BB}) where {AA <: MaskAll, BB <: StaticMask} = A

"""
    A & MaskAll === MaskAll

...
"""
Base.:&(A::Type{AA}, B::Type{BB}) where {AA <: StaticMask, BB <: MaskAll} = B

"""
    MaskNone | MaskAll === MaskNone

...
"""
Base.:|(A::Type{AA}, B::Type{BB}) where {AA <: MaskNone, BB <: MaskAll} = A

"""
    MaskAll | MaskNone === MaskNone

...
"""
Base.:|(A::Type{AA}, B::Type{BB}) where {AA <: MaskAll, BB <: MaskNone} = B

"""
    MaskNone & MaskAll === MaskAll

...
"""
Base.:&(A::Type{AA}, B::Type{BB}) where {AA <: MaskNone, BB <: MaskAll} = B

"""
    MaskAll & MaskNone === MaskAll

...
"""
Base.:&(A::Type{AA}, B::Type{BB}) where {AA <: MaskAll, BB <: MaskNone} = A

"""
    OR{A, B}(::AbstractArray)

...
"""
(::Type{<:OR{A, B}})(array::AbstractArray) where {A <: StaticMask, B <: StaticMask} = broadcast(|, A(array), B(array))

"""
    AND{A, B}(::AbstractArray)

...
"""
(::Type{<:AND{A, B}})(array::AbstractArray) where {A <: StaticMask, B <: StaticMask} = broadcast(&, A(array), B(array))

"""
    ComparisonMask{CMP, V}

...
"""
abstract type ComparisonMask{CMP, V} <: StaticMask{CMP, (V)} end

"""
    EQ{V}

...
"""
const EQ = ComparisonMask{==}

"""
    NEQ{V}

...
"""
const NEQ = ComparisonMask{!=}

"""
    LT{V}

...
"""
const LT = ComparisonMask{<}

"""
    LEQ

...
"""
const LEQ = ComparisonMask{<=}

"""
    GT{V}

...
"""
const GT = ComparisonMask{>}

"""
    GEQ{V}

...
"""
const GEQ = ComparisonMask{>=}

end # baremodule Masks


text:
@testset "6.5.1 (c+d x)^m (a+b sech)^n" begin
    (a, b, c, d, x, ) = @variables a b c d x

    #= ::Package:: =#

    #= ::Section:: =#
    #=Integrands*of*the*form*(c+d*x)^m*sech(a+b*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*sech(a+b*x)^n=#

    @test_int [(c + d*x)^3*sech(a + b*x), x, 9, (2*(c + d*x)^3*atan(ℯ^(a + b*x)))/b - (3*im*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(a + b*x)))/b^2 + (3*im*d*(c + d*x)^2*PolyLog(2, im*ℯ^(a + b*x)))/b^2 + (6*im*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(a + b*x)))/b^3 - (6*im*d^2*(c + d*x)*PolyLog(3, im*ℯ^(a + b*x)))/b^3 - (6*im*d^3*PolyLog(4, (-im)*ℯ^(a + b*x)))/b^4 + (6*im*d^3*PolyLog(4, im*ℯ^(a + b*x)))/b^4]
    @test_int [(c + d*x)^2*sech(a + b*x), x, 7, (2*(c + d*x)^2*atan(ℯ^(a + b*x)))/b - (2*im*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(a + b*x)))/b^2 + (2*im*d*(c + d*x)*PolyLog(2, im*ℯ^(a + b*x)))/b^2 + (2*im*d^2*PolyLog(3, (-im)*ℯ^(a + b*x)))/b^3 - (2*im*d^2*PolyLog(3, im*ℯ^(a + b*x)))/b^3]
    @test_int [(c + d*x)^1*sech(a + b*x), x, 5, (2*(c + d*x)*atan(ℯ^(a + b*x)))/b - (im*d*PolyLog(2, (-im)*ℯ^(a + b*x)))/b^2 + (im*d*PolyLog(2, im*ℯ^(a + b*x)))/b^2]
    @test_int [sech(a + b*x)/(c + d*x)^1, x, 0, Unintegrable(sech(a + b*x)/(c + d*x), x)]

    @test_int [(c + d*x)^3*sech(a + b*x)^2, x, 6, (c + d*x)^3/b - (3*d*(c + d*x)^2*log(1 + ℯ^(2*(a + b*x))))/b^2 - (3*d^2*(c + d*x)*PolyLog(2, -ℯ^(2*(a + b*x))))/b^3 + (3*d^3*PolyLog(3, -ℯ^(2*(a + b*x))))/(2*b^4) + ((c + d*x)^3*tanh(a + b*x))/b]
    @test_int [(c + d*x)^2*sech(a + b*x)^2, x, 5, (c + d*x)^2/b - (2*d*(c + d*x)*log(1 + ℯ^(2*(a + b*x))))/b^2 - (d^2*PolyLog(2, -ℯ^(2*(a + b*x))))/b^3 + ((c + d*x)^2*tanh(a + b*x))/b]
    @test_int [(c + d*x)^1*sech(a + b*x)^2, x, 2, -((d*log(cosh(a + b*x)))/b^2) + ((c + d*x)*tanh(a + b*x))/b]
    @test_int [sech(a + b*x)^2/(c + d*x)^1, x, 0, Unintegrable(sech(a + b*x)^2/(c + d*x), x)]

    @test_int [(c + d*x)^3*sech(a + b*x)^3, x, 15, -((6*d^2*(c + d*x)*atan(ℯ^(a + b*x)))/b^3) + ((c + d*x)^3*atan(ℯ^(a + b*x)))/b + (3*im*d^3*PolyLog(2, (-im)*ℯ^(a + b*x)))/b^4 - (3*im*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(a + b*x)))/(2*b^2) - (3*im*d^3*PolyLog(2, im*ℯ^(a + b*x)))/b^4 + (3*im*d*(c + d*x)^2*PolyLog(2, im*ℯ^(a + b*x)))/(2*b^2) + (3*im*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(a + b*x)))/b^3 - (3*im*d^2*(c + d*x)*PolyLog(3, im*ℯ^(a + b*x)))/b^3 - (3*im*d^3*PolyLog(4, (-im)*ℯ^(a + b*x)))/b^4 + (3*im*d^3*PolyLog(4, im*ℯ^(a + b*x)))/b^4 + (3*d*(c + d*x)^2*sech(a + b*x))/(2*b^2) + ((c + d*x)^3*sech(a + b*x)*tanh(a + b*x))/(2*b)]
    @test_int [(c + d*x)^2*sech(a + b*x)^3, x, 9, ((c + d*x)^2*atan(ℯ^(a + b*x)))/b - (d^2*atan(sinh(a + b*x)))/b^3 - (im*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(a + b*x)))/b^2 + (im*d*(c + d*x)*PolyLog(2, im*ℯ^(a + b*x)))/b^2 + (im*d^2*PolyLog(3, (-im)*ℯ^(a + b*x)))/b^3 - (im*d^2*PolyLog(3, im*ℯ^(a + b*x)))/b^3 + (d*(c + d*x)*sech(a + b*x))/b^2 + ((c + d*x)^2*sech(a + b*x)*tanh(a + b*x))/(2*b)]
    @test_int [(c + d*x)^1*sech(a + b*x)^3, x, 6, ((c + d*x)*atan(ℯ^(a + b*x)))/b - (im*d*PolyLog(2, (-im)*ℯ^(a + b*x)))/(2*b^2) + (im*d*PolyLog(2, im*ℯ^(a + b*x)))/(2*b^2) + (d*sech(a + b*x))/(2*b^2) + ((c + d*x)*sech(a + b*x)*tanh(a + b*x))/(2*b)]
    @test_int [sech(a + b*x)^3/(c + d*x)^1, x, 0, Unintegrable(sech(a + b*x)^3/(c + d*x), x)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(c+d*x)^(m/2)*sech(a+b*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*sech(a+b*x)^(n/2)=#

    @test_int [x/sech(x)^(3/2) - (1/3)*x*sqrt(sech(x)), x, 4, -(4/(9*sech(x)^(3/2))) + (2*x*sinh(x))/(3*sqrt(sech(x)))]
    @test_int [x/sech(x)^(5/2) - (3/5)*x/sqrt(sech(x)), x, 4, -(4/(25*sech(x)^(5/2))) + (2*x*sinh(x))/(5*sech(x)^(3/2))]
    @test_int [x/sech(x)^(7/2) - (5/21)*x*sqrt(sech(x)), x, 5, -(4/(49*sech(x)^(7/2))) - 20/(63*sech(x)^(3/2)) + (2*x*sinh(x))/(7*sech(x)^(5/2)) + (10*x*sinh(x))/(21*sqrt(sech(x)))]
    @test_int [x^2/sech(x)^(3/2) - (1/3)*x^2*sqrt(sech(x)), x, 7, -((8*x)/(9*sech(x)^(3/2))) - (16/27)*im*sqrt(cosh(x))*Elliptic.F((im*x)/2, 2)*sqrt(sech(x)) + (16*sinh(x))/(27*sqrt(sech(x))) + (2*x^2*sinh(x))/(3*sqrt(sech(x)))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(c+d*x)^(m/2)*sech(a+b*x)^(n/2)=#
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.5 Hyperbolic secant/6.5.1 (c+d x)^m (a+b sech)^n.jl


text:
function obtain_const_or_stmt(@nospecialize(x), ci::CodeInfo)
    if x isa SSAValue
        stmt = ci.code[x.id]
        typ = ci.ssavaluetypes[x.id]
        if typ isa Const
            return typ.val, typeof(typ.val)
        else
            return stmt, widenconst(typ)
        end
    elseif x isa QuoteNode
        return x.value, typeof(x.value)
    elseif x isa Const
        return x.val, typeof(x.val)
    elseif x isa GlobalRef
        if isdefined(x.mod, x.name) && isconst(x.mod, x.name)
            val = getfield(x.mod, x.name)
        else
            # TODO: move this to parsing time
            throw(UndefVarError(x.name))
        end

        return val, typeof(val)
    else
        # special value
        return x, typeof(x)
    end
end

obtain_const(@nospecialize(x), ci::CodeInfo) = obtain_const_or_stmt(x, ci)[1]

function obtain_qasm_ifnot_cond(cond::SSAValue, ci::CodeInfo)
    cond_stmt = ci.code[cond.id]

    (cond_stmt isa Expr && cond_stmt.head === :call) || error("invalid cond statement: $cond_stmt")
    x = obtain_const(cond_stmt.args[2], ci)
    if x isa Int
        creg = cond_stmt.args[3]
    else
        x = obtain_const(cond_stmt.args[3], ci)
        if x isa Int
            creg = cond_stmt.args[2]
        else
            error("one of the value in if condition must be constant for QASM compatible program")
        end
    end

    creg isa SlotNumber || creg isa SSAValue || error("invalid statement: $creg")
    return creg, x
end

function obtain_const_measure_stmt(stmt, ci::CodeInfo)
    if stmt.head === :(=)
        cvar, body = _extract_measure(stmt)
    else
        cvar = nothing
        body = stmt
    end

    if body.head === :invoke
        locs = obtain_const(body.args[3], ci)::Locations
    elseif body.head === :call
        locs = obtain_const(body.args[2], ci)::Locations
    end
    return cvar, locs
end

function obtain_const_gate_stmt(stmt, ci::CodeInfo)
    if stmt.head === :invoke
        gate, gt = obtain_const_or_stmt(stmt.args[3], ci)
        locs = obtain_const(stmt.args[4], ci)::Locations
    elseif stmt.head === :call
        gate, gt = obtain_const_or_stmt(stmt.args[2], ci)
        locs = obtain_const(stmt.args[3], ci)::Locations
    end
    return gate, gt, locs
end

function obtain_const_ctrl_stmt(stmt, ci::CodeInfo)
    # NOTE: QASM compatible code should have constant location
    # QASM compatible code should only contain control X gate
    # as control gates, for other control gates, it should be
    # either decomposed or error, so we will assume the gate
    # is a constant here.
    if stmt.head === :invoke
        gate, gt = obtain_const_or_stmt(stmt.args[3], ci)
        locs = obtain_const(stmt.args[4], ci)::Locations
        ctrl = obtain_const(stmt.args[5], ci)::CtrlLocations
    elseif stmt.head === :call
        gate, gt = obtain_const_or_stmt(stmt.args[2], ci)
        locs = obtain_const(stmt.args[3], ci)::Locations
        ctrl = obtain_const(stmt.args[4], ci)::CtrlLocations
    end

    return gate, gt, locs, ctrl
end

function move_to_first_quantum_stmt(ci::CodeInfo, st)
    pc′ = st.pc + 1
    local stmt′
    while pc′ <= length(ci.code)
        stmt′ = ci.code[pc′]
        if stmt′ isa Expr
            is_quantum_statement(stmt′) && break
        else
            # NOTE:
            # we can allow some constant statements
            # but not control flows or other nodes
            error("unexpected statement inside control flow: $stmt′")
        end
        pc′ += 1
    end

    st.pc = pc′
    st.stmt = stmt′
    return stmt′
end


text:
using Documenter, NCDatasets

makedocs(modules = [NCDatasets], sitename = "NCDatasets.jl")

deploydocs(
    repo = "github.com/Alexander-Barth/NCDatasets.jl.git",
)


text:
include("q1/p1.jl");
include("q1/p2.jl");
include("q1/p3.jl");
include("q1/p4.jl");
include("q1/p5.jl");
include("q1/p6.jl");
# include("q1/p7.jl");


text:
export simplify

function simplify(t::Term, rules, maxstep=1000)
    for _ in 1:maxstep
        prev = t
        for r in rules
            t = replace(t, r)
        end
        prev == t && return t
    end
    # fail to simplify
    return t
end


text:
relation=@equ Log(:a,:a)=1
@test Log(:e)&relation==Log(9,9)&relation==Log(:x+:y,:x+:y)&relation==1


text:
###
### Indexing
###

@inline function inbounds_getindex2(x::Kmer{A,K,N}, i::Integer) where {A,K,N}
    # Emulation of BitIndex type
    i′ = i + div(64N - 2K, 2)
    val = (i′ - 1) << 1
    idx = (val >> 6) + 1
    off = 62 - (val & (UInt8(64) - 0x01))
    @inbounds begin
        chunk = x.data[idx]
    end
    bits = (chunk >> off) & UInt64(3)
    return reinterpret(eltype(x), 0x01 << bits)
end

"""
Get the Vector or NTuple of unsigned integers that the encoded
"""
function packed_data end

@inline packed_data(seq::Kmer) = seq.data

@inline function inbounds_getindex(seq::Kmer, i::Integer)
    bidx = bitindex(seq, i)
    elem = extract_encoded_element(bidx, packed_data(seq))
    return decode(Alphabet(seq), elem)
end


text:
struct SampleSize{T<:Real} <: Distributions.DiscreteUnivariateDistribution

  design::Design
  p::T

  function SampleSize{T}(design::Design, p::T) where {T<:Real}

    @checkprob p
    new(design, p)

  end # inner constructor

end # SampleSize

SampleSize(design::Design, p::T) where {T<:Real} = SampleSize{T}(design, p)

design(ss::SampleSize) = ss.design

Base.show(io::IO, d::SampleSize) = print("SampleSize")

Base.size(::SampleSize) = ()

Base.getindex(ss::SampleSize, i) = ss

rand(d::SampleSize) = d.design |> interimsamplesize |> n1 -> Distributions.Binomial(n1, d.p) |> rand |> x1 -> samplesize(d.design, x1)

function pdf(d::SampleSize, n::Int)

  res = 0.0
  for x1 in 0:interimsamplesize(d.design)
    if samplesize(d.design, x1) == n
      res += Distributions.pdf(Distributions.Binomial(interimsamplesize(d.design), d.p), x1)
    end
  end
  return res

end # pdf

minimum(d::SampleSize) = interimsamplesize(d.design)

maximum(d::SampleSize) = maximum(samplesize(d.design))

function quantile(d::SampleSize, q::Real)

    nrange = minimum(d):maximum(d)
    res    = 0
    while Distributions.cdf(d, res) <= q
        res += 1
    end
    return res

end # quantile

function mean(d::SampleSize)

    res = 0.0
    for n in minimum(d):maximum(d)
        res += pdf(d, n)*n
    end
    return(res)

end # mean

function var(d::SampleSize)

    res = 0.0
    m   = mean(d)
    for n in minimum(d):maximum(d)
        res += pdf(d, n)*(n - m)^2
    end
    return(res)

end # var

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
using DataStructures
using Dates
using DecFP
using JSON3
using LightBSON
using Sockets
using StructTypes
using Test
using UUIDs

struct EmptyStruct end

@testset "LightBSON.jl" begin
    include("reader_tests.jl")
    include("index_tests.jl")
    include("indexed_reader_tests.jl")
    include("write_buffer_tests.jl")
    include("writer_tests.jl")
    include("struct_tests.jl")
    include("corpus_tests.jl")
    include("object_id_tests.jl")
    include("convenience_tests.jl")
    include("representation_tests.jl")
end


text:
function benchmark_downconvert()
    results_min = DataFrame(
        Samples = SAMPLES,
        sCPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_16ant = zeros(Float64,length(SAMPLES)),
        GPU_time_1ant = zeros(Float64,length(SAMPLES)),
        GPU_time_4ant = zeros(Float64,length(SAMPLES)),
        GPU_time_16ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_16ant = zeros(Float64,length(SAMPLES))
    )
    results_med = DataFrame(
        Samples = SAMPLES,
        sCPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_16ant = zeros(Float64,length(SAMPLES)),
        GPU_time_1ant = zeros(Float64,length(SAMPLES)),
        GPU_time_4ant = zeros(Float64,length(SAMPLES)),
        GPU_time_16ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_16ant = zeros(Float64,length(SAMPLES))
    )
    results_mean = DataFrame(
        Samples = SAMPLES,
        sCPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_16ant = zeros(Float64,length(SAMPLES)),
        GPU_time_1ant = zeros(Float64,length(SAMPLES)),
        GPU_time_4ant = zeros(Float64,length(SAMPLES)),
        GPU_time_16ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_16ant = zeros(Float64,length(SAMPLES))
    )
    rowpos = Int32(1)
    for N in SAMPLES
        for M in ANTENNA
            # init csv column write position
            M==1 ? columnpos=2 : M==4 ? columnpos=3 : M==16 ? columnpos=4 : throw("Error in antenna num")
            # init signals
            phases = 2π * (1:N) * 1000 / 2.5e6
            cpusignal = Vector{ComplexF32}(cos.(phases) + 1im.*sin.(phases))
            cpucarrier = cpusignal
            gpusignal = CuVector{ComplexF32}(cos.(phases) + 1im.*sin.(phases))
            gpucarrier = gpusignal
            cpudwnsignal = zeros(ComplexF32, N)
            gpudwnsignal = CUDA.zeros(ComplexF32, N)

            # init signals in StructArrays
            scpusignal = StructArray{ComplexF32}(
                (real(cpusignal),
                imag(cpusignal)))
            scpucarrier = copy(scpusignal)
            sgpusignal = StructArray{ComplexF32}(
                (real(gpusignal),
                imag(gpusignal)))
            sgpucarrier = copy(sgpusignal)
            scpudwnsignal = StructArray{ComplexF32}(
                    (real(zeros(Float32, N)),imag(zeros(Float32, N)))
                    )
            sgpudwnsignal = StructArray{ComplexF32}(
                    (real(CUDA.zeros(Float32, N)), imag(CUDA.zeros(Float32, N))))

            # init signals as matrices
            if M > 1
                cpusignal = cpusignal .* ones(ComplexF32, N, M)
                gpusignal = gpusignal .* CUDA.ones(ComplexF32, N, M)
                cpudwnsignal = zeros(ComplexF32, N, M)
                gpudwnsignal = CUDA.zeros(ComplexF32, N, M)

                scpusignal = StructArray{ComplexF32}(
                    (real(cpusignal),
                    imag(cpusignal)))
                sgpusignal = StructArray{ComplexF32}(
                    (real(gpusignal),
                    imag(gpusignal)))
                scpudwnsignal = StructArray{ComplexF32}(
                    (real(zeros(Float32,N,M)),
                    imag(zeros(Float32,N,M))))
                sgpudwnsignal = StructArray{ComplexF32}(
                    (real(CUDA.zeros(Float32,N,M)),
                    imag(CUDA.zeros(Float32,N,M))))
            end

            # do the benchmark
            println("Benchmarking downconvert on CPU: StructArray{ComplexF32} ", N," samples, ", M, " antenna...")
            result = @benchmark cpu_downconvert!(
                $scpudwnsignal.re,
                $scpudwnsignal.im,
                $scpucarrier.re,
                $scpucarrier.im,
                $scpusignal.re,
                $scpusignal.im,
                1,
                $N
            )
            println(minimum(result).time)
            results_min[rowpos, columnpos] = minimum(result).time
            results_med[rowpos, columnpos] = median(result).time
            results_mean[rowpos, columnpos] = mean(result).time
            println("Benchmarking downconvert on GPU: CuArray{ComplexF32} ", N, " samples, ", M, " antenna...")
            result = @benchmark CUDA.@sync gpu_downconvert!(
                $gpudwnsignal,
                $gpucarrier,
                $gpusignal,
                1,
                $N
            )
            println(minimum(result).time)
            results_min[rowpos, columnpos+3] = minimum(result).time
            results_med[rowpos, columnpos+3] = median(result).time
            results_mean[rowpos, columnpos+3] = mean(result).time
            println("Benchmarking downconvert on GPU: StructArray{ComplexF32}(CuArray,CuArray) ", N, " samples, ", M, " antenna...")
            result = @benchmark CUDA.@sync gpu_downconvert!(
                $sgpudwnsignal.re,
                $sgpudwnsignal.im,
                $sgpucarrier.re,
                $sgpucarrier.im,
                $sgpusignal.re,
                $sgpusignal.im,
                1,
                $N
            )
            println(minimum(result).time)
            results_min[rowpos, columnpos+6] = minimum(result).time
            results_med[rowpos, columnpos+6] = median(result).time
            results_mean[rowpos, columnpos+6] = mean(result).time
        end
        rowpos += 1
    end
    CSV.write("data/downconvert_min.csv", results_min)
    CSV.write("data/downconvert_med.csv", results_med)
    CSV.write("data/downconvert_mean.csv", results_mean)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export notificationsprovider

"""
    notificationsprovider(;kwargs...)
    notificationsprovider(children::Any;kwargs...)
    notificationsprovider(children_maker::Function;kwargs...)

A NotificationsProvider component.
Mantine notifications system. For more information, see: https://mantine.dev/others/notifications/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `autoClose` (Real | a value equal to: false; optional): Auto close timeout for all notifications, false to disable auto close, can be overwritten for individual notifications by showNotification function
- `containerWidth` (Real; optional): Notification width in px, cannot exceed 100%
- `limit` (Real; optional): Maximum amount of notifications displayed at a time, other new notifications will be added to queue
- `notificationMaxHeight` (Real; optional): Notification max-height in px, used for transitions
- `position` (a value equal to: "top-left", "top-right", "top-center", "bottom-left", "bottom-right", "bottom-center"; optional): Notifications position
- `zIndex` (Real; optional): Notifications container z-index
"""
function notificationsprovider(; kwargs...)
        available_props = Symbol[:children, :id, :autoClose, :containerWidth, :limit, :notificationMaxHeight, :position, :zIndex]
        wild_props = Symbol[]
        return Component("notificationsprovider", "NotificationsProvider", "dash_mantine_components", available_props, wild_props; kwargs...)
end

notificationsprovider(children::Any; kwargs...) = notificationsprovider(;kwargs..., children = children)
notificationsprovider(children_maker::Function; kwargs...) = notificationsprovider(children_maker(); kwargs...)


text:
struct Point2D{T}
    x::T
    y::T
end

Base.:-(P1::Point2D{T}, P2::Point2D{T}) where T = Point2D{T}(P1.x-P2.x, P1.y-P2.y)
Base.:+(P1::Point2D{T}, P2::Point2D{T}) where T = Point2D{T}(P1.x-P2.x, P1.y-P2.y)

struct Point4D{T}
    x::T
    y::T
    cur::T  # Curvature at this point, positive if dθ/ds > 0
    θ::T  # tangential direction in global coordinates
end

function Point4To2D(P::Point4D{T}) where T
    return Point2D{T}(P.x, P.y)
end

function Distance2D(Sp::Point2D{Float64}, Ep::Point2D{Float64})
    return sqrt((Sp.x-Ep.x)^2 + (Sp.y-Ep.y)^2)
end

function Distance2D(Sp::Array{Float64, 1}, Ep::Array{Float64, 1})
    return sqrt((Sp[1]-Ep[1])^2 + (Sp[2]-Ep[2])^2)
end

function Distance2D(Sp::Point4D{Float64}, Ep::Point4D{Float64})
    return sqrt((Sp.x-Ep.x)^2 + (Sp.y-Ep.y)^2)
end

function round(x::Float64, pre::Union{Float64, Int64}) # pre: precision
    n = floor(Int, x/pre)
    m = x - n*pre
    if m >= pre/2
        return trunc((n+1)*pre; digits=10)
    else
        return trunc(n*pre; digits=10)
    end
end

struct Geometry
    Sp::Point4D{Float64} # Start Point4D
    Ep::Point4D{Float64} # End Point4D
    Length::Float64
    Vstd::Float64 # the standard velocity of this Geometry
    GType::String
end

function Geometry(Sp::Point4D{Float64}, Ep::Point4D{Float64}, Vstd::Float64)
    if Sp.cur == Ep.cur
        if Sp.cur != 0.0
            Length = round(abs((Sp.θ - Ep.θ)/Sp.cur), 1.0e-14)
            GType = "arc"
        else
            return Geometry(Point2D{Float64}(Sp.x, Sp.y), Point2D{Float64}(Ep.x, Ep.y), Vstd) #GType = "line"
        end
    else
        GType = "clothoid"
        Length = 1.0
    end
    return Geometry(Sp, Ep, Length, Vstd, GType)
end

function Geometry(Sp::Point2D{Float64}, Ep::Point2D{Float64}, Vstd::Float64) # constructor for line element
    Length = round(Distance2D(Sp, Ep), 1.0e-14)
    Δx = Ep.x-Sp.x
    Δy = Ep.y-Sp.y
    Theta = 0.0
    if Δx != 0
        Theta = atan(Δy/Δx)
    elseif Δy > 0
        Theta = π/2
    elseif Δy < 0
        Theta = -π/2
    else
        error("Start point and end point are the same")
    end
    return Geometry(Point4D{Float64}(Sp.x, Sp.y, 0.0, Theta), Point4D{Float64}(Ep.x, Ep.y, 0.0, Theta), Length, Vstd, "line") # using default constructor
end

function ArcCenter(arc::Geometry)
    if arc.GType != "arc"
        error("input type incorrect: arc needed.")
    end
    radius = abs(1/arc.Sp.cur)
    if arc.Sp.θ > arc.Ep.θ
        δ = π/2
    else
        δ = -π/2
    end
    Sp2 = [radius; 0]
    Rx = [cos(arc.Sp.θ + δ) -sin(arc.Sp.θ + δ); sin(arc.Sp.θ + δ) cos(arc.Sp.θ + δ)]
    Result = [arc.Sp.x; arc.Sp.y] - Rx*Sp2
    Cpoint = Point2D{Float64}(Result[1], Result[2])
    return Cpoint
end

function PointInSector(Ip::Point2D{Float64}, arc::Geometry)
    Cpoint = ArcCenter(arc)
    Vip = [Ip.x-Cpoint.x, Ip.y-Cpoint.y]   # due to matrix malipulation, self-defiend types should not be used
    Vsp = [arc.Sp.x-Cpoint.x, arc.Sp.y-Cpoint.y]
    Vep = [arc.Ep.x-Cpoint.x, arc.Ep.y-Cpoint.y]
    if Distance2D(Vip, Vsp) <= 1.0e-6
        return (true, 0.0)
    elseif Distance2D(Vip, Vep) <= 1.0e-6
        return (true, arc.Length)
    else
        φs = atan(Vsp[2], Vsp[1]) # [-π, π]
        Rx = [cos(φs) sin(φs); -sin(φs) cos(φs)]
        Vip_r = Rx*Vip  # rotate
        Vep_r = Rx*Vep
        φe_r = atan(Vep_r[2], Vep_r[1])
        φi_r = atan(Vip_r[2], Vip_r[1])

        if arc.Sp.cur > 0 # from Sp to Ep: counter-clockwise
            if (φe_r > 0 && 0 < φi_r < φe_r) || (φe_r < 0 && 0 < φi_r < φe_r+2π)
                return (true, φi_r/arc.Sp.cur)
            end
        else # from Sp to Ep: clockwise
            if (φe_r < 0 && φe_r < φi_r < 0) || (φe_r > 0 && φe_r-2π < φi_r < 0)
                return (true, φi_r/arc.Sp.cur)
            end
        end
    end
    return (false, nothing)
end

function Length2XY(Geo::Geometry, Len::Float64)  # convert local 1D coordinate(:length) to global 2D coordinate. The heading direction θ can also be returned if needed
    Δx = 0.0
    Δy = 0.0

    if Geo.GType == "arc"
        radius = 1/Geo.Sp.cur # no need to calculate radius using (x,y,θ), just get it from curvature
        ΔΘ = Len/radius
        Center = ArcCenter(Geo)
        RotMat = [cos(ΔΘ) -sin(ΔΘ); sin(ΔΘ) cos(ΔΘ)]
        Δx, Δy = ([-1 0; 0 -1]+RotMat)*[Geo.Sp.x - Center.x; Geo.Sp.y - Center.y]
        #Δx = sign*radius*( sin(Geo.Sp.θ+sign*ΔΘi) - sin(Geo.Sp.θ) )
        #Δy = sign*radius*( cos(Geo.Sp.θ) - cos(Geo.Sp.θ+sign*ΔΘi) )
    elseif Geo.GType == "line" # straight line
        Δx = (Geo.Ep.x - Geo.Sp.x) * Len/Geo.Length
        Δy = (Geo.Ep.y - Geo.Sp.y) * Len/Geo.Length
    else # clothoid
        # construct later
    end
    return Point2D{Float64}(Geo.Sp.x+Δx, Geo.Sp.y+Δy)
end

function Length2XYθ(Geo::Geometry, Len::Float64)
    Point = Length2XY(Geo, Len)
    θ = Geo.GType == "arc" ? Geo.Sp.θ + Len*Geo.Sp.cur : Geo.Sp.θ
    return Point, θ
end

# the following code block should be improved using better strcuts' structure
function GeosIntersect(Geo1::Geometry, Geo2::Geometry)    # consider using enumerate to handle names
    if Geo1.GType == Geo2.GType
        if Geo1.GType == "line"
            return LineLine(Geo1, Geo2)
        elseif Geo1.GType == "arc"
            return ArcArc(Geo1, Geo2)
        end
    else
        if Geo1.GType == "line" && Geo2.GType == "arc"
            return LineArc(Geo1, Geo2)
        elseif Geo2.GType == "line" && Geo1.GType == "arc"
            return LineArc(Geo2, Geo1)
        end
    end
end


text:
using GenieDevTools
using Test

@testset "GenieDevTools.jl" begin
    # Write your tests here.
end


text:
using Polyopt
using Base.Test

# unconstrained problem
let
    x, y = variables(["x", "y"])
    obj = 4*x^2 + x*y - 4*y^2 - 21//10*x^4 + 4*y^4 + 1//3*x^6

    prob = momentprob(3, obj)
    y, objval, solsta = solve_mosek(prob)

    @test abs(objval - (-1.0316)) < 1e-4
end

# constrained problems


text:
using Documenter
using FLoops
using Literate
using LiterateTest
using LoadAllPackages

LoadAllPackages.loadall(joinpath((@__DIR__), "Project.toml"))

PAGES = [
    "index.md",
    "Tutorials" => [
        "Sequential loops" => "tutorials/sequential.md",
        "Parallel loops" => "tutorials/parallel.md",
    ],
    "How-to guides" => [
        # TODO: Finer grained pages?
        "How to avoid `Box`" => "howto/avoid-box.md",
        "How to do _X_ in parallel?" => "howto/parallel.md",
    ],
    "Reference" => [
        "API" => "reference/api.md",
        "Syntax" => "reference/syntax.md",
        # "Sequential loop" => "reference/sequential.md",
        "Parallelizable reduction (WIP)" => "reference/reduction.md",
    ],
    "Explanation" => [
        "FAQ" => "explanation/faq.md",
        # TODO: discuss design rationales?
    ],
]

let example_dir = joinpath(dirname(@__DIR__), "examples")
    examples = Pair{String,String}[]

    for subpages in PAGES
        subpages isa Pair || continue
        for (_, mdpath) in subpages[2]::Vector
            stem, _ = splitext(mdpath)
            jlpath = joinpath(example_dir, "$stem.jl")
            if !isfile(jlpath)
                @info "`$jlpath` does not exist. Skipping..."
                continue
            end
            push!(examples, jlpath => joinpath(@__DIR__, "src", dirname(mdpath)))
        end
    end

    @info "Compiling example files" examples
    for (jlpath, dest) in examples
        Literate.markdown(jlpath, dest; config = LiterateTest.config(), documenter = true)
    end
end

makedocs(;
    sitename = "FLoops",
    format = Documenter.HTML(),
    modules = [FLoops],
    pages = PAGES,
    checkdocs = :exports,
    strict = lowercase(get(ENV, "CI", "false")) == "true",
)

deploydocs(; repo = "github.com/JuliaFolds/FLoops.jl", push_preview = true)


text:
using Strategems, Temporal, Indicators, Base.Dates
using Base.Test

# define universe and gather data
assets = ["CHRIS/CME_CL1", "CHRIS/CME_RB1"]
universe = Universe(assets)
function datasource(asset::String; save_downloads::Bool=true)::TS
    savedata_path = Pkg.dir("Strategems", "data/$asset.csv")
    if isfile(savedata_path)
        return Temporal.tsread(savedata_path)
    else
        X = quandl(asset)
        if save_downloads
            Temporal.tswrite(X, savedata_path)
        end
        return X
    end
end
gather!(universe, source=datasource)

# define indicators and parameter space
arg_names = [:fastlimit, :slowlimit]
arg_defaults = [0.5, 0.05]
arg_ranges = [0.01:0.01:0.99, 0.01:0.01:0.99]
paramset = ParameterSet(arg_names, arg_defaults, arg_ranges)
f(x; args...) = Indicators.mama(x; args...)
indicator = Indicator(f, paramset)

# define signals that will trigger trading decisions
siglong = @signal MAMA ↑ FAMA
sigshort = @signal MAMA ↓ FAMA
sigexit = @signal MAMA .== FAMA

# define the trading rules
longrule = @rule siglong → long 100
shortrule = @rule sigshort → short 100
exitrule = @rule sigexit → liquidate 1.0
rules = (longrule, shortrule, exitrule)

# run strategy
strat = Strategy(universe, indicator, rules)
backtest!(strat)
optimize!(strat, samples=10)


text:
## ExoplanetsSysSim/examples/hsu_et_al/sim_param_plots.jl
## (c) 2017 Keir A. Ashby
# Takes any number of SimParams and, calculates observed and physical catalogs with summary statistics and returns:
#   1.A period-radius distribution of the observed, physical, and true catalogs for each sim_param individually
#   2.A period-radius distribution of the observed and true catalogs for each sim_param individually
#   3.An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
#   4.An overlaid comparison period-radius distribution of the observed and true catalogs for every input sim_param
#   5.The percentage of observed planets to true planets in each period-radius bin
#   6.CDF plots for radii and periods

using PyPlot
using ExoplanetsSysSim
using Combinatorics

include(joinpath(Pkg.dir(),"ExoplanetsSysSim","src","constants.jl"))
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2017", "abc_setup_christiansen.jl"))
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2017", "christiansen_func.jl"))
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2017", "calc_sum_stats_pr.jl"))

#see general file definition at the top for a description of this function's purpose.
function sim_param_plots(have_sim_params::Bool; sim_arr::Array{Any,1} = [], num_param::Int64 = 0, print_rates=false) 	#if the user already has SimParams to put in, the first input will be true and the second will be an array (even if there's only one) of your SimParams. If you don't already have one, put false as your first input, and how many SimParams you'd like to generate. If you want to see the ratios of observed to kepler planets in all the bins, set print_rates=true)
  if have_sim_params == false 						#we will make SimParams for the user if she doesn't already have one
    sim_arr = setup_sim_arr(num_param) 					#returns an array of SimParams
  end
  sim_true = setup_sim_arr(1)[1] 					#generates a SimParam that will be used to get the true Kepler catalog.
  ss_arr, ss_true = setup_sum_stats(sim_arr, sim_true) 			#generates physical and observed catalogs and returns summary statistics for each SimParam and for the true Kepler catalog
  num_cats = length(ss_arr)
  n=1									#iterative integer to keep track of figures
  pr_lists = make_pr_lists(num_cats, ss_arr, ss_true)			#Makes lists of all the periods and radii for each catalog
  n = all_cats_pr(n, num_cats, pr_lists)				#A period-radius distribution of the observed, physical, and true catalogs for each sim_param individually
  n = all_cats_pr(n, num_cats, pr_lists, false)				#A period-radius distribution of the observed and true catalogs for each sim_param individually
  if num_cats > 1
    n = compare_2_obs_pr(n, num_cats, pr_lists)				#An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
  end
  if num_cats > 2
    n = compare_all_obs_pr(n, num_cats, pr_lists)			#An overlaid comparison period-radius distribution of the observed and true catalogs for every input sim_param
  end
  ratios = pr_rates(pr_lists["periods_obs_1"], pr_lists["radii_obs_1"], pr_lists["periods_true"], pr_lists["radii_true"], pr_lists["p_lim_arr"], pr_lists["r_lim_arr"]) #The percentage of observed planets to true planets in each period-radius bin
  n = cdf_plot(n, pr_lists["periods_obs_1"], pr_lists["periods_true"], pr_lists["p_lim_arr"], "Periods") 	#CDF plots for periods
  n = cdf_plot(n, pr_lists["radii_obs_1"], pr_lists["radii_true"], pr_lists["r_lim_arr"], "Radii")		#CDF plots for radii
  return(ratios)							#Normally, pr_rates would just be called last so it printed ratios at the end automatically. However, cdf_plot for radii alters the radii_obs_1 list so that it doesn't work in pr_rates. I'm working on how to alter the list in cdf_plot without altering it globally.
end

function setup_sim_arr(num_param::Int64 = 1) 				#takes an integer and returns an array of that many Sim Params
  sim_arr = Any[]
  for i in 1:num_param
    sim_param = setup_sim_param_christiansen() 				#sim-param used in christiansen_func.jl
    hsu_final_rates = [0.039  0.4    2.6   5.1   5.4  3.9  1.6  56  87
		       0.1    0.3    1     3.2   3.6  7.9  1.7  5.5 21
		       0.088  0.18   0.52  1.5   2.8  3.1  4.1  3.2 8.5
		       0.053  0.13   0.63  1.6   2.4  2.4  2.3  4.4 7.6
		       0.047  0.13   0.74  1.1   1.9  1.9  2.3  2.7 2.2
		       0.036  0.057  0.38  0.8   1.6  1.2  2    2.4 3.5
		       0.0096 0.059  0.29  1.1   2.5  3.7  3.5  3.6 6.4
   		       0.007  0.042  0.22  0.93  2.3  2.8  3.5  3.1 3.2
		       0.0044 0.029  0.14  0.6   1.4  1.9  2.3  3.4 5
		       0.0086 0.012  0.13  0.23  0.48 0.82 0.91 1.5 1.8
		       0.0033 0.0068 0.057 0.1   0.19 0.23 0.76 0.7 1.9
       		       0.0045 0.012  0.082 0.16  0.23 0.41 0.32 1.2 2.4
		       0.0045 0.031  0.12  0.075 0.16 0.18 0.33 0.5 0.59]*0.01	#occurence rates over all bins drawn from results in Hsu et al 2017
    hsu_radii = [0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 6, 8, 12, 16]	#limits for the radius bins
    hsu_periods = [0.5, 1.25, 2.5, 5, 10, 20, 40, 80, 160, 320]			#limits for the period bins
    update_param(sim_param, "obs_par", hsu_final_rates)
    update_param(sim_param, "r_lim_arr", hsu_radii*ExoplanetsSysSim.earth_radius)
    update_param(sim_param, "p_lim_arr", hsu_periods)
    update_param(sim_param, "num_targets_sim_pass_one", 150000)
    update_param(sim_param, "num_kepler_targets", 150000)
    push!(sim_arr,sim_param)
  end
  return(sim_arr)
end

#takes an array of Sim Params and returns summary statistics for each. It will also take one Sim Param and make summary statistics for the true Kepler Catalog
function setup_sum_stats(sim_arr::Array{Any,1}, sim_true::ExoplanetsSysSim.SimParam)
  setup_MES()
  cat_true = setup_actual_planet_candidate_catalog(setup_star_table_christiansen(sim_true), sim_true)
  ss_true = calc_sum_stats_pr(cat_true, sim_true, true)
  ss_arr = Array{Any}(length(sim_arr))
  for i in 1:length(sim_arr)
    cat_phys = generate_kepler_physical_catalog(sim_arr[i])
    cat_obs = ExoplanetsSysSim.observe_kepler_targets_single_obs(cat_phys, sim_arr[i])
    ss_arr[i] = calc_sum_stats_pr(cat_obs, sim_arr[i], false, cat_phys)
  end
  return(ss_arr, ss_true)
end

function make_pr_lists(num_cats, ss_arr::Array{Any,1}, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics)	#Makes lists of all the periods and radii for each catalog
  p_lim_arr = [0.5, 1.25, 2.5, 5, 10, 20, 40, 80, 160, 320]				#bin boundaries for periods
  r_lim_arr = [0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 6, 8, 12, 16, 20]		#bin boundaries for radii
  pr_lists = Dict("p_lim_arr"=> p_lim_arr, "r_lim_arr"=> r_lim_arr)
  for i in 1:num_cats
    pr_lists["periods_obs_$i"] = ss_arr[i].stat["period_obs_list"] 			#all the periods in this observable catalog
    pr_lists["radii_obs_$i"] = ss_arr[i].stat["radius_obs_list"]./earth_radius 		#all the radii in this observable catalog
    pr_lists["periods_phys_$i"] = ss_arr[i].stat["period_phys_list"] 			#all the periods in this physical catalog
    pr_lists["radii_phys_$i"] = ss_arr[i].stat["radius_phys_list"]./earth_radius 	#all the radii in this physical catalog
    if length(pr_lists["periods_phys_$i"])==0 || length(pr_lists["periods_obs_$i"])==0 	#test to make sure we have something to plot
      println("One of your catalogs is empty")
      return()
    end
  end
  pr_lists["periods_true"] = ss_true.stat["period_list"] 				#all the periods in this true catalog
  pr_lists["radii_true"] = ss_true.stat["radius_list"]./earth_radius 			#all the radii in this true catalog
  return(pr_lists)
end

#A period-radius distribution of the observed, physical, and true catalogs for each sim_param individually
function all_cats_pr(n::Int64, num_cats::Int64, pr_lists::Dict{String,Array{Float64,1}}, plot_phys = true)
  for i in 1:num_cats
    figure(n)
    if plot_phys == true
      loglog(pr_lists["periods_phys_$i"], pr_lists["radii_phys_$i"], ".", ms = 0.5, color = "orange", label = "Physical")
    end
    loglog(pr_lists["periods_obs_$i"], pr_lists["radii_obs_$i"], ".", ms = 3, label = "Simulated Observed")
    loglog(pr_lists["periods_true"], pr_lists["radii_true"], "x", ms = 2, color = "purple", label = "Kepler Observed")
    title("Periods vs. Radii $i", fontsize = 18)
    axis([0.0, 350.0, 0.4, 22.0])
    xticks(pr_lists["p_lim_arr"], pr_lists["p_lim_arr"])
    yticks(pr_lists["r_lim_arr"], pr_lists["r_lim_arr"])
    xlabel("Periods (in days)", fontsize=14)
    ylabel("Radius (in Earth radii)", fontsize = 14)
    legend(loc = "lower right", fontsize = 10)
    if plot_phys == true
     savefig("all_cats_pr_$i.pdf")
    else
     savefig("obs_cats_pr_$i.pdf")
    end
    n+=1
  end
  return(n)
end

#An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
function compare_2_obs_pr(n::Int64, num_cats::Int64, pr_lists::Dict{String,Array{Float64,1}})
  for (a,b) in combinations(collect(1:num_cats),2)
    figure(n)
    loglog(pr_lists["periods_true"], pr_lists["radii_true"], "x", ms = 2, color = "purple", label = "Kepler Observed")
    loglog(pr_lists["periods_obs_$a"], pr_lists["radii_obs_$a"], ".", ms = 3, label = "Simulated Observed: $a")
    loglog(pr_lists["periods_obs_$b"], pr_lists["radii_obs_$b"], ".", ms = 3, label = "Simulated Observed: $b")
    title("Periods vs. Radii Comparison: $a & $b", fontsize = 18)
    axis([0.0, 350.0, 0.4, 22.0])
    xticks(pr_lists["p_lim_arr"], pr_lists["p_lim_arr"])
    yticks(pr_lists["r_lim_arr"], pr_lists["r_lim_arr"])
    xlabel("Periods (in days)", fontsize = 14)
    ylabel("Radius (in Earth radii)", fontsize = 14)
    legend(loc = "lower right", fontsize = 10)
    n+=1
#    savefig("compare_2_obs_pr_$a_&_$b.pdf")
  end
  return(n)
end

#An overlaid comparison period-radius distribution of the observed and true catalogs for every possible combination of two input sim_params
function compare_all_obs_pr(n::Int64, num_cats::Int64, pr_lists::Dict{String,Array{Float64,1}})
  figure(n)
  loglog(pr_lists["periods_true"], pr_lists["radii_true"], "x", ms = 2, color = "purple", label = "Kepler Observed")
  for i in 1:num_cats
    loglog(pr_lists["periods_obs_$i"], pr_lists["radii_obs_$i"], ".", ms = 3, label = "Simulated Observed: $i")
  end
  title("Periods vs. Radii Comparison: All", fontsize = 18)
  axis([0.0, 350.0, 0.4, 22.0])
  xticks(pr_lists["p_lim_arr"], pr_lists["p_lim_arr"])
  yticks(pr_lists["r_lim_arr"], pr_lists["r_lim_arr"])
  xlabel("Periods (in days)", fontsize=14)
  ylabel("Radius (in Earth radii)", fontsize = 14)
  legend(loc = "lower right", fontsize = 10)
#  savefig("compare_all_obs_pr.pdf")
  n+=1
  return(n)
end

#CDF plots for radii and periods
function cdf_plot(n::Int64, obs_list::Array{Float64,1}, true_list::Array{Float64,1}, lim_arr::Array{Float64,1}, cdf_type::String)
  figure(n)
  obs_list = trim_catalog(obs_list, lim_arr[1], lim_arr[length(lim_arr)])	#we want to restrict our catalog to stay within our bin ranges
  true_list = trim_catalog(true_list, lim_arr[1], lim_arr[length(lim_arr)])	#we want to restrict our catalog to stay within our bin ranges
  semilogx(obs_list, collect(linspace(0.0,1.0,length(obs_list))), label = "Simulated Observed")
  semilogx(true_list, collect(linspace(0.0,1.0,length(true_list))), label = "Kepler Observed")
  title("Cumulative Distribution of $cdf_type")
  xlabel("Periods (in days)")
  ylabel("Probability")
  legend(loc = "lower right", fontsize = 10)
#  savefig("cdf_$cdf_type.pdf")
  n+=1
  return(n)
end

#This is a function for testing the accuracy of our model for how many planets Kepler would see. Returns a 2X2 array  of observed/true rates in fraction or decimal form
function pr_rates(periods_obs::Array{Float64,1}, radii_obs::Array{Float64,1}, periods_true::Array{Float64,1}, radii_true::Array{Float64,1}, p_lim_arr::Array{Float64,1}, r_lim_arr::Array{Float64,1}, percent::Bool = true) #if percent is false, the rates will be printed as fractions
  rates_obs = zeros(length(p_lim_arr)-1,length(r_lim_arr)-1)			#table to hold the number of observed planets in each bin
  rates_true = zeros(length(p_lim_arr)-1,length(r_lim_arr)-1)			#table to hold the number of true planets in each bin
  ratios = Array{String}(length(p_lim_arr)-1,length(r_lim_arr)-1)		#table to hold the ratios of observed/true planets in each bin
  for i in 1:length(p_lim_arr)-1						#for each column of periods
    f(p) = p>p_lim_arr[i] && p<p_lim_arr[i+1]
    p_bins_obs = find(f,periods_obs)						#this will hold all the indices of planets that fall into the ith period bin
    r_bins_obs = zeros(length(p_bins_obs))					#this will hold the radius corresponding to each exoplanet whose index is held in p_bins_obs
    for j in 1:length(p_bins_obs)
      r_bins_obs[j]=radii_obs[p_bins_obs][j]
    end
    for k in 1:length(r_lim_arr)-1						#for each row of radii in the ith column of periods
      g(r) = r>r_lim_arr[k] && r<r_lim_arr[k+1]
      rates_obs[i,length(r_lim_arr)-k] = convert(Int64,length(find(g,r_bins_obs)))
    end
  end
#  for h in 1:length(r_lim_arr)-1 						#shows rates_obs
#    println(rates[:,h])
#  end

  for i in 1:length(p_lim_arr)-1						#for each column of periods
    f(p) = p>p_lim_arr[i] && p<p_lim_arr[i+1]
    p_bins_true = find(f,periods_true)						#this will hold all the indices of planets that fall into the ith period bin
    r_bins_true = zeros(length(p_bins_true))					#this will hold the radius corresponding to each exoplanet whose index is held in p_bins_true
    for j in 1:length(p_bins_true)
      r_bins_true[j]=radii_true[p_bins_true][j]
    end
    for k in 1:length(r_lim_arr)-1						#for each row of radii in the ith column of periods
      g(r) = r>r_lim_arr[k] && r<r_lim_arr[k+1]
      rates_true[i,length(r_lim_arr)-k] = convert(Int64,length(find(g,r_bins_true)))
    end
  end
#  for h in 1:length(r_lim_arr)-1						#shows rates_true
#    println(rates_true[:,h])
#  end

  for i in 1:length(p_lim_arr)-1
    for j in 1:length(r_lim_arr)-1
      if percent == true
	if rates_obs[i,j] == 0.0 && rates_true[i,j] == 0.0
	  ratios[i,j] = "NaN"
        elseif rates_true[i,j] == 0.0
	  ratios[i,j] = "Inf"
	else
	  ratios[i,j] = string(floor(Int64,rates_obs[i,j]/rates_true[i,j]*100), "%")
	end
      else    									#returns the ratios in fraction form
        ratios[i,j] = string(convert(Int64,rates_obs[i,j]), "/", convert(Int64,rates_true[i,j]))
      end
    end
  end
  return ratios
end

function trim_catalog(list::Array{Float64,1},min::Float64,max::Float64)		#removes all values in a list to stay within bin ranges
  sort!(list)
  while list[1] < min
    shift!(list)
  end
  while list[length(list)] > max
    pop!(list)
  end
  return(list)
end

function test_period_radius_plots()
  sim_param_plots(false, num_param = 3)
end