text:
#include("../src/ttv.jl")
#include("/Users/ericagol/Computer/Julia/regress.jl")

@testset "ttv_elements" begin

# This routine takes derivative of transit times with respect
# to the initial orbital elements.
#n = 8
n = 3
n_body = n
t0 = 7257.93115525-7300.0
#h  = 0.12
h  = 0.04
#tmax = 600.0
#tmax = 800.0
#tmax = 1000.0
tmax = 10.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')
elements[:,3] .-= 7300.0
# Make masses of planets bigger
#elements[2,1] *= 10.0
#elements[3,1] *= 10.0

ntt = zeros(Int64,n)

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
dtdq0 = zeros(n,maximum(ntt),7,n)
tt  = zeros(n,maximum(ntt))
tt1 = zeros(n,maximum(ntt))
tt2 = zeros(n,maximum(ntt))
tt3 = zeros(n,maximum(ntt))
tt4 = zeros(n,maximum(ntt))
tt8 = zeros(n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
count1 = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12
dq = ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
# Now call with half the timestep:
count2 = zeros(Int64,n)
count3 = zeros(Int64,n)
dq = ttv_elements!(n,t0,h/10.,tmax,elements,tt2,count2,0.0,0,0,rstar)

mask = zeros(Bool, size(dtdq0))
for jq=1:n_body
  for iq=1:7
    if iq == 7; ivary = 1; else; ivary = iq+1; end  # Shift mass variation to end
    for i=2:n
      for k=1:count2[i]
        # Ignore inclination & longitude of nodes variations:
        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
          mask[i,k,iq,jq] = true
        end
      end
    end
  end
end

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dtdelements0 = zeros(n,maximum(ntt),7,n)
dtdelements0 = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar)
dtdq2 = zeros(n,maximum(ntt),7,n)
dtdelements2 = zeros(n,maximum(ntt),7,n)
dtdelements2 = ttv_elements!(n,t0,h/2.,tmax,elements,tt2,count,dtdq2,rstar)
dtdq4 = zeros(n,maximum(ntt),7,n)
dtdelements4 = zeros(n,maximum(ntt),7,n)
dtdelements4 = ttv_elements!(n,t0,h/4.,tmax,elements,tt4,count,dtdq4,rstar)
dtdq8 = zeros(n,maximum(ntt),7,n)
dtdelements8 = zeros(n,maximum(ntt),7,n)
dtdelements8 = ttv_elements!(n,t0,h/8.,tmax,elements,tt8,count,dtdq8,rstar)
#println("Maximum error on derivative: ",maximum(abs.(dtdelements0-dtdelements2)))
#println("Maximum error on derivative: ",maximum(abs.(dtdelements2-dtdelements4)))
#println("Maximum error on derivative: ",maximum(abs.(dtdelements4-dtdelements8)))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdelements0[mask])-asinh.(dtdelements2[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdq0)-asinh.(dtdq2))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdelements2[mask])-asinh.(dtdelements4[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdq2)-asinh.(dtdq4))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdelements4[mask])-asinh.(dtdelements8[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdq4)-asinh.(dtdq8))))
#read(STDIN,Char)

# Check that this is working properly:
#for i=1:n
#  for j=1:count2[i]
#    println(i," ",j," ",tt[i,j]," ",tt2[i,j]," ",tt[i,j]-tt2[i,j]," ",tt1[i,j]-tt2[i,j])
#  end
#end
#read(STDIN,Char)

# Compute derivatives numerically:
#nq = 15
# This "summarizes" best numerical derivative:
dtdelements0_num = zeros(BigFloat,n,maximum(ntt),7,n)

# Compute derivatives with BigFloat for additional precision:
elements0 = copy(elements)
#delement = big.([1e-15,1e-15,1e-15,1e-15,1e-15,1e-15,1e-15])
#dq0 = big(1e-20)
dq0 = big(1e-10)
t0big = big(t0); tmaxbig = big(tmax); hbig = big(h)
zilch = big(0.0)
# Compute the transit times in BigFloat precision:
tt_big = big.(tt); elementsbig = big.(elements0); rstarbig = big(rstar)
dqbig = ttv_elements!(n,t0big,hbig,tmaxbig,elementsbig,tt_big,count,big(0.0),0,0,rstarbig)

mask = zeros(Bool, size(dtdq0))
for jq=1:n_body
  for iq=1:7
    if iq == 7; ivary = 1; else; ivary = iq+1; end  # Shift mass variation to end
    for i=2:n
      for k=1:count2[i]
        # Ignore inclination & longitude of nodes variations:
        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
          mask[i,k,iq,jq] = true
        end
      end
    end
  end
end

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dtdelements0 = zeros(n,maximum(ntt),7,n)
dtdelements0 = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar)
dtdq2 = zeros(n,maximum(ntt),7,n)
dtdelements2 = zeros(n,maximum(ntt),7,n)
dtdelements2 = ttv_elements!(n,t0,h/2.,tmax,elements,tt2,count,dtdq2,rstar)
dtdq4 = zeros(n,maximum(ntt),7,n)
dtdelements4 = zeros(n,maximum(ntt),7,n)
dtdelements4 = ttv_elements!(n,t0,h/4.,tmax,elements,tt4,count,dtdq4,rstar)
dtdq8 = zeros(n,maximum(ntt),7,n)
dtdelements8 = zeros(n,maximum(ntt),7,n)
dtdelements8 = ttv_elements!(n,t0,h/8.,tmax,elements,tt8,count,dtdq8,rstar)
#println("Maximum error on derivative: ",maximum(abs.(dtdelements0-dtdelements2)))
#println("Maximum error on derivative: ",maximum(abs.(dtdelements2-dtdelements4)))
#println("Maximum error on derivative: ",maximum(abs.(dtdelements4-dtdelements8)))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdelements0[mask])-asinh.(dtdelements2[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdq0)-asinh.(dtdq2))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdelements2[mask])-asinh.(dtdelements4[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdq2)-asinh.(dtdq4))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdelements4[mask])-asinh.(dtdelements8[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtdq4)-asinh.(dtdq8))))
#read(STDIN,Char)

# Check that this is working properly:
#for i=1:n
#  for j=1:count2[i]
#    println(i," ",j," ",tt[i,j]," ",tt2[i,j]," ",tt[i,j]-tt2[i,j]," ",tt1[i,j]-tt2[i,j])
#  end
#end
#read(STDIN,Char)

# Compute derivatives numerically:
#nq = 15
# This "summarizes" best numerical derivative:
dtdelements0_num = zeros(BigFloat,n,maximum(ntt),7,n)

# Compute derivatives with BigFloat for additional precision:
elements0 = copy(elements)
#delement = big.([1e-15,1e-15,1e-15,1e-15,1e-15,1e-15,1e-15])
#dq0 = big(1e-20)
dq0 = big(1e-10)
tt2 = big.(tt2)
tt3 = big.(tt3)
t0big = big(t0); tmaxbig = big(tmax); hbig = big(h)
zilch = big(0.0)
# Compute the transit times in BigFloat precision:
tt_big = big.(tt); elementsbig = big.(elements0)
# Now, compute derivatives numerically:
for jq=1:n_body
  for iq=1:7
    elementsbig = big.(elements0)
#    dq0 = delement[iq]; if jq==1 && iq==7 ; dq0 = big(1e-10); end  # Vary mass of star by a larger factor
    if iq == 7; ivary = 1; else; ivary = iq+1; end  # Shift mass variation to end
    elementsbig[jq,ivary] += dq0
    dq_plus = ttv_elements!(n,t0big,hbig,tmaxbig,elementsbig,tt2,count2,zilch,0,0,big(rstar))
    elementsbig[jq,ivary] -= 2dq0
    dq_minus = ttv_elements!(n,t0big,hbig,tmaxbig,elementsbig,tt3,count2,zilch,0,0,big(rstar))
    #xm,vm = init_nbody(elements,t0,n_body)
    for i=2:n
      for k=1:count2[i]
        # Compute double-sided derivative for more accuracy:
        dtdelements0_num[i,k,iq,jq] = (tt2[i,k]-tt3[i,k])/(2dq0)
        # Ignore inclination & longitude of nodes variations:
        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
          mask[i,k,iq,jq] = true
        end
      end
    end
  end
end

#println("Max diff dtdelements: ",maximum(abs.(dtdelements0[mask]./dtdelements0_num[mask].-1.0)))
println("Max diff asinh(dtdelements): ",maximum(abs.(asinh.(dtdelements0[mask])-asinh.(dtdelements0_num[mask]))))

#ntot = 0
#diff_dtdelements0 = zeros(n,maximum(ntt),7,n)
#for i=1:n, j=1:count[i], k=1:7, l=2:n
#  diff_dtdelements0[i,j,k,l] = abs(dtdelements0[i,j,k,l]-convert(Float64,dtdelements0_num[i,j,k,l]))
#  ntot +=1
#end

#using PyPlot
#
#nderiv = n^2*7*maximum(ntt)
##mask[:,:,2,:] = false
#loglog(abs.(reshape(dtdelements0,nderiv)),abs.(reshape(convert(Array{Float64,4},dtdelements0_num),nderiv)),".")
#axis([1e-6,1e2,1e-12,1e2])
#loglog(abs.(reshape(dtdelements0,nderiv)),abs.(reshape(diff_dtdelements0,nderiv)),".")
#println("Maximum error: ",maximum(diff_dtdelements0))

using PyPlot

# Make a plot of the fractional errors:
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdelements0_num[i,2:count1[i],k,l])) > 0
    loglog(tt[i,2:count1[i]].-tt[i,1],convert(Array{Float64,1},abs.(dtdelements0[i,2:count1[i],k,l]./dtdelements0_num[i,2:count1[i],k,l].-1.)))
  end
end

clf()
# Plot the difference in the TTVs:
for i=2:3
  diff1 = convert(Array{Float64,1},abs.(tt1[i,2:count1[i]]./tt_big[i,2:count1[i]].-1.0));
  loglog(tt[i,2:count1[i]].-tt[i,1],diff1);
end
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdelements0_num[i,2:count1[i],k,l])) > 0 && k != 5 && k != 6
    diff1 = convert(Array{Float64,1},abs.(dtdelements0[i,2:count1[i],k,l]./dtdelements0_num[i,2:count1[i],k,l].-1));
    diff2 = convert(Array{Float64,1},abs.(asinh.(dtdelements0[i,2:count1[i],k,l])-asinh.(dtdelements0_num[i,2:count1[i],k,l])));
    loglog(tt[i,2:count1[i]].-tt[i,1],diff1);
#    loglog(tt[i,2:count1[i]].-tt[i,1],diff2,".");
    println(i," ",k," ",l," frac error: ",convert(Float64,maximum(diff1))," asinh error: ",convert(Float64,maximum(diff2))); #read(STDIN,Char);
  end
end

# Plot a line that scales as time^{3/2}:

loglog([1.0,1024.0],1e-09*[1,2^15],":")
loglog([1.0,1024.0],1e-12*[1,2^15],":")
loglog([1.0,1024.0],1e-15*[1,2^15],":")

#@test isapprox(dtdelements0[mask],dtdelements0_num[mask];norm=maxabs)
@test isapprox(asinh.(dtdelements0[mask]),asinh.(dtdelements0_num[mask]);norm=maxabs)
#unit = ones(dtdelements0[mask])
#@test isapprox(dtdelements0[mask]./dtdelements0_num[mask],unit;norm=maxabs)
#end

## Make a plot of some TTVs:
#
#fig,axes = subplots(4,2)
#
#for i=2:8
#  ax = axes[i-1]
#  fn = zeros(Float64,2,count1[i])
#  sig = ones(count1[i])
#  tti1 = tt1[i,1:count1[i]]
#  tti2 = tt2[i,1:count2[i]]
#  for j=1:count1[i]
#    fn[1,j] = 1.0
#    fn[2,j] = round(Int64,(tti1[j]-elements[i,3])/elements[i,2])
#  end
#  coeff,cov = regress(fn,tti1,sig)
#  tt_ref1 = coeff[1]+coeff[2]*fn[2,:]
#  ttv1 = (tti1.-tt_ref1)*24.*60.
#  coeff,cov = regress(fn,tti2,sig)
#  tt_ref2 = coeff[1]+coeff[2]*fn[2,:]
#  ttv2 = (tti2.-tt_ref2)*24.*60.
#  ax[:plot](tti1,ttv1)
##  ax[:plot](tti2,ttv2)
#  ax[:plot](tti2,((ttv1.-ttv2)-mean(ttv1.-ttv2)))
#  println(i," ",coeff," ",elements[i,2:3]," ",coeff[1]-elements[i,3]," ",coeff[2]-elements[i,2])
#  println(i," ",maximum(ttv1-ttv2-mean(ttv1-ttv2))*60.," sec ", minimum(ttv1-ttv2-mean(ttv1-ttv2))*60.," sec" )
end


text:
using ProgressMeter
using JLD
using LinearAlgebra

function residualInverseIteration(iMax::Int, ϵ::Float64, b::Boundary, n::Int, k0::Complex{Float64}, x0::Array{Complex{Float64},1}; useWiersigKShift::Bool = false)
  residualInverseIteration(iMax, ϵ, b, n, k0, x0, nothing, 0; useWiersigKShift = useWiersigKShift)
end

function residualInverseIteration(iMax::Int, ϵ::Float64, b::Boundary, n::Int, k0::Complex{Float64}, x0::Array{Complex{Float64},1}, progressMeter, iter; useWiersigKShift::Bool = false)
  if progressMeter == nothing
    progressMeter = ProgressThresh(ϵ, "Finding resonance:")
  end

  @debug "Starting iteration, attempt to converge k with ϵ = $ϵ"

  numChoices = min(n,20)
  # Choose from a random large value in x0, don't use the same index each iteration
  # Not sure if this helps in any way
  xMaxInd = sortperm(-abs.(x0[1:n]))[rand(1:numChoices)]
  # xMaxInd = indmax(abs(x0))

  # Construct normalization vector
  e = [j == xMaxInd ? 1.0+0.0im : 0.0+0.0im for j in 1:4n]
  # Normalize incoming x0
  x0 = x0 ./ (e' * x0)[1]
  @debug "max idx = $(xMaxInd), x = $(x0[xMaxInd]), e=$(Array(e)), x = $(x0)"

  @debug "Doing initial precalculation for k=$k0"
  bpc0 = doPrecalculation(b, n, k0)
  @debug "Calculating A matrix"
  A0 = A(b, n, k0, bpc0)
  @debug "Inverting A matrix"
  A0inv = inv(A0)

  xl = x0
  kl = k0
  ks = [kl]

  for l in 1:iMax
    iter += 1
    @debug "l = $l, doing precalculation for k=$kl"
    bpcl = doPrecalculation(b, n, kl)
    @debug "Calculating A' derivative matrix"
    A′l = A′(b, n, kl, bpcl)
    @debug "Calculating A(kl) matrix"
    Akl = A(b, n, kl, bpcl)

    @debug "e' * (A0inv*Akl*xl)" dote = e' * (A0inv*Akl*xl)
    @debug "e' * (A0inv*A′l*xl)" dotA = e' * (A0inv*A′l*xl)

    if useWiersigKShift
      # From Wiersig's paper eq (37), this seems to work slightly more slowly and less
      # accurately than Heider's version. I'm curious if there are ever cases
      # where it could work better?
      Aklinv = inv(Akl)
      dk = - 1/(trace(Aklinv*A′l))
      dkHeider = - ((e' * (A0inv*Akl*xl)) / (e' * (A0inv*A′l*xl)))[1]
      @debug "Heider's dk would have been = $(dkHeider)"
    else
      dk = - ((e' * (A0inv*Akl*xl)) / (e' * (A0inv*A′l*xl)))[1]
    end
    @debug "dk: $(dk)"
    kNext = kl+dk
    @debug "k(l+1): $kNext"
    push!(ks, kNext)
    @debug "Doing l+1 precalculation for k=$kNext"
    bpcNext = doPrecalculation(b, n, kNext)
    @debug "Calculating A matrix"
    Anext = A(b, n, kNext, bpcNext)
    # Calculate residuum
    rl = Anext*xl
    # Use that to determine dx
    # Had this from Heider's paper but since we've already calcuated A0inv, seems like
    # just as good to use it?
    # dxl = A0 \ rl
    dxl = A0inv*rl
    @debug "|dx_l| = $(norm(dxl))"

    #Adjust x vector
    xNextBar = xl - dxl
    xNext = xNextBar ./ ((e' * xNextBar)[1])

    @debug "Checking |A(kl).x_l+1| ~= 0"
    beforeShiftQuality = norm(Anext*xl)/norm(xl)
    afterShiftQuality = norm(Anext*xNext)/norm(xNext)
    @debug "|A(kl+1).x_l|/|x_l|=$beforeShiftQuality"
    @debug "|A(kl+1).x_l+1|/|x_l+1|=$afterShiftQuality"

    kl = kNext
    xl = xNext
    ProgressMeter.update!(progressMeter, abs.(ks[end]-ks[end-1]), showvalues = [
      (:iter,iter),
      (:dk,dk),
      (:kl,kl),
      ("|A(kl+1).x_l|/|x_l|", beforeShiftQuality),
      ("|A(kl+1).x_l+1|/|x_l+1|", afterShiftQuality)
    ])
    if !correctorCondition(ks)
      break
    end
  end

  if length(ks) > 1 && abs.(ks[end]-ks[end-1]) < ϵ
    @debug "k converged to $kl"
    # Clear memory
    bpc0 = 0
    A0 = 0
    A0inv = 0
    x0 = 0
    bpcF = doPrecalculation(b, n, kl)
    Af = A(b, n, kl, bpcF)
    fEigen = eigen(Af)
    fVals = fEigen.values
    fVecs = fEigen.vectors
    perm = sortperm(abs.(fVals))
    @info "Min eigenvalue = $(fVals[perm[1]])"
    xl = Array(xl)
    xfDiff = norm(fVecs[:,perm[1]]-Array(xl))
    @info "Diff in field = $xfDiff"
    (kl, Array(xl))
  elseif (iter < iMax)
    @debug "Restart RII with shift $kl and current xl"
    try
      # Clear memory
      bpc0 = 0
      A0 = 0
      A0inv = 0
      x0 = 0
      residualInverseIteration(iMax, ϵ, b, n, kl, xl, progressMeter, iter; useWiersigKShift = useWiersigKShift)
    catch e
      @error ("Got exception in RII") exception=e
      error("k did not converge after $iMax iterations, last kl = $kl")
    end
  else
    error("k did not converge after $iMax iterations, last kl = $kl")
  end
end

function correctorCondition(ks)
  if size(ks)[1] == 1
    true
  elseif abs.((ks[end]-ks[end-1])/ks[end]) < 1/10
    @debug "Hit corrector condition, kl=$(ks[end])"
    false
  else
    true
  end
end

ext: jl


text:
"""
The `Richardson` module provides a function `extrapolate` that
extrapolates a given function `f(x)` to `f(x0)`, evaluating
`f` only  at a geometric sequence of points `> x0`
(or optionally `< x0`).

The key algorithm is Richardson extrapolation using a Neville—Aitken
tableau, which adaptively increases the degree of an extrapolation
polynomial until convergence is achieved to a desired tolerance
(or convergence stalls due to e.g. floating-point errors).  This
allows one to obtain `f(x0)` to high-order accuracy, assuming
that `f(x0+h)` has a Taylor series or some other power
series in `h`.
"""
module Richardson

using LinearAlgebra

export extrapolate, extrapolate!

"""
    extrapolate(f, h; contract=0.125, x0=zero(h), power=1,
                      atol=0, rtol=atol>0 ? 0 : sqrt(ε),
                      maxeval=typemax(Int), breaktol=2)

Extrapolate `f(x)` to `f₀ ≈ f(x0)`, evaluating `f` only at `x > x0` points
(or `x < x0` if `h < 0`) using Richardson extrapolation starting at
`x=x₀+h`.  It returns a tuple `(f₀, err)` of the estimated `f(x0)`
and an error estimate.

The return value of `f` can be any type supporting `±` and `norm`
operations (i.e. a normed vector space).
Similarly, `h` and `x0` can be in any normed vector space,
in which case `extrapolate` performs Richardson extrapolation
of `f(x0+s*h)` to `s=0⁺` (i.e. it takes the limit as `x` goes
to `x0` along the `h` direction).

On each step of Richardson extrapolation, it shrinks `x-x0` by
a factor of `contract`, stopping when the estimated error is
`< max(rtol*norm(f₀), atol)`, when the estimated error
increases by more than `breaktol` (e.g. due to numerical errors in the
computation of `f`), when `f` returns a non-finite value (`NaN` or `Inf`),
 or when `f` has been evaluated `maxeval` times.   Note that
if the function may converge to zero, you may want
specify a nonzero `atol` (which cannot be set by default
because it depends on the scale/units of `f`); alternatively,
in such cases `extrapolate` will halt when it becomes
limited by the floating-point precision.   (Passing `breaktol=Inf`
can be useful to force `extrapolate` to continue shrinking `h` even
if polynomial extrapolation is initially failing to converge,
possibly at the cost of extraneous function evaluations.)

If `x0 = ±∞` (`±Inf`), then `extrapolate` computes the limit of
`f(x)` as `x ⟶ ±∞` using geometrically *increasing* values
of `h` (by factors of `1/contract`).

In general, the starting `h` should be large enough that `f(x0+h)`
can be computed accurately and efficiently (e.g. without
severe cancellation errors), but small enough that `f` does not
oscillate much between `x0` and `h`.  i.e. `h` should be a typical
scale over which the function `f` varies significantly.

Technically, Richardson extrapolation assumes that `f(x0+h)` can
be expanded in a power series in `h^power`, where the default
`power=1` corresponds to an ordinary Taylor series (i.e. assuming
`f` is analytic at `x0`).  If this is not true, you may obtain
slow convergence from `extrapolate`, but you can pass a different
value of `power` (e.g. `power=0.5`) if your `f` has some different
(Puiseux) power-series expansion.   Conversely, if `f` is
an *even* function around `x0`, i.e. `f(x0+h) == f(x0-h)`,
so that its Taylor series contains only *even* powers of `h`,
you can accelerate convergence by passing `power=2`.
"""
function extrapolate(f, h_::Number; contract::Number=oftype(float(real(h_)), 0.125), x0::Number=zero(h_), power::Number=1,
                     atol::Real=0, rtol::Real = atol > zero(atol) ? zero(one(float(real(x0+h_)))) : sqrt(eps(typeof(one(float(real(x0+h_)))))),
                     maxeval::Integer=typemax(Int), breaktol::Real=2)
    if isinf(x0)
        # use a change of variables x = 1/u
        return extrapolate(u -> f(inv(u)), inv(h_); rtol=rtol, atol=atol, maxeval=maxeval, contract = abs(contract) > 1 ? inv(contract) : contract, x0=inv(x0), power=power)
    end
    (rtol ≥ 0 && atol ≥ zero(atol)) || throw(ArgumentError("rtol and atol must be nonnegative"))
    breaktol > 0 || throw(ArgumentError("breaktol must be positive"))
    0 < abs(contract) < 1 || throw(ArgumentError("contract must be in (0,1)"))
    h::typeof(float(x0+h_*contract)) = h_
    invcontract = inv(contract)^power
    neville = [f(x0+h)] # the current diagonal of the Neville tableau
    f₀ = neville[1]
    err::typeof(float(norm(f₀))) = Inf
    numeval = 1
    while numeval < maxeval
        numeval += 1
        h *= contract
        push!(neville, f(x0+h))
        c = invcontract
        minerr′ = oftype(err, Inf)
        for i = length(neville)-1:-1:1
            old = neville[i]
            neville[i] = neville[i+1] + (neville[i+1] - neville[i]) / (c - 1)
            err′ = norm(neville[i] - old)
            minerr′ = min(minerr′, err′)
            if err′ < err
                f₀, err = neville[i], err′
            end
            c *= invcontract
        end
        (minerr′ > breaktol*err || !isfinite(minerr′)) && break # stop early if error increases too much
        err ≤ max(rtol*norm(f₀), atol) && break # converged
    end
    return (f₀, err)
end

# support non-numeric h as long as it is in a vector space
extrapolate(f, h; x0=zero(h), kws...) =
    extrapolate(s -> f(x0+s*h), one(norm(h)); kws...)

"""
    extrapolate(fh_itr; power=1,
                        atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf)

Similar to `extrapolate(f, h)`, performs Richardson extrapolation of a sequence of
values `f(h)` to `h → 0`, but takes an iterable collection `fh_itr` of a
sequence of `(f(h), h)` tuples (in order of decreasing `|h|`).

There is no `contract` keyword argument since the contraction factors are determined
by the sequence of `h` values (which need not contract by the same amount).  The
tolerances `atol` and `rtol` both default to `0` so that by default it examines
*all* of the values in the `fh_itr` collection.   Otherwise, the keyword arguments
have the same meanings as in `extrapolate(f, h)`.
"""
function extrapolate(fh_itr; power::Number=1, atol::Real=0, rtol::Real = 0,
                             breaktol::Real=Inf, maxeval::Integer=typemax(Int))
    (rtol ≥ 0 && atol ≥ zero(atol)) || throw(ArgumentError("rtol and atol must be nonnegative"))
    breaktol > 0 || throw(ArgumentError("breaktol must be positive"))
    Base.IteratorSize(fh_itr) isa Base.IsInfinite && iszero(atol) && iszero(rtol) && breaktol==Inf && maxeval==typemax(Int) &&
        throw(ArgumentError("stopping criteria must be supplied for infinite iterators"))
    itr = iterate(fh_itr)
    itr === nothing && throw(ArgumentError("(f,h) iterator must be non-empty"))
    (f,h), state = itr
    neville = [f] # the current diagonal of the Neville tableau
    f₀ = neville[1]
    hvals = [h]
    if Base.IteratorSize(fh_itr) isa Base.HasLength
        n = length(fh_itr)
        sizehint!(neville, n)
        sizehint!(hvals, n)
    end
    err::typeof(float(norm(f₀))) = Inf
    numeval = 1
    while numeval < maxeval
        numeval += 1
        itr = iterate(fh_itr, state)
        itr === nothing && break
        (f′,h′), state = itr
        push!(neville, f′)
        push!(hvals, h′)
        abs(h) > abs(h′) || throw(ArgumentError("|$h′| ≥ |$h| is not decreasing"))
        h = h′
        minerr′ = oftype(err, Inf)
        for i = length(neville)-1:-1:1
            c = (hvals[i] / h′)^power
            old = neville[i]
            neville[i] = neville[i+1] + (neville[i+1] - neville[i]) / (c - 1)
            err′ = norm(neville[i] - old)
            minerr′ = min(minerr′, err′)
            if err′ < err
                f₀, err = neville[i], err′
            end
        end
        (minerr′ > breaktol*err || !isfinite(minerr′)) && break # stop early if error increases too much
        err ≤ max(rtol*norm(f₀), atol) && break # converged
    end
    return (f₀, err)
end

"""
    extrapolate!(fh::AbstractVector; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf)

Similar to `extrapolate(fh)`, performs Richardson extrapolation on an array `fh`
of `(f(h), h)` tuples (in order of decreasing `|h|`), but overwrites the array
`fh` in-place with intermediate calculations.

(Thus, the array `fh` must be a vector of `Tuple{T,H}` values, where `H<:Number` is
the type of `h` and `T` is the type of the extrapolated `f(0)` **result**.  This `T`
should be a floating-point type, i.e. `fh` should contain `float(f(h))` if the
function you are extrapolating is not already floating-point-valued.)
"""
function extrapolate!(fh::AbstractVector{<:Tuple{<:Any,<:Number}}; power::Number=1, atol::Real=0, rtol::Real = 0,
                      breaktol::Real=Inf, maxeval::Integer=typemax(Int))
    (rtol ≥ 0 && atol ≥ zero(atol)) || throw(ArgumentError("rtol and atol must be nonnegative"))
    breaktol > 0 || throw(ArgumentError("breaktol must be positive"))
    isempty(fh) && throw(ArgumentError("(f,h) array must be non-empty"))
    (f₀,h) = first(fh)
    err::typeof(float(norm(f₀))) = Inf
    numeval = 1
    maxeval = min(maxeval, length(fh))
    while numeval < maxeval
        numeval += 1
        (f′,h′) = fh[numeval + (firstindex(fh)-1)]
        abs(h) > abs(h′) || throw(ArgumentError("|$h′| ≥ |$h| is not decreasing"))
        h = h′
        minerr′ = oftype(err, Inf)
        fᵢ₊₁ = f′
        for i = numeval-1:-1:1
            fᵢ,hᵢ = fh[i + (firstindex(fh)-1)]
            c = (hᵢ / h′)^power
            fᵢ₊₁ += (fᵢ₊₁ - fᵢ) / (c - 1)
            fh[i + (firstindex(fh)-1)] = (fᵢ₊₁,hᵢ)
            err′ = norm(fᵢ₊₁ - fᵢ)
            minerr′ = min(minerr′, err′)
            if err′ < err
                f₀, err = fᵢ₊₁, err′
            end
        end
        (minerr′ > breaktol*err || !isfinite(minerr′)) && break # stop early if error increases too much
        err ≤ max(rtol*norm(f₀), atol) && break # converged
    end
    return (f₀, err)
end

end # module


text:
"""
    generate_array(
        array_1::Array{Int64},
        array_2::Array{Int64},
        array_merged::Array{Int64},
        i::Int64,
        j::Int64,
        k::Int64,
        l::Int64,
        size::Int64,
        flag::Bool,
        result::Array{Any,1}
    )

`generate_array` is a modified [Breadth-First-Search]https://en.wikipedia.org/wiki/Breadth-first_search]
for generating the result-array with the possible combinations of two merged arrays with
ascending ordering.

# Arguments
- `array_1::Array{Int64}`: First sorted array
- `array_2::Array{Int64}`: Second sorted array
- `array_merged::Array{Int64}`: Merged sorted array
- `i::Int64`: Outer lower increment of first array
- `j::Int64`: Outer higher increment of first array
- `k::Int64`: Outer lower increment of second array
- `l::Int64`: Outer higher increment of second array
- `size::Int64`: Current size
- `flag::Bool`: Switching between first and second array
- `result::Array{Any,1}`: Result list of the merged arrays
"""
function generate_array(
    array_1::Array{Int64},
    array_2::Array{Int64},
    array_merged::Array{Int64},
    i::Int64,
    j::Int64,
    k::Int64,
    l::Int64,
    size::Int64,
    flag::Bool,
    result::Array{Any,1},
    )
    if flag
        if size > 1
            push!(result, array_merged[1:size])
        end
        for m in i:k
            # Define the first element
            if size < 2
                array_merged[size] = array_1[m]
                generate_array(array_1,
                array_2, array_merged, m + 1, j, k, l, size, !flag,
                result)
            else
                # Going to the third and higher element of array_1
                if array_1[m] > array_merged[size]
                    array_merged[size + 1] = array_1[m]
                    generate_array(
                        array_1,
                        array_2,
                        array_merged,
                        m + 1,
                        j,
                        k,
                        l,
                        size + 1,
                        !flag,result
                    )
                end
            end
        end

    else
        # Going to the second of array_2
        for n in j:l
            if array_2[n] > array_merged[size]
                array_merged[size + 1] = array_2[n]
                generate_array(
                    array_1,
                    array_2,
                    array_merged,
                    i,
                    n + 1,
                    k,
                    l,
                    size + 1,
                    !flag,result
                )

            end
        end

    end
    return result
end

"""
    combinations_of_2arrays(array_1::Array{Int64,1}, array_2::Array{Int64,1}))

Provides all combinations of sorted arrays with an increasing number of elements. Original
idea

# Arguments
- `array_1::Array{Int64,1}`: First sorted array
- `array_2::Array{Int64,1}`: Second sorted array

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr_1 = [10, 15, 25]
julia> arr_2 = [5, 20, 30]
julia> ClassicAlgorithmsCollections.combinations_of_2arrays(arr_1, arr_2)
[10, 20]
[10, 20, 25, 30]
[10, 30]
[15, 20]
[15, 20, 25, 30]
[15, 30]
[25, 30]
```
"""
function combinations_of_2arrays(array_1::Array{Int64,1}, array_2::Array{Int64,1})
    length_array_1,  length_array_2 = length(array_1),  length(array_2)
    init_array = zeros(Int64, (length_array_1 +  length_array_2))
     # Create array-list for the searching solution
    result = Array{Any,1}()
    return generate_array(
        array_1,
        array_2,
        init_array,
        1,
        1,
        length_array_1,  length_array_2,
        1,
        true,
        result
    )
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# ? activate Splines
import Pkg
Pkg.activate("Splines")

using PlotlyJS
using Revise
Revise.revise()

using Splines

@info "Computing"
n, m = 8, 8

# cupula
N1 = [
        cat([0:n-1 zeros(n) 4 .* sin.(range(0, stop=π, length=n))], dims=2)',
        cat([zeros(m) 0:m-1 4 .* sin.(range(0, stop=π, length=n))], dims=2)'
    ]

# angled plane
N2 = [
        cat([1:n zeros(n) zeros(n)], dims=2)',
        cat([zeros(m) 1:m 1:m], dims=2)'
     ]

# flat plane
N3 = [
        cat([0:n-1 ones(n).-1 zeros(n)], dims=2)',
        cat([ones(m).-1 0:m-1 zeros(m)], dims=2)'
     ]

# funny shape
N4 = [
        [[-1, -3, 1] [1, 3, -1] [2, 6, -3] [3, 2, -5] [4, -1, 1] [5, 1, 3] [6, 3, 7]],
        [1:7 1:7 1:7]'
     ]

fig = make_subplots(
    rows=2, cols=2,
    specs=fill(Spec(kind="scene"), 2, 2),
)
rows = [1, 1, 2, 2]
cols = [1, 2, 1, 2]

for (n, (nodes_1, nodes_2)) in enumerate([N1, N2, N3, N4])
    # compute beziers
    b1 = Bezier(nodes_1)
    b2 = Bezier(nodes_2)
    bb = b1 ⊗ b2

    # compute and plot points on the bezier curve
    points = hcat(b1(.5), b2(.5), bb(.5, .5))
    add_trace!(fig,
    plot_nodes(points, color="black", node_size=8), row=rows[n], col=cols[n]
    )

    # plot knots
    for (nodes, color) in zip((b1.nodes, b2.nodes, bb.nodes), ("red", "blue", "black"))
        add_trace!(fig,
            plot_nodes(nodes; color=color, node_size=1), row=rows[n], col=cols[n]
        )
    end

    # plot curves
    for (b, color) in zip((b1, b2), ("red", "blue"))
        add_trace!(fig,
            plot_curve(b.coordinates, color=color), row=rows[n], col=cols[n])
    end

    # plot surface
    add_trace!(
        fig,
        plot_surface(bb.coordinates, opacity=1)[1],
        row=rows[n], col=cols[n]
    )

end

display(fig)
@info "Done"

max_stars_repo_name: FedeClaudi/Splines

max_issues_repo_name: FedeClaudi/Splines

max_forks_repo_name: FedeClaudi/Splines


text:
##########################################################################################################################
###                                JULIA FUNCTIONS FOR NETWORK GENERATION                                              ###
##########################################################################################################################

## If not already installed, install required packages
#if !haskey(Pkg.installed(), "StatsBase")
#	Pkg.add("StatsBase")
#end

#= if !haskey(Pkg.installed(), "Combinatorics")
  Pkg.add("Combinatorics")
end
=#

using StatsBase

using Combinatorics

using BioSimulator

# using Combinatorics

# --------------------------------------------------- #
## FUNCTIONS FOR GENERATING QTL VALUES FOR EACH GENE ##
# --------------------------------------------------- #

# ------------------------------------------- #
## FUNCTIONS FOR SAMPLING FROM DISTRIBUTIONS ##
# ------------------------------------------- #

## Sample from an discrete exponential distribution
function sampleexpon(n, lambda, max)
  prob = (1/float(lambda))*exp.(-(1:max)/float(lambda))
  cumprob = cumsum(prob / sum(prob))
  rnb = rand(Int(n))
  res = Int64[]
  for nb in rnb
    push!(res, findfirst(x -> x >= nb, cumprob))
  end
  return res
end

## Sample from an discrete power-law distribution
function samplepowerlaw(n, gamma, max)
  prob = (1:max).^(-float(gamma))
  cumprob = cumsum(prob / sum(prob))
  rnb = rand(Int(n))
  res = Int64[]
  for nb in rnb
    push!(res, findfirst(x -> x >= nb, cumprob))
  end
  return res
end

# ------------------------------------------------------------------------------------------------ #
## FUNCTIONS FOR COMPUTING PROBA FOR PREFERENTIAL ATTACHMENT OR "INVERSE PREFERENTIAL ATTACHMENT" ##
# ------------------------------------------------------------------------------------------------ #

# Linear preferential attachment from Barabasi-Albert
function probaPA(nodes, edg)
  ki = getInDeg(nodes, edg) + 1 # in-degree of each node, add 1 so that nodes with 0 in-degree can still be picked
  res = ki/sum(ki)
  res = res/sum(res)
  return res
end

# "Inverse" preferential attachment from Lachgar
function probaInvPA(nodes, edg)
  ki = getInDeg(nodes, edg) # in-degree of each node
  if sum(ki) == 0
    res = fill(1, length(nodes))
  else
    res = 1 - ki/sum(ki)
  end
  res = res / sum(res)
  return res
end

# ------------------------------------------------------- #
## FUNCTIONS FOR GETTING THE IN- AND OUT-DEGREE OF NODES ##
# ------------------------------------------------------- #
## Edges of a graph are stored in a mx2 array (where m is the number of edges)

function getInDeg(nodes, edg)
  to = sort(edg[:,2])
  res = [length(searchsorted(to, n)) for n in nodes]
  return res
end

function getOutDeg(nodes, edg)
  from = sort(edg[:,1])
  res = [length(searchsorted(from, n)) for n in nodes]
  return res
end

# Function for checking if an edge between two nodes exist
# from is an array of nodes, to is a single nodes
function isEdge(from, to, edg)
  edgS = edg[sortperm(edg[:,2]),:]
  toedg = edgS[searchsorted(edgS[:,2], to),:]
  toedg = sort(toedg[:,1])
  [length(searchsorted(toedg, f))>0 for f in from]
end

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

## MAIN FUNCTION - GENERATE A GRAPH WITH SPECIFIED IN- AND OUT- DEGREE DISTRIBUTION

function nwgeneration(reg, target, indeg, outdeg, outdegexp, autoregproba, twonodesloop, edg = Array{Int64}(0,2)) #, targetweight = [])
## Input:
##    - reg: list of regulator nodes
##    - target: list of target nodes
##    - indeg: string variable (either "exponential" or "powerlaw") specifying the type of preferential attachment used to construct the network
##    - outdeg: string variable (either "exponential" or "powerlaw") specifying the type of distribution from which the out-degree of regulators are sampled
##    - outdegexp: the exponent of the out-degree distribution
##    - autoregproba: probability that a regulatory molecule regulates itself
##    - twonodesloop: do we allow 2-nodes loops? can be true or false
##    - optional argument edg: a 2-D array of edges already existing, to take into account; if not given creates an empty array
##    - optional argument targetweight: an array of weight for the targets, defining the probability of each target to be selected (multiplied to the probability computed from in-degree of target)

## Output:
##    - edg: A 2D array of edges, 1st column: from, 2nd column: to

  # Ensure that reg and target are arrays
  if typeof(reg) == Int64
    reg = [reg]
  end
  if typeof(target) == Int64
    target = [target]
  end

#=
  ## If a target weight vector is provided check that its length matches the number of target
  if length(targetweight)>0 & length(targetweight)!=length(target)
  	error("targetweight must match the length of target")
  end
  ## If no target weight is provided targetweight will simply be one for all targets
  if length(targetweight) == 0
  	targetweight = fill(1.0, length(target))
  end
=#

  ## Get the function for sampling from the desired out- degree distribution
  if outdeg == "exponential"
    foutdeg = getfield(current_module(), Symbol("sampleexpon"))
  elseif outdeg == "powerlaw"
    foutdeg = getfield(current_module(), Symbol("samplepowerlaw"))
  else
    error("Argument outdeg non-valid: must be exponential or powerlaw")
  end

  ## Get the function for computing the probability of each target node to receive a new incoming edge for the desired in-degree distribution
  ##  If in-degree is power law, use the model of preferential attachment from Barabasi-Albert
  ##  If in-degree is exponential, use the model of preferential attachment from Lachgar
  if indeg == "exponential"
    findeg = getfield(current_module(), Symbol("probaInvPA"))
  elseif indeg == "powerlaw"
    findeg = getfield(current_module(), Symbol("probaPA"))
  else
    error("Argument indeg non-valid: must be exponential or powerlaw")
  end

  ## Sample the number of target (out-degree) for each regulator
  out = foutdeg(length(reg), outdegexp, length(target))
  sort!(out, rev = true)

  ## Create the mx2 array of edges, 1st column = from, 2nd column = to
  # edg = Array{Int64}(0,2)

  ## For each regulator, sample from the target list its target according to its number of targets specified in the out variable
  for r in eachindex(reg)

    probTar = findeg(target, edg) # compute for each target the probability of being regulated by r
    # probTar = targetweight .* findeg(target, edg) # compute for each target the probability of being regulated by r, weighted by the weight of each target

    ## we exclude the regulator r from the list of potential targets (autoregulation is addressed later)
    if reg[r] in target
      probTar[findfirst(y -> y == reg[r], target)] = 0
    end

    ## How to deal with loops, i.e. if one or more target(s) already control the regulator r
    ## if twonodesloop == false we don't authorise the regulator to target a node that controls it
    if !twonodesloop
      exEdg = isEdge(target, reg[r], edg)
      probTar[exEdg] = 0
    end

    ## Make sure that the out-degree of regulator r doesn't exceed the number of targets with non-null proba
    out[r] = min(out[r], sum(probTar .> 0))

    ## Sample targets of regulator
    sa = StatsBase.sample(target, Weights(probTar), out[r], replace = false)

    ## Add the created edges in edg
    edg = vcat(edg, [fill(reg[r], out[r]) sa])

    ## Add an autoregulatory edge with probability autoregproba
    if rand() <= autoregproba
      edg = vcat(edg, [reg[r] reg[r]])
    end

  end

  return edg

end

function combreg(edgfrom, edgto, edgsign, p, complexsize, reac)

  edg = [edgfrom edgto edgsign]
  edg = edg[sortperm(edg[:,2]),:] ## sort edges according to the target ids
  complexes = Dict()
  complexesid = 1 ## to give a unique ID to each created complex
  rowstoremove = []
  edgtoadd = Array{Any}(0,3)
  complexsize = Int(complexsize)

  target = unique(edgto)

  for tar in target
    temp = searchsorted(edg[:, 2], tar)
    regs = edg[temp, :] ## look for the edges coming to tar

    # regspos = regs[regs[:,3] .== "1", :] ## look for the edges coming to tar representing positive regulation
    edgpos = temp[regs[:,3] .== "1"] ## stock the row number of the positive edges
    # regsneg = regs[regs[:,3] .== "-1", :] ## look for the edges coming to tar representing negative regulation
    edgneg = temp[regs[:,3] .== "-1"] ## stock the row number of the negative edges

    ntrypos = div(length(edgpos), complexsize) ## max number of complexes you can simultaneously create from the positive regulators of tar
    ntryneg = div(length(edgneg), complexsize) ## max number of complexes you can simultaneously create from the negative regulators of tar

    for i in 1:ntrypos
      if rand()<= p ## with probability p form a complex
        compo = sample(edgpos, complexsize, replace=false) ## sample the components of the complex (in fact sample rows in the edg matrix)
        edgpos = setdiff(edgpos, compo) ## remove the selected rows (regulators) from the list of possible components (future sampling)
        append!(rowstoremove, compo) ## the edges correpsponding to the selected components of the complex will be removed
        compid = string("C", reac, complexesid) ## create the new complex ID
        complexesid +=1
        complexes[compid] = edg[compo, 1] ## in the dictionary of complexes add the composition (ie array of components) of the new complex
        edgtoadd = vcat(edgtoadd, [compid tar "1"]) ## create a new regulatory edge from the complex to the target, with positive regulation
      end
    end

    ## repeat for the negative regulation
    for i in 1:ntryneg
      if rand()<= p ## with probability p form a complex
        compo = sample(edgneg, complexsize, replace=false) ## sample the components of the complex (in fact sample rows in the edg matrix)
        edgneg = setdiff(edgneg, compo) ## remove the selected rows (regulators) from the list of possible components (future sampling)
        append!(rowstoremove, compo) ## the edges correpsponding to the selected components of the complex will be removed
        compid = string("C", reac, complexesid) ## create the new complex ID
        complexesid +=1
        complexes[compid] = edg[compo, 1] ## in the dictionary of complexes add the composition (ie array of components) of the new complex
        edgtoadd = vcat(edgtoadd, [compid tar "-1"]) ## create a new regulatory edge from the complex to the target, with positive regulation
      end
    end

  end

  ## remove the individual edges from regulators forming a complex
  edg = edg[setdiff(collect(1:size(edg)[1]), rowstoremove), :]
  ## add the new edges corresponding to regulation by complexes
  edg = vcat(edg, edgtoadd)

  edg[:,1] = [string(i) for i in edg[:,1]] ## Transform the integer ID of regulators into String ID

  return Dict("newedg" => edg, "Complexes" => complexes)

end

# ------------------------------------------------------------------------------------------------ #
##              FUNCTIONS FOR CREATING THE REACTION LIST AND PROPENSITY LIST                      ##
# ------------------------------------------------------------------------------------------------ #

## Input:
##    - promList: an 1D array, each elemnt being a 2D array. Each elements correspond to the possible active states of a given promoter binding site (1st column) and the associated fold change (2nd column)
function combinactivepromstates(promList)
  proms = [t[:, 1] for t in promList]
  fcs = [t[:, 2] for t in promList]
  elsize = [length(t) for t in proms]
  combproms = Array{String}(prod(elsize), length(elsize))
  combfcs = Array{Float64}(prod(elsize), length(elsize))

  for i in eachindex(elsize)
    combproms[:, i] = repeat(proms[i], inner = prod(elsize[(i+1):end]), outer = prod(elsize[1:(i-1)]))
    combfcs[:, i] = repeat(fcs[i], inner = prod(elsize[(i+1):end]), outer = prod(elsize[1:(i-1)]))
  end

  return Dict("proms" => combproms, "fcs" => combfcs)

end

## Input:
##    - promList: an 1D array, each elemnt being an array. Each elements correspond to all possible states of the binding site
function combinallpromstates(proms)
  elsize = [length(t) for t in proms]
  combproms = Array{String}(prod(elsize), length(elsize))

  for i in eachindex(elsize)
    combproms[:, i] = repeat(proms[i], inner = prod(elsize[(i+1):end]), outer = prod(elsize[1:(i-1)]))
  end

  return combproms

end

## Gives all the possible combinations of length elem of the vector vect
function allposscomb(vect, elem)

  nb = length(vect)
  comb = Array{eltype(vect)}(nb^elem, elem)

  for i in 1:elem
    comb[:,i] = repeat(vect, inner = [nb^(elem-i)], outer = [nb^(i-1)])
  end

  return comb
end

## Generate all possible combinations of promoter states given the number of promoter and their regulatory effect (activtor or repressor)
## Input:
##    - edgSign: array with for each regulator the sign of its regulation ("1" for activation and "-1" for repression)
function allactivepromstates(edgSign)

  act = find(edgSign .== "1")
  R = length(act)
  Rtot = length(edgSign)
  comb = zeros(Int64, 2^R, Rtot)

  temp = [0, 1]
  for i in 1:R
    comb[:,act[i]] = repeat(temp, inner = [2^(R-i)], outer = [2^(i-1)])
  end

  return comb
end

function reactBioSim(reacList, prodList)
  return join(reacList, " + ") * " --> " * join(prodList, " + ")
end

function createbindingpromreactions(edg, exprstep, promPrefix, nod, functform, complexes, complexeskinetics, complexsize, complexvariants, gcnList)

  spec = []
  inicond = []
  react = []
  reactnames = []
  propens = []

  regsingl = find( x -> !in('C', x), edg["from"]) ## identify single-molecule regulators (i.e. not regulatory complexes)
  regcompl = find( x -> in('C', x), edg["from"]) ## identify regulatory complexes

  promActiveStates = Dict(string(i)*j => [] for i in nod["id"], j in gcnList) ## dictionary, 1 element for each allele version of each gene. Value is an array with m elements, 1 for each regulator.
                                                                              ## The j-th element of this array is a matrix: rows: all possible active states of the binding site of regulator j (column 1: name of the binding site state, column 2 fold-change associated with this state)
  promAllStates = Dict(string(i)*j => [] for i in nod["id"], j in gcnList) ## dictionary, 1 element for each allele version of each gene. Value is an array with m elements, 1 for each regulator.
                                                                              ## The j-th element of this array is an array, listing all the possible states of the binding site of regulator j (even if not active)

  ## Generate the binding of single regulators on binding sites of targets
  for r in regsingl, gcn in gcnList
    tarid = edg["to"][r]
    tar = string(tarid) * gcn
    reg = edg["from"][r]

    ## Create the binding site for regulator
    prom = promPrefix * tar * "reg" * reg

    ## tempactivestates is a matrix where each row correspond to an active state of the promoter:
    ##    1st column = name of species corresponding to the promoter state, 2nd column = fold change associated with the promoter state
    tempactivestates = [prom*"F" 1]
    boundactive = edg["RegSign"][r] == "1" ## are the bound states of the promoter actives i.e. able to transcribe? No if the regulator is a repressor
    tempallstates = [prom*"F"]

    ## Add the free form of the promoter to the list of species (the bound form is specific to the bound molecule)
    push!(spec, prom*"F")
    ## Add the initial abundance of the free form of the promoter to the list of initial conditions.
    ## The initial abundance of free promoter set to 1 if it is a binding site on the DNA (exprstep = "TC")
    ## If we are working on translation (exprstep = "TL"), each free binding site has an initial abundance of TCrate*qtlTCrate/(RDrate*qtlRDrate) (RNA initial abundance)
    if exprstep == "TC"
      #push!(inicond, :(1))
      push!(inicond, "1")
    else
      #push!(inicond, :(($(nod["TCrate"][tarid]/nod["RDrate"][tarid])*QTLeffects[$(gcn)]["qtlTCrate"][$(tarid)]/QTLeffects[$(gcn)]["qtlRDrate"][$(tarid)])*InitVar[$(gcn)]["R"][$(tarid)]))
        #push!(inicond, """($(nod["TCrate"][tarid]/nod["RDrate"][tarid])*QTLeffects["$(gcn)"]["qtlTCrate"][$(tarid)]/QTLeffects["$(gcn)"]["qtlRDrate"][$(tarid)])*InitVar["$(gcn)"]["R"][$(tarid)]""")
        push!(inicond, """$(nod["TCrate"][tarid]/nod["RDrate"][tarid])*InitVar["$(gcn)"]["R"][$(tarid)]""")
    end

    for gcnreg in gcnList
      prombound = prom*gcnreg*"B"
      ## add the binding reaction to the list
      push!(spec, prombound) ## add the bound promoter to the list of species
      #push!(inicond, :(0)) ## add its initial abundance to the list of initial conditions. Initial abundance of bound promoter set to 0
      push!(inicond, "0") ## add its initial abundance to the list of initial conditions. Initial abundance of bound promoter set to 0
      if boundactive
        tempactivestates = vcat(tempactivestates, [prombound edg[exprstep*"foldchange"][r]])
      end
      push!(tempallstates, prombound)
      push!(react, reactBioSim([prom*"F", functform[reg]*gcnreg], [prombound])) ## promF + reg -> promregB
      push!(reactnames, "binding"*prom*gcnreg)
      #push!(propens, :($(edg[exprstep*"bindingrate"][r])*QTLeffects[$(gcn)][$("qtl"*exprstep*"regbind")][$(tarid)]*QTLeffects[$(gcnreg)]["qtlactivity"][$(parse(Int64, reg))]))
      push!(propens, """$(edg[exprstep*"bindingrate"][r])*QTLeffects["$(gcn)"]["$("qtl"*exprstep*"regbind")"][$(tarid)]*QTLeffects["$(gcnreg)"]["qtlactivity"][$(parse(Int64, reg))]""")
      ## add the unbinding reaction to the list
      push!(react, reactBioSim([prombound], [prom*"F", functform[reg]*gcnreg])) ## promregB -> promF + reg
      push!(reactnames, "unbinding"*prom*gcnreg)
      #push!(propens, :($(edg[exprstep*"unbindingrate"][r])))
      push!(propens, """$(edg[exprstep*"unbindingrate"][r])""")
    end

    push!(promActiveStates[tar], tempactivestates)
    push!(promAllStates[tar], tempallstates)

  end

  ## Generate the binding of TL regulatory complexes on RNA-binding sites of targets
  for r in regcompl, gcn in gcnList
    tarid = edg["to"][r]
    tar = string(tarid) * gcn
    compl = edg["from"][r]

    ## Create the binding site for regulator
    prom = promPrefix * tar* "reg" * compl

    ## tempactivestates is a matrix where each lign correspond to an active state of the promoter:
    ##    1st column = name of species corresponding to the promoter state, 2nd column = fold change associated with the promoter state
    tempactivestates = [prom*"F" 1]
    boundactive = edg["RegSign"][r] == "1" ## are the bound states of the promoter actives i.e. able to transcribe? No if the regulator is a repressor
    tempallstates = [prom*"F"]

    ## Add the free and bound forms of the promoter to the list of species
    push!(spec, prom*"F")
    ## Add the initial abundance of the free form of the promoter to the list of initial conditions.
    ## The initial abundance of free promoter set to 1 if it is a binding site on the DNA (exprstep = "TC")
    ## If we are working on translation (exprstep = "TL"), each free binding site has an initial abundance of TCrate*qtlTCrate/(RDrate*qtlRDrate) (RNA initial abundance)
    if exprstep == "TC"
      #push!(inicond, :(1))
      push!(inicond, "1")
    else
      #push!(inicond, :(($(nod["TCrate"][tarid]/nod["RDrate"][tarid])*QTLeffects[$(gcn)]["qtlTCrate"][$(tarid)]/QTLeffects[$(gcn)]["qtlRDrate"][$(tarid)])*InitVar[$(gcn)]["R"][$(tarid)]))
      #push!(inicond, """($(nod["TCrate"][tarid]/nod["RDrate"][tarid])*QTLeffects["$(gcn)"]["qtlTCrate"][$(tarid)]/QTLeffects["$(gcn)"]["qtlRDrate"][$(tarid)])*InitVar["$(gcn)"]["R"][$(tarid)]""")
      push!(inicond, """$(nod["TCrate"][tarid]/nod["RDrate"][tarid])*InitVar["$(gcn)"]["R"][$(tarid)]""")
    end

    for t in 1:size(complexvariants)[1]
        complvar = compl*"comp"*join([string(complexes[compl][i])*complexvariants[t, i] for i in 1:complexsize])
        ## Create the binding of the complex on the binding site
        prombound = prom*join(complexvariants[t, :])*"B"
        push!(spec, prombound)
        #push!(inicond, :(0)) ## add its initial abundance to the list of initial conditions. Initial abundance of bound promoter set to 0
        push!(inicond, "0") ## add its initial abundance to the list of initial conditions. Initial abundance of bound promoter set to 0
        if boundactive
            tempactivestates = vcat(tempactivestates, [prombound edg[exprstep*"foldchange"][r]])
        end
        push!(tempallstates, prombound)
        push!(react, reactBioSim([prom*"F", complvar], [prombound])) ## promF + complvar -> promB
        push!(reactnames, "binding"*prom*complvar)
        prop = """$(edg[exprstep*"bindingrate"][r])*QTLeffects["$(gcn)"]["$("qtl"*exprstep*"regbind")"][$(tarid)]"""*join(["""*QTLeffects["$(complexvariants[t, i])"]["qtlactivity"][$(complexes[compl][i])]""" for i in 1:complexsize])
        #push!(propens, parse(prop))
        push!(propens, prop)
        ## add the unbinding of the complex from the binding site
        push!(react, reactBioSim([prombound], [prom*"F", complvar])) ## promB -> promF + complvar
        push!(reactnames, "unbinding"*prom*complvar)
        #push!(propens, :($(edg[exprstep*"unbindingrate"][r])))
        push!(propens, """$(edg[exprstep*"unbindingrate"][r])""")
    end

    push!(promActiveStates[tar], tempactivestates)
    push!(promAllStates[tar], tempallstates)

  end

  return Dict("species" => spec, "initialconditions" => inicond, "reactions" => react, "reactionsnames" => reactnames, "propensities" => propens, "promActiveStates" => promActiveStates, "promAllStates" => promAllStates)

end

function generateReactionList(nod, edgTCRN, edgTLRN, edgRDRN, edgPDRN, edgPTMRN, complexes, complexeskinetics, complexsize, gcnList)

    ## ensures that gcnList is a vector (not true if there is only one allele for each gene)
    if typeof(gcnList) == String
      gcnList = [gcnList]
    end

    ## Output of the function
    species = []
    initialconditions = []
    reactions = []
    reactionsnames = []
    propensities = []

    ## Specify the active form of each regulator, ie the species that performs the regulation
    functform = Dict(zip(map(string, nod["id"]), nod["ActiveForm"]))

    ## Create all possible allele combinations of the components of a complex
    complexvariants = allposscomb(gcnList, complexsize)

    ## Generate the formation and dissociation reactions of the different regulatory complexes
    for compl in keys(complexes), t in 1:size(complexvariants)[1]
        complvar = compl*"comp"*join([string(complexes[compl][i])*complexvariants[t, i] for i in 1:complexsize])
        push!(species, complvar) ## add the complex form to the list of species
        #push!(initialconditions, :(0)) ## add the initial abundance of the complex form to the list of initial conditions. For complexes, initial abundance set to 0
        push!(initialconditions, "0") ## add the initial abundance of the complex form to the list of initial conditions. For complexes, initial abundance set to 0
        ## Create the reaction of complex formation
        push!(reactions, reactBioSim([functform[string(complexes[compl][i])]*complexvariants[t, i] for i in 1:complexsize], [complvar])) ## sum of complex components -> compl
        push!(reactionsnames, "formation"*complvar)
        #push!(propensities, :($(complexeskinetics[compl]["formationrate"])))
        push!(propensities, """$(complexeskinetics[compl]["formationrate"])""")
        ## Create the reaction of complex dissociation
        push!(reactions, reactBioSim([complvar], [functform[string(complexes[compl][i])]*complexvariants[t, i] for i in 1:complexsize])) ## compl -> sum of complex components
        push!(reactionsnames, "dissociation"*complvar)
        #push!(propensities, :($(complexeskinetics[compl]["dissociationrate"])))
        push!(propensities, """$(complexeskinetics[compl]["dissociationrate"])""")
    end

    TLbinding = createbindingpromreactions(edgTLRN, "TL", "RBS", nod, functform, complexes, complexeskinetics, complexsize, complexvariants, gcnList)
    promactiveTL = TLbinding["promActiveStates"]
    promallTL = TLbinding["promAllStates"]

    TCbinding = createbindingpromreactions(edgTCRN, "TC", "Pr", nod, functform, complexes, complexeskinetics, complexsize, complexvariants, gcnList)
    promactiveTC = TCbinding["promActiveStates"]
    promallTC = TCbinding["promAllStates"]

    species  = vcat(species, TCbinding["species"], TLbinding["species"])
    initialconditions  = vcat(initialconditions, TCbinding["initialconditions"], TLbinding["initialconditions"])
    reactions  = vcat(reactions, TCbinding["reactions"], TLbinding["reactions"])
    reactionsnames  = vcat(reactionsnames, TCbinding["reactionsnames"], TLbinding["reactionsnames"])
    propensities  = vcat(propensities, TCbinding["propensities"], TLbinding["propensities"])

    ## CREATE TRANSCRIPTION AND RNA DECAY REACTIONS

    for g in nod["id"], gcn in gcnList
        gname = string(g) * gcn
        TCreg = promactiveTC[gname]
        TLreg = promallTL[gname]

        ## What is the RNA form of the gene?
        ##    R[gene] if the gene is not controlled at the translation level
        ##    otherwise transcription produces the free (unbound) form of each binding site present on the RNA (1 site per TL regulator)
        if length(TLreg) == 0
            RNAform = ["R"*gname]
            push!(species, "R"*gname)
            #push!(initialconditions, :(($(nod["TCrate"][g]/nod["RDrate"][g])*QTLeffects[$(gcn)]["qtlTCrate"][$(g)]/QTLeffects[$(gcn)]["qtlRDrate"][$(g)])*InitVar[$(gcn)]["R"][$(g)]))
            #push!(initialconditions, """($(nod["TCrate"][g]/nod["RDrate"][g])*QTLeffects["$(gcn)"]["qtlTCrate"][$(g)]/QTLeffects["$(gcn)"]["qtlRDrate"][$(g)])*InitVar["$(gcn)"]["R"][$(g)]""")
            push!(initialconditions, """$(nod["TCrate"][g]/nod["RDrate"][g])*InitVar["$(gcn)"]["R"][$(g)]""")

            transcriptforms = ["R"*gname] ## for RNA decay
        else
            RNAform = [t[1] for t in TLreg] ## The TLreg dictionary is constructed such that each free binding site name is at the position 1
            transcriptforms = combinallpromstates(TLreg) ## for RNA decay
        end

        if length(TCreg) == 0
            ## Create transcription of the gene
            push!(reactions, reactBioSim([0], RNAform))
            push!(reactionsnames, "transcription"*gname)
            #push!(propensities, :($(nod["TCrate"][g])*QTLeffects[$(gcn)]["qtlTCrate"][$(g)]))
            push!(propensities, """$(nod["TCrate"][g])*QTLeffects["$(gcn)"]["qtlTCrate"][$(g)]""")
        else
            promcomb = combinactivepromstates(TCreg) ## gives all possible active combinations of the different binding site states
            for t in 1:size(promcomb["proms"])[1]
                push!(reactions, reactBioSim(promcomb["proms"][t,:], vcat(promcomb["proms"][t,:], RNAform)))
                push!(reactionsnames, "transcription"*join(promcomb["proms"][t,:]))
                #push!(propensities, :($(nod["TCrate"][g])*QTLeffects[$(gcn)]["qtlTCrate"][$(g)]*$(prod(promcomb["fcs"][t,:]))))
                push!(propensities, """$(nod["TCrate"][g])*QTLeffects["$(gcn)"]["qtlTCrate"][$(g)]*$(prod(promcomb["fcs"][t,:]))""")
            end
        end

        ## Generate RNA decay reactions
        edgreg = find(edgRDRN["to"] .== g)
        edgsingl = edgreg[map(x -> !in('C',x), edgRDRN["from"][edgreg])]
        edgcompl = edgreg[map(x -> in('C',x), edgRDRN["from"][edgreg])]

        for t in 1:size(transcriptforms)[1]
            ## Basal decay rate
            push!(reactions, reactBioSim(transcriptforms[t,:], [0]))
            push!(reactionsnames, "RNAdecay"*join(transcriptforms[t,:]))
            #push!(propensities, :($(nod["RDrate"][g])*QTLeffects[$(gcn)]["qtlRDrate"][$(g)]))
            push!(propensities, """$(nod["RDrate"][g])*QTLeffects["$(gcn)"]["qtlRDrate"][$(g)]""")

            for r in edgsingl, gcnreg in gcnList
                reg = edgRDRN["from"][r]
                push!(reactions, reactBioSim(vcat(transcriptforms[t,:], functform[reg]*gcnreg), [functform[reg]*gcnreg]))
                push!(reactionsnames, "RNAdecay"*join(transcriptforms[t,:])*"reg"*reg*gcnreg)
                #push!(propensities, :($(edgRDRN["RDbindingrate"][r])*QTLeffects[$(gcnreg)]["qtlactivity"][$(parse(Int64, reg))]))
                push!(propensities, """$(edgRDRN["RDbindingrate"][r])*QTLeffects["$(gcnreg)"]["qtlactivity"][$(parse(Int64, reg))]""")
            end

            for r in edgcompl, j in 1:size(complexvariants)[1]
                compl = edgRDRN["from"][r]
                complvar = compl*"comp"*join([string(complexes[compl][i])*complexvariants[j, i] for i in 1:complexsize])
                push!(reactions, reactBioSim(vcat(transcriptforms[t,:], complvar), [complvar]))
                push!(reactionsnames, "RNAdecay"*join(transcriptforms[t,:])*"reg"*complvar)
                propens = """$(edgRDRN["RDbindingrate"][r])"""*join(["*"*"""QTLeffects["$(complexvariants[j, i])"]["qtlactivity"][$(complexes[compl][i])]""" for i in 1:complexsize])
                #push!(propensities, parse(propens))
                push!(propensities, propens)
            end
        end
    end

    ## CREATE TRANSLATION AND PROTEIN DECAY REACTIONS
    for g in nod["id"][nod["coding"] .== "PC"], gcn in gcnList
        gname = string(g) * gcn
        TLreg = promactiveTL[gname]

        push!(species, "P"*gname)
        ## Initial abundance of a protein = TCrate/RDrate * TLrate/PDrate (given the QTL effects as well)
        #push!(initialconditions, :(($(nod["TCrate"][g]*nod["TLrate"][g]/(nod["RDrate"][g]*nod["PDrate"][g]))*QTLeffects[$(gcn)]["qtlTCrate"][$(g)]*QTLeffects[$(gcn)]["qtlTLrate"][$(g)]/(QTLeffects[$(gcn)]["qtlRDrate"][$(g)]*QTLeffects[$(gcn)]["qtlPDrate"][$(g)]))*InitVar[$(gcn)]["P"][$(g)]))
        #push!(initialconditions, """($(nod["TCrate"][g]*nod["TLrate"][g]/(nod["RDrate"][g]*nod["PDrate"][g]))*QTLeffects["$(gcn)"]["qtlTCrate"][$(g)]*QTLeffects["$(gcn)"]["qtlTLrate"][$(g)]/(QTLeffects["$(gcn)"]["qtlRDrate"][$(g)]*QTLeffects["$(gcn)"]["qtlPDrate"][$(g)]))*InitVar["$(gcn)"]["P"][$(g)]""")
        push!(initialconditions, """($(nod["TCrate"][g]*nod["TLrate"][g]/(nod["RDrate"][g]*nod["PDrate"][g])))*InitVar["$(gcn)"]["P"][$(g)]""")

        if length(TLreg) == 0
            ## Create transltion of the gene
            push!(reactions, reactBioSim(["R"*gname], ["R"*gname, "P"*gname]))
            push!(reactionsnames, "translation"*gname)
            #push!(propensities, :($(nod["TLrate"][g])*QTLeffects[$(gcn)]["qtlTLrate"][$(g)]))
            push!(propensities, """$(nod["TLrate"][g])*QTLeffects["$(gcn)"]["qtlTLrate"][$(g)]""")

        else
            promcomb = combinactivepromstates(TLreg) ## gives all possible active combinations of the different binding site states
            for t in 1:size(promcomb["proms"])[1]
                push!(reactions, reactBioSim(promcomb["proms"][t,:], vcat(promcomb["proms"][t,:], "P"*gname)))
                push!(reactionsnames, "translation"*join(promcomb["proms"][t,:]))
                #push!(propensities, :($(nod["TLrate"][g])*QTLeffects[$(gcn)]["qtlTLrate"][$(g)]*$(prod(promcomb["fcs"][t,:]))))
                push!(propensities, """$(nod["TLrate"][g])*QTLeffects["$(gcn)"]["qtlTLrate"][$(g)]*$(prod(promcomb["fcs"][t,:]))""")
            end
        end

        ## Create basal protein decay rate
        push!(reactions, reactBioSim(["P"*gname], [0]))
        push!(reactionsnames, "proteindecay"*"P"*gname)
        #push!(propensities, :($(nod["PDrate"][g])*QTLeffects[$(gcn)]["qtlPDrate"][$(g)]))
        push!(propensities, """$(nod["PDrate"][g])*QTLeffects["$(gcn)"]["qtlPDrate"][$(g)]""")
        if nod["PTMform"][g] =="1"

            ## Add to the list of species the PTM form of the protein
            push!(species, "Pm"*string(g)*gcn)
            #push!(initialconditions, :(0)) ## initial abundance of modified proteins set to 0
            push!(initialconditions, "0") ## initial abundance of modified proteins set to 0

            push!(reactions, reactBioSim(["Pm"*gname], [0]))
            push!(reactionsnames, "proteindecay"*"Pm"*gname)
            #push!(propensities, :($(nod["PDrate"][g])*QTLeffects[$(gcn)]["qtlPDrate"][$(g)]))
            push!(propensities, """$(nod["PDrate"][g])*QTLeffects["$(gcn)"]["qtlPDrate"][$(g)]""")
        end
    end

      ## CREATE REGULATED PROTEIN DECAY

    regsingl = find( x -> !in('C', x), edgPDRN["from"]) ## identify single-molecule regulators (i.e. not regulatory complexes)
    regcompl = find( x -> in('C', x), edgPDRN["from"]) ## identify regulatory complexes

    for r in regsingl, gcn in gcnList
        tarid = edgPDRN["to"][r]
        tar = string(tarid) * gcn
        reg = edgPDRN["from"][r]

        pform = ["P"*tar]
        if nod["PTMform"][tarid] == "1"
          push!(pform, "Pm"*tar)
        end

        for p in pform, gcnreg in gcnList
            push!(reactions, reactBioSim([p, functform[reg]*gcnreg], [functform[reg]*gcnreg]))
            push!(reactionsnames, "proteindecay"*p*"reg"*reg*gcnreg)
            #push!(propensities, :($(edgPDRN["PDbindingrate"][r])*QTLeffects[$(gcnreg)]["qtlactivity"][$(parse(Int64, reg))]))
            push!(propensities, """$(edgPDRN["PDbindingrate"][r])*QTLeffects["$(gcnreg)"]["qtlactivity"][$(reg)]""")
        end
    end

    for r in regcompl, gcn in gcnList
        tarid = edgPDRN["to"][r]
        tar = string(tarid) * gcn
        compl = edgPDRN["from"][r]

        pform = ["P"*tar]
        if nod["PTMform"][tarid] == "1"
            push!(pform, "Pm"*tar)
        end

        for p in pform, t in 1:size(complexvariants)[1]
            complvariant = compl*"comp"*join([string(complexes[compl][i])*complexvariants[t, i] for i in 1:complexsize])
            push!(reactions, reactBioSim([p, complvariant], [complvariant]))
            push!(reactionsnames, "proteindecay"*p*"reg"*complvariant)
            propens = """$(edgPDRN["PDbindingrate"][r])"""*join(["*"*"""QTLeffects["$(complexvariants[t, i])"]["qtlactivity"][$(complexes[compl][i])]""" for i in 1:complexsize])
            #push!(propensities, parse(propens))
            push!(propensities, propens)
        end
    end

  ## CREATE PROTEIN POST-TRANSLATIONAL MODIFICATION

    regsingl = find( x -> !in('C', x), edgPTMRN["from"]) ## identify single-molecule regulators (i.e. not regulatory complexes)
    regcompl = find( x -> in('C', x), edgPTMRN["from"]) ## identify regulatory complexes

    for r in regsingl, gcn in gcnList
        tarid = edgPTMRN["to"][r]
        tar = string(tarid) * gcn
        reg = edgPTMRN["from"][r]
        ispos = edgPTMRN["RegSign"][r] == "1" ## is the regulator transforming the orginal protein into its modified form (RegSign = "1") or the opposite (RegSign = "-1")

        for gcnreg in gcnList
            if ispos
                push!(reactions, reactBioSim([ "P"*tar, functform[reg]*gcnreg], ["Pm"*tar, functform[reg]*gcnreg]))
                push!(reactionsnames, "PTM"*tar*"reg"*reg*gcnreg)
                #push!(propensities, :($(edgPTMRN["PTMbindingrate"][r])*QTLeffects[$(gcnreg)]["qtlactivity"][$(parse(Int64, reg))]))
                push!(propensities, """$(edgPTMRN["PTMbindingrate"][r])*QTLeffects["$(gcnreg)"]["qtlactivity"][$(reg)]""")
            else
                push!(reactions, reactBioSim([ "Pm"*tar, functform[reg]*gcnreg], ["P"*tar, functform[reg]*gcnreg]))
                push!(reactionsnames, "de-PTM"*tar*"reg"*reg*gcnreg)
                #push!(propensities, :($(edgPTMRN["PTMbindingrate"][r])*QTLeffects[$(gcnreg)]["qtlactivity"][$(parse(Int64, reg))]))
                push!(propensities, """$(edgPTMRN["PTMbindingrate"][r])*QTLeffects["$(gcnreg)"]["qtlactivity"][$(reg)]""")
            end
        end
    end

    for r in regcompl, gcn in gcnList
        tarid = edgPTMRN["to"][r]
        tar = string(tarid) * gcn
        compl = edgPTMRN["from"][r]
        ispos = edgPTMRN["RegSign"][r] == "1" ## is the regulator transforming the orginal protein into its modified form (RegSign = "1") or the opposite (RegSign = "-1")

        for t in 1:size(complexvariants)[1]
            complvariant = compl*"comp"*join([string(complexes[compl][i])*complexvariants[t, i] for i in 1:complexsize])

            if ispos
                push!(reactions, reactBioSim([ "P"*tar, complvariant], ["Pm"*tar, complvariant]))
                push!(reactionsnames, "PTM"*tar*"reg"*complvariant)
                propens = """$(edgPTMRN["PTMbindingrate"][r])"""*join(["*"*"""QTLeffects["$(complexvariants[t, i])"]["qtlactivity"][$(complexes[compl][i])]""" for i in 1:complexsize])
                #push!(propensities, parse(propens))
                push!(propensities, propens)
            else
                push!(reactions, reactBioSim([ "Pm"*tar, complvariant], ["P"*tar, complvariant]))
                push!(reactionsnames, "de-PTM"*tar*"reg"*complvariant)
                propens = """$(edgPTMRN["PTMbindingrate"][r])"""*join(["*"*"""QTLeffects["$(complexvariants[t, i])"]["qtlactivity"][$(complexes[compl][i])]""" for i in 1:complexsize])
                #push!(propensities, parse(propens))
                push!(propensities,propens)
            end
        end
    end

    return Dict("species" => species, "initialconditions" => initialconditions, "reactions" => reactions, "reactionsnames" => reactionsnames, "propensities" => propensities, "TCproms" => promallTC, "TLproms" => promallTL)
end

function getDictfromKey(mydict, mykey)
    return mydict[mykey]
end

# ------------------------------------------------------------------------------------------------ #
##              FUNCTIONS FOR CREATING THE REACTION LIST AND PROPENSITY LIST                      ##
# ------------------------------------------------------------------------------------------------ #

function res2df(output :: BioSimulator.SimData)
    t, data = get_data(output)
    id2ind  = output.id2ind

    d1, d2, d3 = size(data)

    ids = collect(keys(id2ind))

    df = DataFrame()
    df[:time] = repeat(collect(t), outer=[d3])
    for i = ids
        df[i] = vec(convert(Array, view(data, id2ind[i], :, :)))
    end
    df[:trial] = repeat(collect(1:d3), inner=[d2])
    return df
end

#=
function smallmodel()

    model = BioSimulator.Network("coucou")

    model <= BioSimulator.Species("X", 25)
    model <= BioSimulator.Reaction("transcription", 0.02, "0 --> X")
    model <= BioSimulator.Reaction("decay", 0.00069, "X --> 0")

    result = BioSimulator.simulate(model, algorithm = SSA, time = 1000.0, epochs = 1000)
    resdf = res2df(result)

    return resdf
end
=#

function allequal(x)
    all(y->y==x[1], x)
end

function stochasticsimulation(stochmodel, QTLeffects, InitVar, nod, simtime; modelname = "MySimulation", ntrials = 1, nepochs = -1, simalgorithm = "SSA")

    try
        ## If no nepochs provided, record abundance at each time units of the simulation
        if nepochs == -1
            nepochs = simtime
        end

        ## Convert the String name of the simulator to use into a BioSimulator object of type inheriting from Algorithm
        if in(simalgorithm, ["SSA", "FRM", "NRM", "ODM", "SAL"])
            simalgorithm = eval(parse("BioSimulator."*simalgorithm))
        else
            error("Specified algorithm is not implemented in module BioSimulator")
        end

        model = BioSimulator.Network(modelname)

        ## Add the species in the model, with their initial abundance
        for i in 1:length(stochmodel["species"])
            i0 = replace(stochmodel["initialconditions"][i], "QTLeffects", "$QTLeffects")
            i0 = replace(i0, "InitVar", "$InitVar")
            i0 = eval(parse(i0))
            #println(stochmodel["species"][i]* "\t"*string(i0))
            model <= BioSimulator.Species(stochmodel["species"][i], round(Int, i0))

            if !isa(i0, Number)
                println(stochmodel["initialconditions"][i])
                error("Pb i0")
            end
            if typeof(stochmodel["species"][i]) != String
                println(stochmodel["species"][i])
                error("Pb species")
            end
        end

        ## Add the reactions in the model, with their name and rate
        for i in eachindex(stochmodel["reactions"])
            prop = replace(stochmodel["propensities"][i], "QTLeffects", "$QTLeffects")
            #println(prop)
            prop = eval(parse(prop))
            model <= BioSimulator.Reaction(stochmodel["reactionsnames"][i], prop, stochmodel["reactions"][i])

            if !isa(prop, Number)
                println(stochmodel["propensities"][i])
                error("Pb prop")
            end
            if typeof(stochmodel["reactionsnames"][i]) != String
                println(stochmodel["reactionsnames"][i])
                error("Pb reactionsnames")
            end
            if typeof(stochmodel["reactions"][i]) != String
                println(stochmodel["reactions"][i])
                error("Pb reactions")
            end

        end

#=
        for i in collect(keys(model.reaction_list))
          print(model.reaction_list[i])
          print("\t")
          print(model.reaction_list[i].rate, "\n")
        end

=#

        #println("JULIA> Running simulation ...")
        #tic()
        #JSON.print(getfield(Main, :RJuliaSocket), toR("please wait"))
        result = simulate(model, algorithm = simalgorithm, time = convert(Float64, simtime), epochs = round(Int64, nepochs), trials = convert(Int64, ntrials))
        #toc()
        #println("JULIA> Done.")

        resultdf = res2df(result)

        abundancedf = resultdf[:, [:time, :trial]]

        for g in collect(keys(stochmodel["TCproms"]))

            gid = parse(Int64, replace(g, r"GCN.+$","")) ## gives the gene id

            ## Check that for each binding site on the promoter of each gene, at each time the sum of the abundance of all species corresponding to all possible binding site states equals 1 (bc only 1 binding site per gene)
            for i in eachindex(stochmodel["TCproms"][g])
                prabundance = resultdf[:, map(Symbol, stochmodel["TCproms"][g][i])]
                sumprom = [sum(convert(Array, row)) for row in eachrow(prabundance)]
                if any(sumprom .!=1)
                    error("The sum of promoter states for "*stochmodel["TCproms"][g][i][1]*"not equal to 1.")
                end
            end

            if length(stochmodel["TLproms"][g]) > 0

                ## Check for each RNA that the different binding sites on the RNA are in equal abundance at each time of the simulation
                rbsabundance = [sum(convert(Array, resultdf[t, map(Symbol, x)]))for t in 1:size(resultdf)[1], x in stochmodel["TLproms"][g]]

                if !all(mapslices(allequal, rbsabundance, 2))
                    error("The abundance of the different binding sites on the RNA "*g*"are not equal.")
                end

                ## Add to abundancedf a column corresponding to the abundance of the RNA associated with g
                abundancedf[Symbol("R"*g)] = rbsabundance[:,1]
            else
                abundancedf[Symbol("R"*g)] = resultdf[:, Symbol("R"*g)]
            end

            ## MAYBE to change if we don't make the disctinction between original and modified protein
            if nod["coding"][gid] == "PC"
                abundancedf[Symbol("P"*g)] = resultdf[:, Symbol("P"*g)]

                if nod["PTMform"][gid] == "1"
                    abundancedf[Symbol("Pm"*g)] = resultdf[:, Symbol("Pm"*g)]
                end

            end
        end

        return abundancedf

    catch err
        isa(err, InterruptException) || rethrow(err)
        return nothing
    end

end

#=
function runstochsim(stochmodel, QTLeffects, InitVar, nod, simtime; modelname = "MySimulation", ntrials = 1, nepochs = -1, simalgorithm = "SSA")
    try
        res = stochasticsimulation(stochmodel, QTLeffects, InitVar, nod, simtime, modelname = modelname, ntrials = ntrials, nepochs = nepochs, simalgorithm = simalgorithm)
        return res
    catch err
        isa(err, InterruptException) || rethrow(err)
        return nothing
    end
end
=#

function whatisit(x)
  println(x)
  println(typeof(x))
end

# ------------------------------------------------------------------------------------------------ #
##              SOLUTION TO RUN A SIMULATION AND STOP IT IF RUNNING TIME IS                       ##
##                         LARGER THAN A DEFINED TIME LIMIT                                       ##
##                    CURRENTLY ONLY WORKS FROM THE JULIA TERMINAL                                ##
# ------------------------------------------------------------------------------------------------ #

##   !!!WARNING!!! DO NOT UNCOMMENT IN THE julia_functions.jl FILE OR IT WILL BUG

#=
p = addprocs(1)[1]
@everywhere include("winData/multiomics_networks_simulation/julia_functions.jl")

@everywhere function stochsimtimelimit(p, stochmodel, QTLeffects, InitVar, nod, simtime; modelname = "MySimulation", ntrials = 1, nepochs = -1, simalgorithm = "SSA", time_limit = Inf)

    output = Channel(1) ## create a Channel, that will store the simulation results
    @async put!(output, remotecall_fetch(stochasticsimulation, p, stochmodel, QTLeffects, InitVar, nod, 1)) ## start the simulation on the process 2 of the Julia evaluator

    start=time() ## start the timer

    while !isready(output) & (time() - start < time_limit) ## While the channel is not ready (i.e. computation not done) and the timer is < to the defined value (default no time limit)
      sleep(0.1)
    end

    if !isready(output) ## As soon as the Channel is ready or the time is out, check if the computation is done
      interrupt(2) ## if not interrupt the second process i.e. the computation
      data = "timeout" ## return a timeout message
      #close(output)
    else
      data = take!(output) ## if the computation is done take the result from the Channel
    end

    return data
end

## Loading an example
using JLD
gentil = load("/home/oangelin/Documents/goodmodel.jld")
QTLeffects = gentil["QTLeffects"]
InitVar = gentil["InitVar"]
stochmodel = gentil["stochmodel"]
nod = load("/home/oangelin/Documents/goodmodelnod.jld")["nod"]

## Testing on the example

@time stochsimtimelimit(p, stochmodel, QTLeffects, InitVar, nod, 1, time_limit = 10)
@time stochsimtimelimit(p, stochmodel, QTLeffects, InitVar, nod, 1, time_limit = 70)
@time stochsimtimelimit(p, stochmodel, QTLeffects, InitVar, nod, 1)
=#


text:
"""
    A = JopHighpass(sp)

Build a 2D Highpass operator, which first apply a center weighted three point
smoothing operator [0.25, 0.5, 0.25] to each dimension with a given iterations
number (nit[1:2]) to get the low frequency component of the data, then subtract
it from the original data to return the high pass component.  It works with
`sp::JetSpace`, `nit::Array{Integer, 1}` and `A::JopHighpass`.  It can specify
different number of iterations (nit) in x and z direction. For a seismic image,
more iterations in z direction is often applied to obtain more vertical
resolution.

For example,
```julia
A = JopHighpass(JetSpace(Float64,128,256), nit=(3,1)])
m = rand(domain(A))
d = A*m
```
"""
JopHighpass(spc::JetSpace; nit=(1,1)) = JopLn(dom = spc, rng = spc, df! = JopHighpass_df!, s = (nit=nit,))

export JopHighpass

function JopHighpass_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; nit, kwarg...) where {T}
    if all(nit .== 0)
        copyto!(d,m)
        return d
    end

    nz,nx = size(m)

    temp1 = zeros(T, nz+2, nx)
    temp2 = zeros(T, nz, nx)
    padding_z!(temp1, m)
    for it = 1:nit[1]
        for ix = 1:nx
            for iz = 1:nz
                temp2[iz, ix] = 0.25*(temp1[iz,ix] + temp1[iz+2,ix]) + 0.5 * temp1[iz+1,ix]
            end
        end
        padding_z!(temp1, temp2)
    end
    temp3 = zeros(T,nz,nx+2)
    temp4 = zeros(T,nz,nx)
    padding_x!(temp3, temp1[2:nz+1,:])
    for it = 1:nit[2]
        for ix = 1:nx
            for iz = 1:nz
                temp4[iz, ix] = 0.25*(temp3[iz,ix] + temp3[iz,ix+2]) + 0.5 * temp3[iz,ix+1]
            end
        end
        padding_x!(temp3, temp4)
    end
    for ix = 1:nx
        for iz = 1:nz
            d[iz,ix] = m[iz,ix] - temp3[iz,ix+1]
        end
    end
    d
end

function padding_x!(dout::Array{T,2}, din::Array{T,2}) where T<:AbstractFloat
    dout[:,2:end-1] = din
    dout[:,1] = dout[:,2]
    dout[:,end] = dout[:,end-1]
end

function padding_z!(dout::Array{T,2}, din::Array{T,2}) where T<:AbstractFloat
    dout[2:end-1,:] = din
    dout[1,:] = dout[2,:]
    dout[end,:] = dout[end-1,:]
end


text:
using AbstractPlotting

using Test

include("quaternions.jl")

# TODO write some AbstractPlotting specific tests... So far functionality is tested in Makie.jl


text:
export MPSSurveillanceUnit

struct MPSSurveillanceUnit <: SurveillanceUnit
  sd::SerialDevice
  CommandStart::String
  CommandEnd::String
  delim::String
end

neededDependencies(::MPSSurveillanceUnit) = []
optionalDependencies(::MPSSurveillanceUnit) = []

Base.close(su::MPSSurveillanceUnit) = close(su.sd.sp)

function MPSSurveillanceUnit(params::Dict)
  # Here we could put more parameters into the TOML file
  su = MPSSurveillanceUnit(params["connection"])
 # DisableWatchDog(su)
  #DisableSurveillance(su)
  return su
end

function ArduinoCommand(Arduino::MPSSurveillanceUnit, cmd::String)
    cmd=Arduino.CommandStart*cmd*Arduino.CommandEnd*Arduino.delim;
   return query(Arduino.sd,cmd);
end

function MPSSurveillanceUnit(portAdress::AbstractString)
    pause_ms::Int=30
    timeout_ms::Int=500
    delim::String="#"
    delim_read::String="#"
    delim_write::String="#"
    baudrate::Integer = 115200
    CommandStart="!"
    CommandEnd="*"
    ndatabits::Integer=8
    parity::SPParity=SP_PARITY_NONE
    nstopbits::Integer=1
    sp = SerialPort(resolvedSymlink(portAdress))
    open(sp)
	set_speed(sp, baudrate)
	set_frame(sp,ndatabits=ndatabits,parity=parity,nstopbits=nstopbits)
	#set_flow_control(sp,rts=rts,cts=cts,dtr=dtr,dsr=dsr,xonxoff=xonxoff)
    sleep(2)
    flush(sp)
    write(sp, "!ConnectionEstablished*#")
    response = readuntil(sp, Vector{Char}(delim_read), timeout_ms);
    @show response
    if(response == "MPSSurveillanceV1")
            @info "Connection to MPSSurveillanceUnit established"
            return MPSSurveillanceUnit(SerialDevice(sp,pause_ms, timeout_ms, delim_read, delim_write),CommandStart,CommandEnd,delim)
    else
            @warn "Connected to wrong Device" response
            return nothing;
    end
end

function getTemperatures(Arduino::MPSSurveillanceUnit)
    Temps = ArduinoCommand(Arduino, "GET:TEMP")
    TempDelim = "T"

    temp =  tryparse.(Float64,split(Temps,TempDelim))

    # We hardcode this to four here
    tempFloat = zeros(2)

    for i=1:min(length(tempFloat),length(temp))
      if temp[i] != nothing
          tempFloat[i] = temp[i]
      end
    end

    return tempFloat
end

function enablePMPS(scanner::MPIScanner, rp::RedPitaya)
  DIO(rp,scanner.params["DAQ"]["pinRSDSM1850"],true) #activate SM18-50
  DIO(rp,scanner.params["DAQ"]["pinInterlockHubert"],true) #activate Hubert
end

function disablePMPS(scanner::MPIScanner, rp::RedPitaya)
  DIO(rp,scanner.params["DAQ"]["pinRSDSM1850"],false) #remoteshutdown SM18-50
  DIO(rp,scanner.params["DAQ"]["pinInterlockHubert"],false) #interlock Hubert
end

function enableACPower(su::MPSSurveillanceUnit, scanner::MPIScanner)
  @info "Enable AC Power"
  daq = getDAQ(scanner)
  rp = master(daq.rpc)

  temps = getTemperatures(su)
  tempMax = scanner.params["SurveillanceUnit"]["maxTemperatures"]

  #code for pMPS
  if scanner.params["General"]["scannerName"] == "pMPS"
    if iszero(temps .> tempMax) != false
      @info "pMPS: Hubert interlock and RSD off. Ready to measure."
      enablePMPS(scanner, rp)
    else
      @warn "pMPS: maxTemp exceeded. Hubert interlock and RSD stay on." temps tempMax
      disablePMPS(scanner, rp)
    end
  end
end

function disableACPower(su::MPSSurveillanceUnit, scanner::MPIScanner)
  @info "Disable AC Power"
  daq = getDAQ(scanner)
  rp = master(daq.rpc)

  #code for pMPS
  if scanner.params["General"]["scannerName"] == "pMPS"
      @info "pMPS: Hubert interlock and RSD on."
      disablePMPS(scanner, rp)
  end
end


text:
"""
    build_data(; kwargs...)
"""
function build_data(; save = true, overwrite = false)
    set = read_set(joinpath("src","readfiles","setlist.yml"))

    disagg = read_build("disagg"; save = save, overwrite = overwrite);
    !isempty(disagg) && (return (disagg, set))

    io = read_from(joinpath("src","readfiles","build","partitioninp.yml"))

    io = partition!(io, set; save = save, overwrite = overwrite)
    cal = calibrate(copy(io), set; save = save, overwrite = overwrite)

    shr = Dict(:va0 => cal[:va0])
    shr = share!(shr, set; save = save, overwrite = overwrite)

    disagg = merge(copy(shr),copy(cal),Dict(
        :r => fill_with((r = set[:r],), 1.0),
        (:yr,:r,:g) => fill_with((yr = set[:yr], r = set[:r], g = set[:g]), 1.0)))

    disagg = disagg!(disagg, set; save = save, overwrite = overwrite)
    return (disagg, set)
end

"""
    build_path(build_step::String; kwargs...)
"""
function build_path(build_step::String; save = true)
    path = save == true ? convert_type(String, today()) : save
    path = joinpath(SLIDE_DIR, "data", path, "build", build_step)
    return path
end

"""
    read_from(path::String)
"""
function read_from(path::String)
    println("  Reading from $path")
    y = read_file(joinpath(SLIDE_DIR, path));
    d = Dict((length(ensurearray(k)) == 1 ? Symbol(k) : Tuple(Symbol.(k))) =>
        read_file(joinpath(y["Path"]..., ensurearray(v)...)) for (k,v) in y["Input"])
    return d
end

"""
    read_set(path::String)
"""
function read_set(path::String)
    cp(joinpath(SLIDE_DIR,"temp","gdpcat.csv"), joinpath(SLIDE_DIR,"data","coresets","gdpcat.csv"), force = true)
    cp(joinpath(SLIDE_DIR,"temp","oth_use.csv"), joinpath(SLIDE_DIR,"data","coresets","i","oth_use.csv"), force = true)
    cp(joinpath(SLIDE_DIR,"temp","cross_fd.csv"), joinpath(SLIDE_DIR,"data","coremaps","crosswalk","fd.csv"), force = true)
    cp(joinpath(SLIDE_DIR,"temp","parse_fd.csv"), joinpath(SLIDE_DIR,"data","coremaps","parse","fd.csv"), force = true)

    d = read_from(path)
    d = Dict(k => sort(v[:,1]) for (k,v) in d)
    return d
end

"""
    write_build(build_step::String, d::Dict; kwargs...)
"""
function write_build(build_step::String, d::Dict; save = false)
    if save !== false
        save_path = build_path(build_step; save = save)

        !isdir(save_path) && mkpath(save_path)
        println("Saving $build_step in $save_path")

        for (k,df) in d
            println("  Saving $k")
            CSV.write(joinpath(save_path, "$k.csv"), df)
        end
    end
end

"""
    read_build(build_step::String; kwargs...)
"""
function read_build(build_step::String; save = true, overwrite::Bool = false)

    save_path = build_path(build_step; save = save)

    if overwrite == true
        d = Dict()
    else
        if !isdir(save_path)
            @warn("Cannot read data from $save_path, as it does not exist.
                \nCalculating $build_step data.")
            d = Dict()
        else
            println("Reading $build_step data from $save_path")
            files = readdir(save_path)
            d = Dict(Symbol(file[1:end-4]) => read_file(joinpath(save_path,file)) for file in files)
        end
    end
    return d
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


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

mutable struct OrchestrationServiceSummary <: SwaggerModel
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName
    serviceState::Any # spec type: Union{ Nothing, String } # spec name: serviceState

    function OrchestrationServiceSummary(;serviceName=nothing, serviceState=nothing)
        o = new()
        validate_property(OrchestrationServiceSummary, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        validate_property(OrchestrationServiceSummary, Symbol("serviceState"), serviceState)
        setfield!(o, Symbol("serviceState"), serviceState)
        o
    end
end # type OrchestrationServiceSummary

const _property_map_OrchestrationServiceSummary = Dict{Symbol,Symbol}(Symbol("serviceName")=>Symbol("serviceName"), Symbol("serviceState")=>Symbol("serviceState"))
const _property_types_OrchestrationServiceSummary = Dict{Symbol,String}(Symbol("serviceName")=>"String", Symbol("serviceState")=>"String")
Base.propertynames(::Type{ OrchestrationServiceSummary }) = collect(keys(_property_map_OrchestrationServiceSummary))
Swagger.property_type(::Type{ OrchestrationServiceSummary }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OrchestrationServiceSummary[name]))}
Swagger.field_name(::Type{ OrchestrationServiceSummary }, property_name::Symbol) =  _property_map_OrchestrationServiceSummary[property_name]

const _allowed_OrchestrationServiceSummary_serviceName = ["AutomaticRepairs"]

const _allowed_OrchestrationServiceSummary_serviceState = ["NotRunning", "Running", "Suspended"]

function check_required(o::OrchestrationServiceSummary)
    true
end

function validate_property(::Type{ OrchestrationServiceSummary }, name::Symbol, val)
    if name === Symbol("serviceName")
        Swagger.validate_param(name, "OrchestrationServiceSummary", :enum, val, _allowed_OrchestrationServiceSummary_serviceName)
    end
    if name === Symbol("serviceState")
        Swagger.validate_param(name, "OrchestrationServiceSummary", :enum, val, _allowed_OrchestrationServiceSummary_serviceState)
    end
end


text:
##################################################
##
## infix logical operators

## XXX Experimental! Not sure these are such a good idea ...
## but used with piecewise
Base.:&(x::Sym, y::Sym) = x.__and__(y)
Base.:|(x::Sym, y::Sym) = x.__or__(y)
Base.:!(x::Sym)         = x.__invert__()

## use ∨, ∧, ¬ for |,&,! (\vee<tab>, \wedge<tab>, \neg<tab>)
∨(x::Sym, y::Sym) = x | y
∧(x::Sym, y::Sym) = x & y
¬(x::Sym) = !x
export ∨, ∧, ¬

## In SymPy, symbolic equations are not represented by `=` or `==`
## rather ther function `Eq` is used. Here we use the unicode
## `\Equal<tab>` for an infix operator. There are also unicode values to represent analogs of `<`, `<=`, `>=`. `>`. These are

## * `<`  is `\ll<tab>`
## * `<=` is `\leqq<tab>
## * `==` is `\Equal<tab>`
## * `>=` is `\geqq<tab>`
## * `>`  is `\gg<tab>`

export ≪,≦,⩵,≧,≫, ≶, ≷ # from CommonEq

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
using PyCall
@pyimport matplotlib.pyplot as plt
using Gen

#######################################
# mixture of two normal distributions #
#######################################

import Distributions
import Gen: random, logpdf, get_static_argument_types
struct TwoNormals <: Distribution{Float64} end
const two_normals = TwoNormals()

function logpdf(::TwoNormals, x, w1, mu1, mu2, sigma1, sigma2)
    if sigma1 < 0 || sigma2 < 0
        return -Inf
    end
    l1 = Distributions.logpdf(Distributions.Normal(mu1, sigma1), x) + log(w1)
    l2 = Distributions.logpdf(Distributions.Normal(mu2, sigma2), x) + log(1 - w1)
    m = max(l1, l2)
    m + log(exp(l1 - m) + exp(l2 - m))
end

function logpdf_grad(::TwoNormals, x, w1, mu1, m2, sigma1, sigma2)
    l1 = Distributions.logpdf(Distributions.Normal(mu1, sigma1), x) + log(w1)
    l2 = Distributions.logpdf(Distributions.Normal(mu2, sigma2), x) + log(1 - w1)
    (deriv_x_1, deriv_mu_1, deriv_sigma_1) = logpdf_grad(normal, x, mu1, sigma1)
    (deriv_x_2, deriv_mu_2, deriv_sigma_2) = logpdf_grad(normal, x, mu2, sigma2)
    w1 = 1.0 / (1.0 + exp(l1 - l2))
    w2 = 1.0 / (1.0 + exp(l2 - l1))
    @assert isapprox(w1 + w2, 1.0)
    deriv_x = deriv_x_1 * w1 + deriv_x_2 * w2
    (deriv_x, NaN, w1 * deriv_mu_1, w2 * deriv_mu_2, w1 * deriv_std_1, w2 * deriv_std_2)
end

function random(::TwoNormals, w1, mu1, mu2, sigma1, sigma2)
    if rand() < w1
        mu1 + sigma1 * randn()
    else
        mu2 + sigma2 * randn()
    end
end

has_output_grad(::TwoNormals) = true
has_argument_grads(::TwoNormals) = (false, true, true, true, true)
get_static_argument_types(::TwoNormals) = [Float64, Float64, Float64, Float64, Float64]

#########
# model #
#########

@gen function model(n::Int)

    # model selection

    if @addr(bernoulli(0.5), :branch)

        # one cluster
        mu = @addr(normal(0, 100), :mu)
        std = @addr(gamma(1, 1), :std)
        for i=1:n
            @addr(normal(mu, std), "y-$i")
        end

    else

        # two clusters
        w1 = @addr(beta(1, 1), :w1)
        mu1 = @addr(normal(0, 100), :mu1)
        mu2 = @addr(normal(0, 100), :mu2)
        std1 = @addr(gamma(1, 1), :std1)
        std2 = @addr(gamma(1, 1), :std2)
        for i=1:n
            @addr(two_normals(w1, mu1, mu2, std1, std2), "y-$i")
        end

    end
end

##################################
# generic trans-dimensional move #
##################################

function resimulation_mh(selection, trace)
    model_args = get_call_record(trace).args
    (new_trace, weight) = regenerate(model, model_args, NoChange(), trace, selection)
    if log(rand()) < weight
        # accept
        return (new_trace, true)
    else
        # reject
        return (trace, false)
    end
end

branch_selection = DynamicAddressSet()
push_leaf_node!(branch_selection, :branch)

function generic_transdim_move(trace)
    resimulation_mh(branch_selection, trace)
end

one_cluster_params = DynamicAddressSet()
push_leaf_node!(one_cluster_params, :mu)
push_leaf_node!(one_cluster_params, :std)

two_cluster_params = DynamicAddressSet()
push_leaf_node!(two_cluster_params, :mu1)
push_leaf_node!(two_cluster_params, :mu2)
push_leaf_node!(two_cluster_params, :std1)
push_leaf_node!(two_cluster_params, :std2)

w1_selection = DynamicAddressSet()
push_leaf_node!(w1_selection, :w1)

function fixed_dim_move(trace)
    if get_assignment(trace)[:branch]
        (trace, accept) = resimulation_mh(one_cluster_params, trace)
        mala(model, one_cluster_params, trace, 0.01)
    else
        (trace, accept) = resimulation_mh(two_cluster_params, trace)
        mala(model, two_cluster_params, trace, 0.01)
    end
    trace
end

###############
# RJMCMC move #
###############

@gen function split_proposal(prev)
    u1 = @addr(beta(2, 2), :u1)
    u2 = @addr(beta(2, 2), :u2)
    u3 = @addr(beta(1, 1), :u3)
end

@gen function merge_proposal(prev) end

const MODEL = :model
const PROPOSAL = :proposal

@inj function split_injection(n::Int)
    @assert @read(:model => :branch)
    @write(false, :model => :branch)
    mu = @read(:model => :mu)
    std = @read(:model => :std)
    u1 = @read(:proposal => :u1)
    u2 = @read(:proposal => :u2)
    u3 = @read(:proposal => :u3)
    w1 = u1
    w2 = 1.0 - u1

    # always make mu1 < mu2
    mu1 = mu - u2 * std * sqrt(w2 / w1)
    mu2 = mu + u2 * std * sqrt(w1 / w2)

    var = std * std
    std1 = sqrt(u3 * (1 - u2 * u2) * var / w1)
    std2 = sqrt((1 - u3) * (1 - u2 * u2) * var / w2)

    @write(mu1, :model => :mu1)
    @write(std1, :model => :std1)
    @write(mu2, :model => :mu2)
    @write(std2, :model => :std2)
    @write(w1, :model => :w1)
    for i=1:n
        @copy(:model => "y-$i", :model => "y-$i")
    end
end

@inj function merge_injection(n::Int)
    @assert !@read(:model => :branch)
    @write(true, :model => :branch)
    w1 = @read(:model => :w1)
    w2 = 1 - w1
    mu1 = @read(:model => :mu1)
    mu2 = @read(:model => :mu2)
    std1 = @read(:model => :std1)
    std2 = @read(:model => :std2)
    u1 = w1

    # make mu1 <= mu2
    if mu1 > mu2
        tmp = mu1
        mu1 = mu2
        mu2 = tmp
    end

    @write(u1, :proposal => :u1)
    mu = w1 * mu1 + w2 * mu2
    @write(mu, :model => :mu)
    var = w1 * (mu1 * mu1 + std1 * std1) + w2 * (mu2 * mu2 + std2 * std2) - mu * mu
    std = sqrt(var)
    @write(std, :model => :std)
    u2 = (mu - mu1) / (std * sqrt(w2 / w1))
    std_ratio = std1 / std2
    weighted_var_ratio = std_ratio * std_ratio * w1 / w2
    u3 = weighted_var_ratio / (1 + weighted_var_ratio)
    @write(u2, :proposal => :u2)
    @write(u3, :proposal => :u3)
    for i=1:n
        @copy(:model => "y-$i", :model => "y-$i")
    end
end

correction = (new_trace) -> 0.
function rjmcmc_transdim_move(trace)
    branch = get_assignment(trace)[:branch]
    n = get_call_record(trace).args[1]
    if branch
        rjmcmc(model,
            split_proposal, (),
            merge_proposal, (),
            split_injection, (n,),
            trace, correction)
    else
        rjmcmc(model,
            merge_proposal, (),
            split_proposal, (),
            merge_injection, (n,),
            trace, correction)
    end
end

########################
# inference experiment #
########################

Gen.load_generated_functions()

import Random
Random.seed!(1)

# load data set

import CSV
function load_data_set()
    df = CSV.read("coal.csv")
    ys = df[1]
    ys = ys - mean(ys)
    ys = ys ./ std(ys)
    return ys
end

#const ys = load_data_set()
ys = Float64[-1.0, 1.0]
const n = length(ys)

function plot_trace_plot()
    plt.figure(figsize=(6, 5))

    observations = DynamicAssignment()
    for (i, y) in enumerate(ys)
        observations["y-$i"] = y
    end

    # RJMCMC
    (trace, _) = generate(model, (n,), observations)
    model_choice_vec = Bool[]
    burn_in = 10000
    scores = Float64[]
    for iter=1:burn_in + 2000
        (trace, accept) = rjmcmc_transdim_move(trace)
        trace = fixed_dim_move(trace)
        score = get_call_record(trace).score
        if iter > burn_in
            push!(model_choice_vec, get_assignment(trace)[:branch])
            push!(scores, score)
        end
        println("iter: $iter, score: $score")
    end
    plt.subplot(2, 1, 1)
    plt.plot(convert(Vector{Int}, model_choice_vec) .+ 1, "r")
    plt.ylabel("Selected Model")
    plt.title("Custom Reversible Jump Trans-Dimensional Moves")

    # generic
    (trace, _) = generate(model, (n,), observations)
    model_choice_vec = Bool[]
    burn_in = 10000
    scores = Float64[]
    for iter=1:burn_in + 2000
        (trace, accept) = generic_transdim_move(trace)
        trace = fixed_dim_move(trace)
        score = get_call_record(trace).score
        if iter > burn_in
            push!(model_choice_vec, get_assignment(trace)[:branch])
            push!(scores, score)
        end
        println("iter: $iter, score: $score")
    end
    plt.subplot(2, 1, 2)
    plt.plot(convert(Vector{Int}, model_choice_vec) .+ 1, "r")
    plt.ylabel("Selected Model")
    plt.title("Generic Trans-Dimensional Moves")

    ax = plt.gca()
    plt.tight_layout()
    plt.savefig("mixture_trace_plot.png")
end

plot_trace_plot()

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
# dependency package

include("../packages/DependencyTree/src/DependencyTree.jl")

@test all(cur_char ->
  isdefined(
    DependencyTree,
    Symbol("Type$(cur_char)")
  ),
  'A':'Z'
)


text:
function sqlite3_close(dbptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_close, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
    @unix_only return ccall(
        (:sqlite3_close, sqlite3_lib),
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
end

function sqlite3_close_v2(dbptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_close_v2, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
    @unix_only return ccall(
        (:sqlite3_close_v2, sqlite3_lib),
        Cint,
        (Ptr{Void}, ),
        dbptr
    )
end