text:
function ergodic_distribution(model, sol)
    return ergodic_distribution(model, sol.dr, sol.dr.grid_exo, sol.dr.grid_endo, sol.dprocess)
end

function ergodic_distribution(model, dr, exo_grid:: UnstructuredGrid, endo_grid:: CartesianGrid, dp)
    N_m = n_nodes(exo_grid)
    N_s = n_nodes(endo_grid)
    N = N_m*N_s
    Π = zeros(N_m, N_s, N_m, endo_grid.n...)
    s = nodes(endo_grid)
    p = SVector(model.calibration[:parameters]...)
    a = SVector(endo_grid.min...)
    b = SVector(endo_grid.max...)
    for i_m in 1:n_nodes(exo_grid)
        x = dr(i_m, s)
        m = node(exo_grid, i_m)
        for i_M in 1:n_inodes(dp, i_m)
            M = inode(Point, dp, i_m, i_M)
            w = iweight(dp, i_m, i_M)
            S = transition(model, m, s, x, M, p)
            S = [(S[n]-a)./(b-a) for n=1:length(S)]
            if typeof(Π) <: AbstractArray{Float64,4}
                trembling_hand!(view(Π,i_m, :, i_M, :), S, w)
            elseif typeof(Π) <: AbstractArray{Float64,5}
                trembling_hand!(view(Π,i_m, :, i_M, :, :), S, w)
            end
        end
    end
    Π0 = (reshape(Π,N,N))' - Matrix(I, N, N)
    Π0[end,:] .= 1
    B = zeros(N)
    B[end] = 1.0
    μ = Π0\B

    μ = reshape(μ, N_m, N_s)
    return Π, μ
end

function ergodic_distribution(model, dr, exo_grid:: CartesianGrid, endo_grid:: CartesianGrid, dp)
    N_m = n_nodes(exo_grid)
    N_s = n_nodes(endo_grid)
    N = N_m*N_s
    Π = zeros(N_m, N_s, exo_grid.n..., endo_grid.n...)
    s = nodes(endo_grid)
    p = SVector(model.calibration[:parameters]...)
    a = SVector(exo_grid.min..., endo_grid.min...)
    b = SVector(exo_grid.max..., endo_grid.max...)
    for i_m in 1:n_nodes(exo_grid)
        x = dr(i_m, s)
        m = node(exo_grid, i_m)
        for i_M in 1:n_inodes(dp, i_m)
            M = inode(Point, dp, i_m, i_M)
            w = iweight(dp, i_m, i_M)
            S = transition(model, m, s, x, M, p)
            V = [(SVector(M..., el...)-a)./(b.-a) for el in S]
            # return Π, i_m, V, w
            if typeof(Π) <: AbstractArray{Float64,4}
                trembling_hand!(view(Π,i_m, :, :, :), V, w)
            elseif typeof(Π) <: AbstractArray{Float64,3}
                trembling_hand!(view(Π,i_m, :, :), V, w)
            end
        end
    end

    Π0 = (reshape(Π,N,N))' - Matrix(I, N, N)
    Π0[end,:] .= 1
    B = zeros(N)
    B[end] = 1.0
    μ = Π0\B

    μ = reshape(μ, N_m, N_s)
    return reshape(Π, N, N), μ
end

function ergodic_distribution(model, dr, exo_grid:: EmptyGrid, endo_grid:: CartesianGrid, dp)
    N_m = 1
    N_s = n_nodes(endo_grid)
    N = N_m*N_s
    Π = zeros(N_s, endo_grid.n...)
    s = nodes(endo_grid)
    p = SVector(model.calibration[:parameters]...)
    a = SVector(endo_grid.min...)
    b = SVector(endo_grid.max...)
    i_m = 1
    x = dr(s)
    m = SVector(model.calibration[:exogenous]...)
    for i_M in 1:n_inodes(dp, i_m)
        M = inode(Point, dp, i_m, i_M)
        w = iweight(dp, i_m, i_M)
        S = transition(model, m, s, x, M, p)
        S = [(S[n]-a)./(b-a) for n=1:length(S)]
        trembling_hand!(Π, S, w)
    end
    Π0 = (reshape(Π,N,N))' - Matrix(I, N, N)
    Π0[end,:] .= 1
    B = zeros(N)
    B[end] = 1.0
    μ = Π0\B

    μ = reshape(μ, endo_grid.n...)
    return reshape(Π, N, N), μ
end

function trembling_hand!(A::AbstractArray{Float64,2}, x, w)
    N,n0 = size(A)
    δ0 = 1.0./(n0-1.0)
    for n in 1:N
        x0 = x[n][1]
        x0 = min.(max.(x0, 0.0),1.0)
        q0 = div.(x0, δ0)
        q0 = max.(0, q0)
        q0 = min.(q0, n0-2)
        λ0 = (x0./δ0-q0) # ∈[0,1[ by construction
        q0_ = round.(Int,q0) + 1
        A[n, q0_]   += (1-λ0)*w
        A[n, q0_+1] += λ0*w
    end
end

function trembling_hand!(A::AbstractArray{Float64,3}, x, w)
    N,n0,n1 = size(A)
    δ0 = 1.0./(n0-1.0)
    δ1 = 1.0./(n1-1.0)
    for n in 1:N
        x0 = x[n][1]
        x0 = min.(max.(x0, 0.0),1.0)
        q0 = div.(x0, δ0)
        q0 = max.(0, q0)
        q0 = min.(q0, n0-2)
        λ0 = (x0./δ0-q0) # ∈[0,1[ by construction
        q0_ = round.(Int,q0) + 1

        x1 = x[n][2]
        x1 = min.(max.(x1, 0.0),1.0)
        q1 = div.(x1, δ1)
        q1 = max.(0, q1)
        q1 = min.(q1, n1-2)
        λ1 = (x1./δ1-q1) # ∈[0,1[ by construction
        q1_ = round.(Int,q1) + 1

        A[n, q0_ ,  q1_] += (1-λ0)*(1-λ1)*w
        A[n, q0_+1, q1_] += λ0*(1-λ1)*w
        A[n, q0_, q1_+1] += (1-λ0)*λ1*w
        A[n, q0_+1, q1_+1] += λ0*λ1*w
    end
end


text:
import StatsFuns: log2π

@node MvNormalWeightedMeanPrecision Stochastic [ out, (ξ, aliases = [ xi, weightedmean ]), (Λ, aliases = [ invcov, precision ]) ]

conjugate_type(::Type{ <: MvNormalWeightedMeanPrecision }, ::Type{ Val{ :Λ } })   = Wishart

@average_energy MvNormalWeightedMeanPrecision (q_out::Any, q_ξ::Any, q_Λ::Any) = begin
    (m_mean, v_mean) = mean_cov(q_ξ)
    (m_out, v_out)   = mean_cov(q_out)
    0.5 * (ndims(q_out) * log2π + logdet(cholinv(mean(q_Λ))) + tr(mean(q_Λ)*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)')))
end

@average_energy MvNormalWeightedMeanPrecision (q_out_ξ::Any, q_Λ::Any) = begin
    (m, V) = mean_cov(q_out_ξ)
    d = Int64(ndims(q_out_ξ)/2)
    @views 0.5*(d*log2π + logdet(cholinv(mean(q_Λ))) + tr(mean(q_Λ)*( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end,1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ) ))
end

ext: jl


text:
# Julia 1.0.0
# 字符、字符串的使用

ch = 'k'
println(ch)
# 转换为数值
println(Int(ch))

str = "THIS a test str"
# 换行、不换行输出
println(str)
print(str)
print(str, "\n")
# 大小写转换
println("uppercase: ", uppercase(str))
println("lowercase: ", lowercase(str))
# 子串
show(str[1:4]); println()
println(str[6:7])
println(str[end-3:end])

# julia 支持字符串插值用法
a = "sb"
b = "jiakong"
println("$a $b")
# 插值用法可以扩展到数学计算
println("2 + 22 = $(2 + 22)")

# 字符串拼接
c = string(a, " ", b)
println(c)
# 可以用 *（非 +）实现字符串拼接
d = a * " + " * b
print(d)

max_stars_repo_name: NopSmile/Blog

max_issues_repo_name: NopSmile/Blog

max_forks_repo_name: NopSmile/blog


text:
using Parameters
include("common.jl")

# Background current parameters
@with_kw struct BGParams
    G_CAB = 3.22E-3  # Conductance of Background Calcium current
    G_NAB = 5.45E-4  # Conductance of Background Sodium current
end

# Background Ca current (iCaB) (μA/cm²)
icab(vm, eCa, G_CAB) = G_CAB * (vm - eCa)

# Background Na current (iNaB) (μA/cm²)
inab(vm, eNa, G_NAB) = G_NAB * (vm - eNa)

# Non-specific Ca-activated current parameters
@with_kw struct NSParams
    # Na permebility of non-specific current (cm/s)
    P_NA = 1.75E-7
    # K permebility of non-specific current (cm/s)
    P_K = 0.
    # Ca Michaelis constant for non-specific current (mM)
    KM_CA = 1.2E-3
end

# Non-specific Ca-activated Na current (iNsNa) (uA/cm^2)"
function insna(na_i, ca_i, vfrt, evfrtm1, na_o, KM_CA, P_NA)
    imax = 0.75 * P_NA * F * vfrt * (na_i * (evfrtm1 + 1) - na_o) / evfrtm1
    insna = imax * _hills(ca_i, KM_CA, 3)
end

function insna(na_i, ca_i, vfrt, evfrtm1, na_o, pNS::NSParams)
    @unpack KM_CA, P_NA = pNS
    insna(na_i, ca_i, vfrt, evfrtm1, na_o, KM_CA, P_NA)
end


text:
################################################################################
#
# NfOrd/Ideal/Arithmetic.jl : Arithmetic for ideals in orders of absolute
#                             number fields
#
# This file is part of Hecke.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann
#  Copyright (C) 2015, 2016, 2017 Claus Fieker
#
################################################################################

function check_parent(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
   if order(x) !== order(y)
     error("Ideals do not have the same order.")
   end
end

################################################################################
#
#  Ideal addition / GCD
#
################################################################################

#x has a princ gen special
function sum_princ_gen_special(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  OK = order(x)
  genx = x.princ_gen_special[2]+x.princ_gen_special[3]
  if has_2_elem(y)
    if has_minimum(y)
      gen1 = gcd(genx, minimum(y, copy = false))
      if gen1 == minimum(y, copy = false)
        return y
      end
    else
      gen1 = gcd(genx, y.gen_one)
      if gen1 == y.gen_one
        return y
      end
    end
    res = ideal(OK, gen1, y.gen_two)
    fl = false
    if has_minimum(y)
      fl = gen1 == ppio(minimum(y, copy = false), genx)[1]
    else
      fl = gen1 == ppio(y.gen_one, genx)[1]
    end
    if fl
      if has_norm(y)
        res.norm = ppio(norm(y), genx)[1]
      end
      if has_minimum(y)
        res.minimum = ppio(minimum(y), genx)[1]
      end
    end
  else
    M1 = _hnf_modular_eldiv(basis_matrix(y, copy = false), genx, :lowerleft)
    res = ideal(OK, M1, false, true)
  end
  @hassert :NfOrd 1 res == sum_via_basis_matrix(x, y)
  return res
end

function sum_via_basis_matrix(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  OK = order(x)
  d = degree(OK)
  g = gcd(minimum(x, copy = false), minimum(y, copy = false))
  H = vcat(basis_matrix(x, copy = false), basis_matrix(y, copy = false))
  hnf_modular_eldiv!(H, g, :lowerleft)
  H = view(H, (d + 1):2*d, 1:d)
  res = ideal(OK, H, false, true)
  if isone(basis(OK, copy = false)[1])
    res.minimum = H[1, 1]
  end
  res.norm = prod_diagonal(H)
  return res
end

#x is a principal ideal
function sum_princ_gen(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  OK = order(x)
  d = degree(OK)
  g = gcd(minimum(y, copy = false), minimum(x, copy = false))
  M = basis_matrix(y, copy = false)
  N = representation_matrix_mod(x.princ_gen, g)
  L = vcat(M, N)
  hnf_modular_eldiv!(L, g, :lowerleft)
  H = view(L, (d + 1):2*d, 1:d)
  res = ideal(OK, H, false, true)
  if isone(basis(OK, copy = false)[1])
    res.minimum = H[1, 1]
  end
  res.norm = prod_diagonal(H)
  @hassert :NfOrd 1 res == sum_via_basis_matrix(x, y)
  return res
end

@doc Markdown.doc"""
    +(x::NfOrdIdl, y::NfOrdIdl)

Returns $x + y$.
"""
function +(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  check_parent(x, y)
  OK = order(x)

  if iszero(x)
    return y
  end
  if iszero(y)
    return x
  end

  if isdefined(x, :gen_one) && isdefined(y, :gen_one) && isone(gcd(x.gen_one, y.gen_one))
    return ideal(OK, 1)
  end
  if isdefined(x, :norm) && isdefined(y, :norm) && isone(gcd(x.norm, y.norm))
    return ideal(OK, 1)
  end
  if has_princ_gen_special(x)
    return sum_princ_gen_special(x, y)
  end
  if has_princ_gen_special(y)
    return sum_princ_gen_special(y, x)
  end
  g = gcd(minimum(x, copy = false), minimum(y, copy = false))
  if isone(g)
    return ideal(order(x), g)
  end
  d = degree(OK)
  if issimple(nf(OK)) && isdefining_polynomial_nice(nf(OK)) && contains_equation_order(OK) && isprime(g) && !isindex_divisor(OK, g) && has_2_elem(x) && has_2_elem(y)
    #I can use polynomial arithmetic
    if fits(Int, g)
      R1 = ResidueRing(FlintZZ, Int(g), cached = false)
      R1x = PolynomialRing(R1, "x", cached = false)[1]
      ggp_small = gcd(R1x(x.gen_two.elem_in_nf), R1x(y.gen_two.elem_in_nf))
      if isone(ggp_small)
        return ideal(OK, 1)
      end
      Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
      ggZ = lift(Zx, ggp_small)
    else
      R = ResidueRing(FlintZZ, g, cached = false)
      Rx = PolynomialRing(R, "x", cached = false)[1]
      ggp_large = gcd(Rx(x.gen_two.elem_in_nf), Rx(y.gen_two.elem_in_nf))
      if isone(ggp_large)
        return ideal(OK, 1)
      end
      Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
      ggZ = lift(Zx, ggp_large)
    end
    gen_2 = OK(nf(OK)(ggZ))
    return ideal(OK, g, gen_2)
  end
  if has_princ_gen(x) && !has_basis_matrix(x)
    return sum_princ_gen(x, y)
  end
  if has_princ_gen(y) && !has_basis_matrix(y)
    return sum_princ_gen(y, x)
  end
  return sum_via_basis_matrix(x, y)
end

################################################################################
#
#  Intersection / LCM
#
################################################################################

@doc Markdown.doc"""
    intersect(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl
    lcm(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl

Returns $x \cap y$.
"""
function intersect(x::NfOrdIdl, y::NfOrdIdl)
  check_parent(x, y)
  g = gcd(minimum(x), minimum(y))
  if isone(g)
    #The ideals are coprime, the intersection is equal to the product
    return x*y
  end
  d = degree(order(x))
  H = vcat(basis_matrix(x, copy = false), basis_matrix(y, copy = false))
  K = left_kernel(H)[2]
  g = lcm(minimum(x),minimum(y))
  return ideal(order(x), _hnf_modular_eldiv(view(K, 1:d, 1:d)*basis_matrix(x, copy = false), g, :lowerleft), false, true)
end

@doc Markdown.doc"""
    intersect(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl
    lcm(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl

Returns $x \cap y$.
"""
lcm(x::NfOrdIdl, y::NfOrdIdl) = intersect(x, y)

################################################################################
#
#  Ideal multiplication
#
################################################################################

@doc Markdown.doc"""
    *(x::NfOrdIdl, y::NfOrdIdl)

Returns $x \cdot y$.
"""
function *(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  check_parent(x, y)
  @hassert :NfOrd 1 isconsistent(x)
  @hassert :NfOrd 1 isconsistent(y)
  OK = order(x)
  if ismaximal_known_and_maximal(OK)
    z = mul_maximal(x, y)
  else
    z = mul_gen(x, y)
  end
  @hassert :NfOrd 1 isconsistent(z)
  return z
end

function mul_gen(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  check_parent(x, y)
  if iszero(x) || iszero(y)
    I1 = ideal(order(x), 0)
    I1.iszero = 1
    return I1
  end
  O = order(x)
  d = degree(O)
  l = minimum(x, copy = false)*minimum(y, copy = false)
  if has_2_elem(x) && has_basis_matrix(y)
    M1 = representation_matrix_mod(x.gen_two, l)
    Mf = vcat(minimum(x, copy = false)*basis_matrix(y, copy = false), basis_matrix(y, copy = false)*M1)
    hnf_modular_eldiv!(Mf, l, :lowerleft)
    J = ideal(O, view(Mf, (d+1):2*d, 1:d), false, true)
    if iscoprime(minimum(x, copy = false), minimum(y, copy = false))
      J.minimum = minimum(x, copy = false)*minimum(y, copy = false)
    end
    return J
  end
  if has_2_elem(y) && has_basis_matrix(x)
    M1 = representation_matrix_mod(y.gen_two, l)
    Mf = vcat(minimum(y, copy = false)*basis_matrix(x, copy = false), basis_matrix(x, copy = false)*M1)
    hnf_modular_eldiv!(Mf, l, :lowerleft)
    J = ideal(O, view(Mf, (d+1):2*d, 1:d), false, true)
    if iscoprime(minimum(x, copy = false), minimum(y, copy = false))
      J.minimum = minimum(x, copy = false)*minimum(y, copy = false)
    end
    return J
  end
  z = zero_matrix(FlintZZ, 2*degree(O), degree(O))
  z1 = zero_matrix(FlintZZ, 2*degree(O), degree(O))
  X = basis(x, copy = false)
  Y = basis_matrix(y, copy = false)
  for i in 1:d
    M1 = representation_matrix_mod(X[i], l)
    _copy_matrix_into_matrix(z1, 1, 1, M1)
    hnf_modular_eldiv!(z1, minimum(x, copy = false), :lowerleft)
    mul!(M1, Y, M1)
    _copy_matrix_into_matrix(z, 1, 1, M1)
    hnf_modular_eldiv!(z, l, :lowerleft)
    if view(z1, d+1:2*d, 1:d) == basis_matrix(x, copy = false)
      break
    end
  end
  # This is a d^2 x d matrix
  J = ideal(O, view(hnf_modular_eldiv!(z, l, :lowerleft),
                      (d+1):2*d, 1:d), false, true)
  if iscoprime(minimum(x, copy = false), minimum(y, copy = false))
    J.minimum = minimum(x, copy = false)*minimum(y, copy = false)
  end
  return J
end

# using the 2-normal representation
function prod_via_2_elem_normal(a::NfAbsOrdIdl, b::NfAbsOrdIdl)
  check_parent(a, b)
  @hassert :NfOrd 1 has_2_elem_normal(a)
  @hassert :NfOrd 1 has_2_elem_normal(b)
  O = order(a)
  a1 = a.gen_one
  if has_minimum(a)
    a1 = minimum(a, copy = false)
  end
  b1 = b.gen_one
  if has_minimum(b)
    b1 = minimum(b, copy = false)
  end
  m = lcm(a1, b1)
  e, f = ppio(m, a1)
  if f == 1
    a2 = a.gen_two
  else
    # we need to become normal at m, we are normal at a
    # higher powers in a are fine
    # CRT: the 2nd gen of a needs to stay the same at a
    # and should be  1 at f

    #a2 = a.gen_two*f + a1^2
    mul!(e, a1, a1)
    a2 = f*a.gen_two
    add!(a2, a2, e)
    # now (a1, a2) should be m-normal for a
  end

  e, f = ppio(m, b1)
  if f == 1
    b2 = b.gen_two
  else
    #b2 = b.gen_two*f + b1^2
    mul!(e, b1, b1)
    b2 = f*b.gen_two
    add!(b2, b2, e)
  end
  C = ideal(O, a1*b1, a2*b2)
  C.norm = norm(a) * norm(b)

  if has_minimum(a) && has_minimum(b)
    ma = minimum(a, copy = false)
    mb = minimum(b, copy = false)

    if gcd(ma, mb) == 1
      C.minimum = ma * mb
    end

    # Otherwise, I don't know the
    # correct value.
  end

  C.gens_normal = m
  C.gens_weakly_normal = true #for the 1-ideal: it will not be normal...

  return C
end

# using the 2-weak-normal representation
function prod_via_2_elem_weakly(a::NfAbsOrdIdl, b::NfAbsOrdIdl)
  check_parent(a, b)
  @hassert :NfOrd 1 has_2_elem(a)
  @hassert :NfOrd 1 has_2_elem(b)

  O = order(a)
  K = nf(O)
  bas = basis(O)
  n = degree(O)

  norm_c = norm(a, copy = false) * norm(b, copy = false)        # we ARE in the maximal order case
  first_gen_new = fmpz(1)
  if has_minimum(a)
    first_gen_new *= minimum(a, copy = false)
  else
    first_gen_new *= gcd(a.gen_one, norm(a, copy = false))
  end
  if has_minimum(b)
    first_gen_new *= minimum(b, copy = false)
  else
    first_gen_new *= gcd(b.gen_one, norm(b, copy = false))
  end
  norm_int_c = norm_c
  mod_c = fmpz(1)

  if has_minimum(a)
    mod_c *= minimum(a)
  else
    mod_c *= norm(a)
  end

  if has_minimum(b)
    mod_c *= minimum(b)
  else
    mod_c *= norm(b)
  end

  # Is this a good idea? Bad magic constants

  if mod_c > 250
    r = -500:500  # should still have enough elements...
  else
    r = -Int(div(mod_c, 2)):Int(div(mod_c, 2))
  end

  @vprint :NfOrd 1 "a: $a \nb: $b"
  @vprint :NfOrd 1 "using basis: $bas"

  gen = O()
  gen2 = O()
  t = O()
  s = O()
  u = O()

  cnt = 0
  while true
    cnt += 1
    if cnt % 20 == 0
      assure_2_normal(a)
      assure_2_normal(b)
      return a*b
    end
    rand!(t, O, r) # Nemo.rand_into!(bas, r, t)
    r2 = rand(r)
    rand!(s, O, r) # Nemo.rand_into!(bas, r, s)
    r4 = rand(r)
    mul!(t, t, a.gen_two)      # Nemo.mult_into!(t, a.gen_two, t)
    add!(t, t, r2*a.gen_one)   # Nemo.add_into!(t, r2*a.gen_one, t)
    mul!(s, s, b.gen_two)      # Nemo.mult_into!(s, b.gen_two, s)
    add!(s, s, r4*b.gen_one)   # Nemo.add_into!(s, r4*b.gen_one, s)
    mul!(u, s, t)              # Nemo.mult_into!(s, t, u)
    add!(gen, u, gen)          # Nemo.add_into!(u, gen, gen)
#    gen2 += (r1*K(a.gen_two) + r2*a.gen_one) *
#           (r3*K(b.gen_two) + r4*b.gen_one)
    gen = mod(gen, mod_c^2)    # = element_reduce_mod(gen, O, FlintZZ(mod_c)^2)

    if gcd(norm(gen), norm_int_c^2) == norm_int_c # should be ^n, but for
                                                  # the test ^2 is sufficient
      break
    end
  end

  @vprint :NfOrd 1 "prod_via_2_elem: used $cnt tries\n"

  c = ideal(O, first_gen_new, gen)

  c.norm = norm_c

  if has_minimum(a) && has_minimum(b) && gcd(minimum(a), minimum(b)) == 1
    c.minimum = minimum(a) * minimum(b)
                    # otherwise, I don't know the correct value
  end

  c.gens_weakly_normal = true

  return c
end

# dispatching
@doc Markdown.doc"""
    *(x::NfOrdIdl, y::NfOrdIdl)

Returns the ideal x*y.
"""
function mul_maximal(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  check_parent(x, y)
  if iszero(x) || iszero(y)
    z = ideal(order(x), 0)
    z.iszero = 1
    return z
  end
  if isone(x)
    return y
  elseif isone(y)
    return x
  end
  if has_princ_gen_special(x)
    return princ_gen_special(x)*y
  end
  if has_princ_gen_special(y)
    return princ_gen_special(y)*x
  end
  if has_2_elem_normal(x) && has_2_elem_normal(y)
    return prod_via_2_elem_normal(x, y)
  end
  if has_2_elem(x) && has_2_elem(y)
    return prod_via_2_elem_weakly(x, y)
  end
  # Fall back to the generic algorithm _mul(::NfOrdIdl, ::NfOrdIdl)
  # Could also use invoke
  return mul_gen(x, y)
end

################################################################################
#
#  Addition
#
################################################################################

# Falls back to generic case +(::NfOrd, ::NfOrd)
#for ideals in the maximal order, the gcd is well defined...

@doc Markdown.doc"""
    gcd(A::NfOrdIdl, B::NfOrdIdl) -> NfOrdIdl

The gcd or sum (A+B).
"""
function gcd(A::NfAbsOrdIdl, B::NfAbsOrdIdl)
  check_parent(A, B)
  return A+B
end

function gcd_into!(A::NfAbsOrdIdl, B::NfAbsOrdIdl, C::NfAbsOrdIdl)
  return C+B
end

#TODO: write a ppio version that allows for p-powers as well
@doc Markdown.doc"""
    gcd(A::NfOrdIdl, p::fmpz) -> NfOrdIdl

The gcd or sum (A + pO).
"""
function gcd(A::NfAbsOrdIdl, p::fmpz)
  if isdefined(A, :minimum)
    if gcd(A.minimum, p) == 1
      return ideal(order(A), fmpz(1))
    end
  end
  if isdefined(A, :norm)
    if gcd(A.norm, p) == 1
      return ideal(order(A), fmpz(1))
    end
  end
  if has_2_elem(A)
    g = gcd(A.gen_one, p)
    return ideal(order(A), g, A.gen_two)
  else
    @assert isdefined(A, :basis_matrix)
    return A + ideal(order(A), p)
  end
end

################################################################################
#
#  Powering
#
################################################################################

function Base.:(^)(A::NfAbsOrdIdl, e::Int)
  @hassert :NfOrd 1 isconsistent(A)
  OK = order(A)
  if e == 0
    return ideal(OK, 1)
  elseif e == 1
    return A
  end
  if ismaximal_known_and_maximal(OK) && issimple(nf(OK)) && has_2_elem(A)
    return pow_2_elem(A, e)
  else
    return Base.power_by_squaring(A, e)
  end
end

function pow_2_elem(A::NfAbsOrdIdl, e::Int)
  OK = order(A)
  if A.is_principal == 1 && isdefined(A, :princ_gen)
    gen = (A.princ_gen)^e
    I = ideal(OK, gen)
    if isdefined(A, :norm)
      I.norm = norm(A)^e
    end
    if isprime_known(A) && isprime(A)
      eA = A.splitting_type[1]
      I.minimum = minimum(A)^(div(e-1, eA)+1)
    end
    I.is_prime = 2
    return I
  elseif isprime_known(A) && isprime(A)
    eA = A.splitting_type[1]
    minim = minimum(A)^(div(e-1, eA)+1)
    gen1 = minim
    gen2 = A.gen_two^e
    I = ideal(OK, gen1, gen2)
    I.minimum = deepcopy(minim)
    if isdefined(A, :norm)
      I.norm = norm(A)^e
    end
    if has_2_elem_normal(A)
      I.gens_normal = A.gens_normal
    end
    I.gens_weakly_normal = has_weakly_normal(A)
    I.is_prime = 2
    return I
  else
    gen1 = A.gen_one^e
    gen2 = A.gen_two^e
    I = ideal(OK, gen1, gen2)
    if isdefined(A, :norm)
      I.norm = norm(A)^e
    end
    if isdefined(A, :minimum) && isone(gcd(minimum(A, copy = false), discriminant(OK)))
      I.minimum = A.minimum^e
    end
    if has_2_elem_normal(A)
      I.gens_normal = A.gens_normal
    end
    I.gens_weakly_normal = has_weakly_normal(A)
    I.is_prime = 2
    return I
  end
end

# To stop the wrong julia behavior for I^2 and I^-2
Base.literal_pow(::typeof(^), A::NfAbsOrdIdl, ::Val{p}) where {p} = A^p

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

# multiplication by fmpz, using two normal presentation
function prod_by_int_2_elem_normal(A::NfOrdIdl, a::fmpz)
  @assert has_2_elem(A) && has_2_elem_normal(A)

  # <a,a> is a a-normal presentation
  # we need to have a common presentation for A and a though

  m = lcm(a, A.gen_one)

  e, f = ppio(m, A.gen_one)
  if f == 1
    a2 = A.gen_two
  else
                           # we need to become normal at m, we are normal at a
                           # higher powers in a are fine
                           # CRT: the 2nd gen of a needs to stay the same at a
                           # and should be  1 at f
    a2 = A.gen_two*f + A.gen_one^2 # now (a1, a2) should be m-normal for a
  end

  B = NfOrdIdl(A.gen_one*a, a2*a)
  if has_princ_gen_special(A)
    B.princ_gen_special = (2, 0, (A.princ_gen_special[2] + A.princ_gen_special[3])*a)
  end

  if has_princ_gen(A)
    B.is_principal = 1
    B.princ_gen = a*A.princ_gen
  end
  B.gens_normal = m

  if has_minimum(A)
    B.minimum = A.minimum * a
  end

  if has_norm(A)
    B.norm = A.norm * a^degree(order(A))
  end

  @assert has_2_elem(B) && has_2_elem_normal(B)
  return B
end

function prod_by_int_2_elem(A::NfOrdIdl, a::fmpz)
  @assert has_2_elem(A)

  B = NfOrdIdl(A.gen_one*a, A.gen_two*a)
  if has_princ_gen(A)
    B.is_principal = 1
    B.princ_gen = A.princ_gen*a
  end

  if has_princ_gen_special(A)
    B.princ_gen_special = (2, 0, (A.princ_gen_special[2]+A.princ_gen_special[3])*a)
  end

  if has_minimum(A)
    B.minimum = A.minimum * a
  end

  if has_norm(A)
    B.norm = A.norm * a^degree(order(A))
  end

  @assert has_2_elem(B)
  return B
end

@doc Markdown.doc"""
    *(x::NfOrdIdl, y::fmpz) -> NfOrdIdl

Returns the ideal $x \cdot y$.
"""
function *(x::NfOrdIdl, y::fmpz)
  if ismaximal_known_and_maximal(order(x))
    return mul_maximal(x, y)
  else
    return mul_gen(x, y)
  end
end

function mul_maximal(x::NfOrdIdl, y::fmpz)
  if iszero(y)
    z = ideal(order(x), 0)
    z.iszero = 1
    return z
  end

  if isone(y) || y == -1
    return x
  end

  if has_2_elem(x)
    if has_2_elem_normal(x)
      return prod_by_int_2_elem_normal(x,y)
    else
      return prod_by_int_2_elem(x,y)
    end
  end

  @hassert :NfOrd 1 has_basis_matrix(x)
  return ideal(order(x), basis_matrix(x, copy = false)*y)
end

*(x::fmpz, y::NfOrdIdl) = y * x

@doc Markdown.doc"""
    *(x::NfOrdIdl, y::Integer) -> NfOrdIdl

Returns the ideal $x \cdot y$.
"""
*(x::NfOrdIdl, y::Integer) = x * fmpz(y)

*(x::Integer, y::NfOrdIdl) = y * x

function *(x::NfOrdElem, y::NfOrdIdl)
  parent(x) !== order(y) && error("Orders of element and ideal must be equal")
  return ideal(parent(x), x) * y
end

*(x::NfOrdIdl, y::NfOrdElem) = y * x

function mul_gen(x::NfOrdIdl, y::fmpz)
  if y == 0
    z = ideal(order(x), zero_matrix(FlintZZ, degree(order(x)), degree(order(x))))
    z.iszero = 1
    return z
  end

  z = ideal(order(x), basis_matrix(x, copy = false)*y)
  if isdefined(x, :princ_gen)
    z.princ_gen = x.princ_gen * y
  end
  if isdefined(x, :princ_gen_special)
    z.princ_gen_special = (2, 0, x.princ_gen_special[x.princ_gen_special[1] + 1] * y)
  end
  if isdefined(x, :norm)
    z.norm = x.norm*y^degree(order(x))
  end
  if isdefined(x, :minimum)
    z.minimum = x.minimum*y
  end
  return z
end

################################################################################
#
#  Idempotents
#
################################################################################

@doc Markdown.doc"""
    idempotents(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdElem, NfOrdElem

Returns a tuple `(e, f)` consisting of elements `e in x`, `f in y` such that
`1 = e + f`.

If the ideals are not coprime, an error is raised.
"""
function idempotents(x::NfAbsOrdIdl, y::NfAbsOrdIdl)
  check_parent(x, y)
  O = order(x)
  d = degree(O)

 if has_2_elem(x) && has_2_elem(y)
    g, ux, vx = gcdx(x.gen_one, y.gen_one)
    if isone(g)
      z = O(ux*x.gen_one)
      @hassert :NfOrd 2 z in x
      @hassert :NfOrd 2 (1 - z) in y
      return z, 1 - z
    end
  end

  mx = minimum(x, copy = false)
  my = minimum(y, copy = false)

  g, ux, vy = gcdx(mx, my)
  if isone(g)
    z = O(ux*mx)
    @hassert :NfOrd 2 z in x
    @hassert :NfOrd 2 (1 - z) in y
    return z, 1 - z
  end
  return _idempotents_via_matrices(x, y)
end

function _idempotents_via_matrices(x::NfAbsOrdIdl, y::NfAbsOrdIdl)

  O = order(x)
  d = degree(O)
  mx = minimum(x, copy = false)
  my = minimum(y, copy = false)
  # form the matrix
  #
  # ( 1 |  1  | 0 )
  # ( 0 | M_x | I )
  # ( 0 | M_y | 0 )

  V = O.tidempotents

  u = coordinates(one(O))

  V[1, 1] = 1

  for i in 1:d
    V[1, i + 1] = u[i]
  end

  _copy_matrix_into_matrix(V, 2, 2, basis_matrix(x, copy = false))
  _copy_matrix_into_matrix(V, 2 + d, 2, basis_matrix(y, copy = false))

  for i in 1:d
    V[1 + i, d + 1 + i] = 1
  end

  H = hnf_modular_eldiv!(V, lcm(mx, my)) # upper right

  for i in 2:(1 + d)
    if H[1, i] != 0
      error("Ideals are not coprime")
    end
  end

  basisofx = basis(x, copy = false)

  aux = O()
  z = basisofx[1]*H[1, d + 2]
  for i in 2:d
    mul!(aux, basisofx[i], H[1, d + 1 + i])
    add!(z, z, aux)
  end

  @hassert :NfOrd 2 -z in x
  @hassert :NfOrd 2 1 + z in y

  ccall((:fmpz_mat_zero, libflint), Nothing, (Ref{fmpz_mat}, ), V)

  return -z, 1 + z
end

################################################################################
#
#  crt
#
################################################################################

@doc Markdown.doc"""
    crt(r1::NfOrdElem, i1::NfOrdIdl, r2::NfOrdElem, i2::NfOrdIdl) -> NfOrdElem

Find $x$ s.th. $x \equiv r1 \bmod i1$ and $x \equiv r2 \bmod i2$
using `idempotents`.
"""
function crt(r1::S, i1::T, r2::S, i2::T) where { S <: Union{NfOrdElem, NfRelOrdElem, AlgAssAbsOrdElem}, T <: Union{NfOrdIdl, NfRelOrdIdl, AlgAssAbsOrdIdl} }
  u, v = idempotents(i1, i2)
  return r1*v + r2*u
end

function crt(a::Vector{S}, I::Vector{T}) where { S <: Union{NfOrdElem, NfRelOrdElem, AlgAssAbsOrdElem}, T <: Union{NfOrdIdl, NfRelOrdIdl, AlgAssAbsOrdIdl} }
  if length(a) == 1
    return a[1]
  end
  if length(a) == 2
    return crt(a[1], I[1], a[2], I[2])
  end
  A = [crt(a[2*i-1], I[2*i-1], a[2*i], I[2*i]) for i=1:div(length(a), 2)]
  B = [I[2*i-1]*I[2*i] for i=1:div(length(a), 2)]
  if isodd(length(a))
    push!(A, a[end])
    push!(B, I[end])
  end
  return crt(A, B)
end

################################################################################
#
#  Division
#
################################################################################

divexact(A::NfAbsOrdIdl, b::Integer) = divexact(A, fmpz(b))

#TODO: write a divexact! to change the ideal?
#  difficult due to Julia's inability to unset entries...

@doc Markdown.doc"""
    divexact(A::NfOrdIdl, y::fmpz) -> NfOrdIdl
    divexact(A::NfOrdIdl, y::Integer) -> NfOrdIdl

Returns $A/y$ assuming that $A/y$ is again an integral ideal.
"""
function divexact(A::NfAbsOrdIdl, b::fmpz)
  if iszero(A)
    return A
  end
  zk = order(A)
  b = abs(b)
  if has_2_elem(A)
    B = ideal(zk, divexact(A.gen_one, b), divexact(A.gen_two, b))
    if isdefined(A, :gens_normal)
      B.gens_normal = A.gens_normal
    end
    B.gens_weakly_normal = A.gens_weakly_normal
    if has_basis_matrix(A)
      B.basis_matrix = divexact(A.basis_matrix, b)
    end
    if false && has_basis_mat_inv(A)
      error("not defined at all")
      B.basis_mat_inv = b*A.basis_mat_inv
    end
  else
    B = ideal(zk, divexact(A.basis_matrix, b))
    if false && has_basis_mat_inv(A)
      error("not defined at all")
      B.basis_mat_inv = b*A.basis_mat_inv
    end
  end
  if has_minimum(A)
    B.minimum = divexact(A.minimum, b)
  end
  if has_norm(A)
    B.norm = divexact(A.norm, b^degree(zk))
  end
  if has_princ_gen(A)
    B.princ_gen = divexact(A.princ_gen, b)
  end
  #TODO princ_gen_special missing
  return B
end

@doc Markdown.doc"""
    divexact(A::NfOrdIdl, B::NfOrdIdl) -> NfOrdIdl

Returns $AB^{-1}$ assuming that $AB^{-1}$ is again an integral ideal.
"""
function divexact(A::NfOrdIdl, B::NfOrdIdl)
  check_parent(A, B)
  # It is assumed that B divides A, that is, A \subseteq B
  t_prod = 0.0
  t_simpl = 0.0
  t_b_mat = 0.0
  t_2_elem_weak = 0.0
  t_2_elem = 0.0

  if norm(A, copy = false) == norm(B, copy = false)
    return ideal(order(A), one(FlintZZ), order(A)(1))
  else
    t_prod += @elapsed I = A*inv(B)
    t_simpl += @elapsed simplify_exact!(I)
    #t_b_mat += @elapsed B = basis_matrix(I)
    I.den != 1 && error("Division not exact")
    #I = ideal(order(A), B.num)
    #t_2_elem_weak += @elapsed _assure_weakly_normal_presentation(I)
    #t_2_elem += @elapsed assure_2_normal(I)

    #println("  computation for product: $t_prod")
    #println("  simplification         : $t_simpl")
    #println("  basis matrix           : $t_b_mat")
    #println("  2 weak presentation    : $t_2_elem_weak")
    #println("  2 elem presentation    : $t_2_elem")

    return I.num

  end
end

################################################################################
#
#  Extend/contract
#
################################################################################

function extend(A::NfOrdIdl, O::NfOrd)
  if order(A) === O
    return A
  end
  # Assumes order(A) \subseteq O

  if iszero(A)
    B = ideal(O, zero_matrix(FlintZZ, degree(O), degree(O)))
    B.iszero = 1
    return B
  end

  d = degree(O)
  M = zero_matrix(FlintZZ, d^2, d)
  X = basis(O, copy = false)
  Y = map(O, basis(A, copy = false))
  t = O()
  for i = 1:d
    for j = 1:d
      mul!(t, X[i], Y[j])
      for k = 1:d
        M[(i - 1)*d + j, k] = coordinates(t, copy = false)[k]
      end
    end
  end
  M = sub(_hnf_modular_eldiv(M, minimum(A), :lowerleft), d*(d - 1) + 1:d^2, 1:d)
  return ideal(O, M, false, true)
end

*(A::NfOrdIdl, O::NfOrd) = extend(A, O)
*(O::NfOrd, A::NfOrdIdl) = extend(A, O)

function contract(A::NfOrdIdl, O::NfOrd)
  if order(A) === O
    return A
  end
  # Assumes O \subseteq order(A)

  if iszero(A)
    B = ideal(O, zero_matrix(FlintZZ, degree(O), degree(O)))
    B.iszero = 1
    return B
  end

  d = degree(O)
  M = basis_matrix(O, copy = false)*basis_mat_inv(order(A), copy = false)
  @assert M.den == 1
  H = vcat(basis_matrix(A, copy = false), M.num)
  K = left_kernel(H)[2]
  M = view(K, 1:d, 1:d)*basis_matrix(A, copy = false)
  M = M*basis_matrix(order(A), copy = false)*basis_mat_inv(O, copy = false)
  @assert M.den == 1
  M = _hnf_modular_eldiv(M.num, minimum(A), :lowerleft)
  res = ideal(O, M, false, true)
  if A.is_prime == 1
    res.is_prime = 1
  end
  return res
end

intersect(O::NfOrd, A::NfOrdIdl) = contract(A, O)

intersect(A::NfOrdIdl, O::NfOrd) = contract(A, O)


text:
# helpers for wrapping the library

function unsafe_message(ptr, args...)
    str = unsafe_string(ptr, args...)
    API.LLVMDisposeMessage(ptr)
    str
end

## defining types in the LLVM type hierarchy

# llvm.org/docs/doxygen/html/group__LLVMCSupportTypes.html

# macro that adds an inner constructor to a type definition,
# calling `refcheck` on the ref field argument
macro checked(typedef)
    # decode structure definition
    if Meta.isexpr(typedef, :macrocall)
        # handle `@compat` prefixing 0.6-style type declarations
        typedef = macroexpand(typedef)
    end
    if Meta.isexpr(typedef, :struct)
        structure = typedef.args[2]
        body = typedef.args[3]
    else
        error("argument is not a structure definition")
    end
    if isa(structure, Symbol)
        # basic type definition
        typename = structure
    elseif Meta.isexpr(structure, :<:)
        # typename <: parentname
        all(e->isa(e,Symbol), structure.args) ||
            error("typedef should consist of plain types, ie. not parametric ones")
        typename = structure.args[1]
    else
        error("malformed type definition: cannot decode type name")
    end

    # decode fields
    field_names = Symbol[]
    field_defs = Union{Symbol,Expr}[]
    for arg in body.args
        if isa(arg, LineNumberNode)
            continue
        elseif isa(arg, Symbol)
            push!(field_names, arg)
            push!(field_defs, arg)
        elseif Meta.isexpr(arg, :(::))
            push!(field_names, arg.args[1])
            push!(field_defs, arg)
        end
    end
    :ref in field_names || error("structure definition should contain 'ref' field")

    # insert checked constructor
    push!(body.args, :(
        $typename($(field_defs...)) = (refcheck($typename, ref); new($(field_names...)))
    ))

    return esc(typedef)
end

# the most basic check is asserting that we don't use a null pointer
@inline function refcheck(::Type, ref::Ptr)
    ref==C_NULL && throw(UndefRefError())
end


text:
using StructOpt
using DataStructures, LinearAlgebra, PGFPlotsX

function main()
    ## build Lasso pball reg problem, with solution on partly smooth manifold ||x||=1, dimension 2.
    problems = []

    ## Lasso 2d
    sf = 1.0
    A = sf * Matrix{Float64}([1. 1.; 0. 1.])
    y = sf * Vector([1., 2.])
    λ = 1.0
    pb = LassoPb{regularizer_l1}(A, y, λ, size(A, 2), Vector([0, 1]))
    objopt, xopt = get_optimsol(pb)
    @show objopt, xopt
    objopt = 1.5

    xstart = Vector([4, 2]) + 50*Vector([1, 1])

    push!(problems, (
        # name = "Lasso-2d-IstaFista",
        name = "Lasso-2d",
        pb = pb,
        x0 = xstart,
        xopt = xopt,
        objopt = objopt,
    ))

    ## Build algorithms
    algorithms = []

    αuser = 0.1
    push!(algorithms, (
        name="ISTA",
        updatefunc=extra_ISTA,
        params=Dict(
            :αuser => αuser,
            :printstep => 1,
        ),
        # params=Dict(),
    ))

    push!(algorithms, (
        name="FISTA",
        updatefunc=extra_FISTA,
        # params=Dict(:αuser=>αuser),
        params=Dict(
            :αuser => αuser,
            :printstep => 1,
        ),
    ))

    push!(algorithms, (
        name="T1",
        updatefunc=extra_CondInertia,
        params=Dict(
            :id_testkind=>:default,
            :αuser => αuser,
            :printstep => 1,
        ),
    ))

    push!(algorithms, (
        name="T2",
        updatefunc=extra_CondPredInertia,
        params=Dict(
            :printstep => 1,
            :αuser => αuser,
            :printstep => 1,
            :proxgrad_per_it => 2,
        ),
    ))

    push!(algorithms, (
        name="MFISTA",
        updatefunc=extra_MFISTA,
        params=Dict(
            :αuser => αuser,
            :printstep => 1,
        ),
    ))

    ## Start from interesting point and optimize with ISTA, FISTA.
    problem_to_algstats = OrderedDict()
    for (pb_id, pb_data) in enumerate(problems)
        printstyled("\n----- Solving problem $pb_id, $(pb_data.name)\n", color=:light_blue)

        algo_to_stats = OrderedDict()
        for (algo_id, algo) in enumerate(algorithms)
            printstyled("\n---- Algorithm $algo_id, $(algo.name)\n", color=:light_blue)

            xopt, hist = solve_proxgrad(pb_data.pb, pb_data.x0, algo.updatefunc; algo.params...)
            algo_to_stats[algo] = (hist=hist, xopt=xopt)
        end

        problem_to_algstats[pb_data] = algo_to_stats
    end

    ####################################################
    ## Display suboptimality, iterates position.
    FIGS_FOLDER = "./figs"
    basename(pwd()) == "src" && (FIGS_FOLDER = joinpath("..", FIGS_FOLDER))
    !ispath(FIGS_FOLDER) && mkpath(FIGS_FOLDER)

    for (pb, algo_to_stats) in problem_to_algstats
        println("Generating graphs for ", pb.name)

        ## Suboptimality plot
        build_save_suboptimality(pb, algo_to_stats, FIGS_FOLDER)

        ## Steplength plot
        # build_save_steplength(pb, algo_to_stats, FIGS_FOLDER)

        if pb.pb.n == 2
            reg = typeof(pb.pb).parameters[1]

            ## Iterates position
            ps = []

            xmin, xmax = -1, 4
            ymin, ymax = -1, 2

            ## Plot contour
            x = xmin:(xmax - xmin) / 100:xmax
            y = ymin:(ymax - ymin) / 100:ymax
            f = (x,y) -> 0.5 * norm(pb.pb.A * [x, y] - pb.pb.y)^2
            add_contour!(ps, f, x, y)

            for (algo, stats) in algo_to_stats
                coords = [ (point[1], point[2]) for (iter, point) in stats.hist[:iter_x]]

                add_traj!(ps, coords, algo.name, params=get_iterates_algoparams(algo.name))
            end

            ## Plot optimal point
            add_point!(ps, xopt)

            ## Plot manifolds
            coords = [(xmin, 0), (xmax, 0)]; add_manifold!(ps, coords)
            coords = [(0, ymin), (0, ymax)]; add_manifold!(ps, coords)

            fig = pgf_build_iteratesfig(ps, xmin, xmax, ymin, ymax)

            pgfsave(joinpath(FIGS_FOLDER, "$(pb.name)_iterates.pdf"), fig)
            pgfsave(joinpath(FIGS_FOLDER, "$(pb.name)_iterates.tikz"), fig)
        end
    end

    println("\n")
    for (pb, algo_to_stats) in problem_to_algstats
        @show pb.name
        for (algo, stats) in algo_to_stats
            @show algo.name
            display(stats[:xopt])
        end
    end

    return
end

main()


text:
using Revise
using Pkg

Base.exit_on_sigint(false)

try
    Revise.async_steal_repl_backend()
    includet("Jus.jl")
    dir=joinpath(dirname(Pkg.project().path), "src")
    map(f-> Revise.track(Jus, joinpath(dir, f)),
        ["types.jl", "vars.jl", "Jus.jl", "server.jl"])
    Jus.exec(ARGS) do config, ws
        Revise.revise()
        Base.invokelatest(Jus.serve, config, ws)
    end
catch err
    !(err isa InterruptException) && @warn "Error" exception=(err, catch_backtrace())
    exit(1)
end


text:
# This file currently contains both the objective function and preprocessing
# steps, plus helper functions SSE and norm

using Debug

# some global parameters -
BIG = 1e10
num_objs = 11 # changes for cross validation
SMALL = 1e-3

###################
try
@show start_random
catch
global start_random = false
@show start_random
end
###################
# How we calculate error
function SSE(model_values,data_values)
  y = model_values
  m = data_values
  # (max(y)-10)^2
  E = sum((m-y).^2)
end

# This doesn't work properly because the SSEs are being calculated on normalized data
# function SSE(model_values,data_values,Mscale=250)
#   y = model_values
#   m = data_values
#   # (max(y)-10)^2
#   E = sum((m-y).^2) + ((Mscale-maximum(y))/Mscale).^2
# end

# How we normalize the data
function norm_dat(data)
  min = minimum(data)
  max = maximum(data)

#  if max == 0
  #  max = 1
#  end
if sum(data) > SMALL
  normalized_data = (data-min)/(max-min)
else
  normalized_data = zeros(data)
end
end

#     #    #   # Prep objective data
# Load and Normalize each blot individually
fpath="../poets/medici_blot_combined.jl"
include(fpath)
#Fig1B = norm_dat(Fig1B)
#Fig1F = norm_dat(Fig1F)
#Fig2A = norm_dat(Fig2A)
#Fig2E = norm_dat(Fig2E)
#Fig3D = norm_dat(Fig3D)
#Fig3F = norm_dat(Fig3F)
#Fig4D = norm_dat(Fig4D)
#Fig4E = norm_dat(Fig4E)
#Fig4G = norm_dat(Fig4G)
#Fig5A = norm_dat(Fig5A)
#Fig5C = norm_dat(Fig5C)

#     #    #   # Initialize and  Run model to steady state
 # EXECUTE RunAllSimulations UP TO LOOP

 # Load data dictionary off disk

 # Initialize simulation functions
 fpath="../src/new_SimulationFunctions.jl" # Loads up DataFile , SolveBalances , Balances, Kinetics, Control
 include(fpath)

 # Establish simulation parameters
 TSTART = 0.0;
 Ts = 1.0; # 0.1
 TSTOP = 250.0; # 100 ; 500

 # Initialize intial data dictionary
 data_dictionary = DataFile(TSTART,TSTOP,Ts)

 # # # #

 # update control parameter array here for now
 if start_random == false
 GAIN_MASTER = 0.3
 num_regs = 19;
 control_parameter_array=data_dictionary["CONTROL_PARAMETER_ARRAY"]
 for i in 1:num_regs
   control_parameter_array[i,1]=GAIN_MASTER # Gain for transcriptional activity
 end
 # TGFB3_inhibition: TGFB3_Ab inhibition; TGFB3_binding
 control_parameter_array[1,1]=30.0 # GAIN # _tgfb3_Ab
 # NFATc_inhibition: vivit inhibition; Calcineurin_activation_NFATc
 control_parameter_array[2,1]=3.0 # GAIN # NFATc_inhibition: vivit inhibition; Calcineurin_activation_NFATc
 # GSK3_P_inhibition_Binding_APC_AXIN_forw: GSK3_P inhibition; Binding_APC_AXIN_forw
 control_parameter_array[3,1]=0.5 # GAIN # GSK3_P_inhibition_Binding_APC_AXIN_forw
 control_parameter_array[3,2]=1.0 # ORDER# GSK3_P_inhibition_Binding_APC_AXIN_forw

 # Active_TCF4_inhibition_Ecad: Ecadherin_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
 control_parameter_array[4,1]=0.1 # GAIN # Active_TCF4_inhibition_Ecad
 control_parameter_array[4,1]=0.9 # GAIN # Active_TCF4_inhibition_Ecad
 # Active_TCF4_inhibition_APC: APC_AXIN_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
 control_parameter_array[5,1]=0.2 # GAIN # Active_TCF4_inhibition_APC
 control_parameter_array[5,2]=1.0 # ORDER # Active_TCF4_inhibition_APC

 # activation_Ecadherin_by_Active_NFATc: Active_NFATc induction induction_gene_Ecadherin
 control_parameter_array[6,1]=0.3 # GAIN # activation_Ecadherin_by_Active_NFATc
 # repression_Ecadherin_by_pSmad: pSmad repression induction_gene_Ecadherin
 control_parameter_array[7,1]=5.0 # GAIN # REP_ECAD_PSMAD
 control_parameter_array[7,2]=2.0 # ORDER # REP_ECAD_PSMAD
 # repression_Ecadherin_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_Ecadherin
 control_parameter_array[8,1]=0.5 # GAIN_SNAIL_SLUG_ECAD #2.0
 # repression_Ecadherin_by_Active_LEF1: Active_LEF1 repression induction_gene_Ecadherin
 control_parameter_array[9,1]=0.6# GAIN_ACTIVE_LEF1_ECAD #2.0 0.5
 control_parameter_array[9,2]=1.0# ORDER_ACTIVE_LEF1_ECAD
 #  control_parameter_array[6,1]=0.3 # GAIN_BCAT_TGFB3
 # activation_Ecadherin_by_ecad_virus: ecad_virus induction induction_gene_Ecadherin
 control_parameter_array[10,1]=5.0 # ecad_virus ***
 # activation_TGFB3_by_TCF4: Active_TCF4 induction induction_gene_TGFB3
 control_parameter_array[11,1]=0.15 # GAIN_TCF4_TGFB3
 control_parameter_array[11,2]=1.5 # GAIN_TCF4_TGFB3

 # repression_TGFB3_by_YREG1: YREG1 repression induction_gene_TGFB3
 control_parameter_array[12,1]=1.0 # YREG_TGFB3 #1.0
 control_parameter_array[12,2]=1.0 # YREG_TGFB3
 # repression_YREG1_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_YREG1
 control_parameter_array[13,1]=0.1 # GAIN_SNAIL_YREG1 # 0.1
 control_parameter_array[13,2]=2.0 # ORDER_SNAIL_YREG1 # 2.0
 # activation_LEF1_by_pSmad: pSmad induction induction_gene_LEF1
 control_parameter_array[14,1]=0.8# GAIN # _PSMAD_LEF1 # 0.5 # 0.8
 control_parameter_array[14,2]=2.0# ORDER# _PSMAD_LEF1       # 1.0, 2.0
 # repression_LEF1_by_YREG1: YREG1 repression induction_gene_LEF1
 control_parameter_array[15,1]=5.0 # GAIN YREG1_LEF1 # 0.1
 control_parameter_array[15,2]=0.5 # ORDER YREG1_LEF1 # 2.0
 # activation_vimentin_by_Active_LEF1: Active_LEF1 induction induction_gene_vimentin
 control_parameter_array[16,1]=0.3# GAIN_LEF1_VIM
 control_parameter_array[16,2]=2.0# ORDER_LEF1_VIM
 # activation_SNAIL_SLUG_by_AP1_SP1_P: AP1_SP1_P induction induction_gene_SNAIL_SLUG
 control_parameter_array[18,1]=1.0# SS by AP1 SP1
 # activation_SNAIL_SLUG_by_ss_virus: ss_virus induction induction_gene_SNAIL_SLUG
 control_parameter_array[19,1]=5.0# SS_VIRUS ***

# # # #

 data_dictionary["CONTROL_PARAMETER_ARRAY"]=control_parameter_array
 # Also use edit_file_global here

 fpath = "../src/edit_file_global.jl"
 include(fpath)
 UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
 UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
end

 println("Running to steady state ... ")

@debug function stop()
  @bp
end
#@show data_dictionary["RATE_CONSTANT_ARRAY"][15]

 # # # STEADY STATE
 println(now())
 data_dict_cpa_global = deepcopy(data_dictionary)
 XSS = EstimateSteadyState(data_dict_cpa_global) #!! scaling
 println(now())
 initial_condition_array = vec(XSS);
 data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
 #xss_data_dict = deepcopy(data_dictionary) # may not be necessary
#@show data_dictionary["RATE_CONSTANT_ARRAY"][15]

#     #    #   # Prep data dictionary for rapid deconstruction / reassembly
len_cpa = length(data_dictionary["CONTROL_PARAMETER_ARRAY"])
len_rca = length(data_dictionary["RATE_CONSTANT_ARRAY"])
len_sca = length(data_dictionary["SATURATION_CONSTANT_ARRAY"])

# We cannot simply dump each array into the parameter array, since each
# sub array as a different shape (especially the SCA)
 # get lengths of data dictoinary components to be included in simulation

 # SCA is sparse, so let's gather linear indices for its non-zero terms
 i = 1 #index to start at
 sca_indices,sca_values = [],[]
 for i in 1:len_sca
   val = data_dictionary["SATURATION_CONSTANT_ARRAY"][i]
 if val !=0
   push!(sca_indices,i)
   push!(sca_values,val)
 end
 end
 # package for quick reconstruction
 # this needs to be reassembled at runtime anyway :/
 sca_pack = zip(sca_indices,sca_values)

  # CPA, let's give access to gain AND power
  cpa_gain = data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] # Gain
  cpa_order= data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] # Power

# Initialize parameter_array
 parameter_array = []
 append!(parameter_array,data_dictionary["RATE_CONSTANT_ARRAY"])
 append!(parameter_array,cpa_gain)
 append!(parameter_array,cpa_order)
 append!(parameter_array,sca_values)
original_parameter_array = deepcopy(parameter_array)
 # get lengths of new arrays
 len_cpa_gain = length(cpa_gain)
 len_cpa_order = length(cpa_order)
 len_sca_values = length(sca_values)
 # array order RCA,CPAGAIN,CPAORDER,SCA
   # need indices for reconstruction
len_sca_prime = length(sca_values)
const  A = len_rca + len_cpa_gain
const  B = A + len_cpa_order
const  C = B + len_sca_prime

# record ribosome and rnap terms
  # pa97=parameter_array[97]
  # pa98=parameter_array[98]
  # pa99=parameter_array[99]
  # pa100=parameter_array[100]
  # pa101=parameter_array[101]
  # pa102=parameter_array[102]

# record and lock addition terms to their initial values

  const addition_terms_saved = parameter_array[85:100]

 #  const pa117=parameter_array[117]
 #  const pa118=parameter_array[118]
#println(parameter_array)

# parameter_array gets sent to estimate estimate_ensemble as initial_state

# Evaluate blot values (objective targets) before loadin objective function

include("../src/spin_up_edits.jl")

#Evaluates the objective function values -
function objective_function(parameter_array)
  #@show data_dictionary["RATE_CONSTANT_ARRAY"][15]

 # Alias the species -
 # in order of appearance below!

# Descriptions include the CURRENT indexed value,
# NOT necessarily the correct one
# Indices are from the DataFile and are NOT automously updated!!!
 # LEF 1 protein (technically this not not total, due to active amt)
 index_lef = 58
 # LEF 1 active
 index_lef_act = 59
 # Snail_slug total protein
 index_ss = 62
 # ecadherin protein (unbound)
 index_ecad = 51
 # ecad_bcat index
 index_ecad_bcat = 53
 # vimentin protein
 index_vim = 63
 # tgfb3 protein
 index_tgfb3 = 92
 # This is fit to ACTIVE_TCF4, data is IP from bcat and tcf-4
 index_bcat_tcf = 57 # i.e. this index is for Active_TCF4
 # Should probably be total protein, so sum index_bcat with index_tcf
 index_tcf = 56
# GSK3_P
 index_gsk3_p = 44
 #
 #

 # Alias species as they will appear in the output array

# indices for simulation data array sim_data
# 1 "edit_file_none.jl"
# 2 "edit_file_tgfb.jl"
# 3 "edit_file_mapkin.jl"
# 4 "edit_file_ss.jl"
# 5 "edit_file_ssecad.jl"
# 6 "edit_file_ss_tgfb3ab.jl"
# 7 "edit_file_ss_dnlef1.jl"
# 8 "edit_file_tgfb3.jl"
# 9 "edit_file_dnsmad.jl"
# 10 "edit_file_tgfb3_dnlef1.jl"
# 11 "edit_file_tgfb_tgfb3ab.jl"# !!! THIS HAS BEEN MODIFIED TO BE TGFB3 +PI3K INHIBITOR

#  MOD,BLOT = norm_dat(MOD),norm_dat(BLOT)

 # Execute simulations

 # Calculate the objective function array -
 obj_array = BIG*ones(num_objs,1)

 # Refactor data dictionary for model execution
 # need to specify implicit strucutre of parameter_array
 # structure: rca;cpa(1);cpa(2);sca
 data_dictionary["RATE_CONSTANT_ARRAY"] = parameter_array[1:len_rca]
 data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] = parameter_array[len_rca+1:A]
 data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] = parameter_array[A+1:B]
# Rebuild SCA
 sca_pack = zip(sca_indices,parameter_array[B+1:C]) # new values
 for (i,val) in sca_pack
 data_dictionary["SATURATION_CONSTANT_ARRAY"][i]=val
 end
 #@show data_dictionary["RATE_CONSTANT_ARRAY"][15]

### RUN TO NEW STEADY STATE!
# Make new xss_data_dict

# # # STEADY STATE
# WAIT, what initial condiiton array do we use???
# use last steady state, we should still be stored in the data_dictionary object
# as long as it is global
# # # # THIS DATA DICTIONARY CONTAINS THE MASTER STATE
# # # # FOR A GIVEN SET OF PARAMETERS

data_dict_from_param_array = deepcopy(data_dictionary)
XSS = EstimateSteadyState(data_dict_from_param_array) # BE CAREFUL WITH DATA DICTIONARY NAME
initial_condition_array = vec(XSS);
data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
xss_data_dict = deepcopy(data_dictionary) # BE CAREFUL WITH DATA DICTIONARY NAME

### EXECUTE RunAllSimulations loop (without file writes)

# ??? This should be read off of ram SolveModelGeneric needs to be updated
#simulations = ["edit_file_none.jl","edit_file_tgfb.jl","edit_file_mapkin.jl","edit_file_ss.jl","edit_file_ssecad.jl","edit_file_ss_tgfb3ab.jl","edit_file_ss_dnlef1.jl","edit_file_tgfb3.jl","edit_file_dnsmad.jl","edit_file_tgfb3_dnlef1.jl","edit_file_tgfb_tgfb3ab.jl"]
sim_data = [] # Keep data in RAM
for sim in simulations
  ##### What data_dictionary should be loaded here??? clearly the active one
  data_dictionary_active = deepcopy(xss_data_dict) # RENAME THIS TO RETAIN ORIGINAL call it data_dictionary_active
  (TSIM,X) =  SolveModelGenericPOETS_ram(TSTART,TSTOP,Ts,data_dictionary_active,sim)
  push!(sim_data,X)
  # now reset the data dictionary
end

###
t48 = 149 # 48 hours(24hrs = 264)

 #  sim_data[1][526,55] -> [t48,index_species]
 # specific objectives
 # EXP [conditions] RO [readouts]
 #1 EXP SS RO SS LEF # Fig1B

 MOD  = [sim_data[1][t48,index_ss],sim_data[1][t48,index_lef]+sim_data[1][t48,index_lef_act],
         sim_data[4][t48,index_ss],sim_data[4][t48,index_lef]+sim_data[4][t48,index_lef_act]]
 BLOT = [Fig1B[1],Fig1B[7],
         Fig1B[2],Fig1B[9]]
#@show MOD
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[1]  = SSE(MOD,BLOT)

#@show BLOT,MOD #SHOWCODEOBJ
 #@bp
#@show MOD

 #2 EXP SS DNLEF RO ECAD VIM # Fig1F
 MOD  = [sim_data[1][t48,index_ecad]+sim_data[1][t48,index_ecad_bcat],sim_data[1][t48,index_vim],
         sim_data[4][t48,index_ecad]+sim_data[4][t48,index_ecad_bcat],sim_data[4][t48,index_vim],
         sim_data[7][t48,index_ecad]+sim_data[7][t48,index_ecad_bcat],sim_data[7][t48,index_vim]]
 BLOT = [Fig1F[1],Fig1F[5],
         Fig1F[2],Fig1F[6],
         Fig1F[4],Fig1F[8]]

 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
#@show MOD
 obj_array[2]  = SSE(MOD,BLOT)
#@show BLOT,MOD #SHOWCODEOBJ
 #@bp
 #3 EXP SS RO TGFB3 # Fig2A
 MOD  = [sim_data[1][t48,index_tgfb3],
         sim_data[4][t48,index_tgfb3],
         sim_data[5][t48,index_tgfb3]]
 BLOT = [Fig2A[1],
         Fig2A[2],
         Fig2A[4]]
#@show MOD
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[3]  = SSE(MOD,BLOT)
#@show MOD

 # println(MOD,' ',BLOT,' ',SSE(MOD,BLOT))
 # if MOD == [1.0,0.0]
 #   #@bp
 #   println("WARNING Obj3 is BROKEN")
 #   MOD = [0.0,1.0]
 #   @show parameter_array
 # end
 #@show BLOT,MOD #SHOWCODEOBJ
 #4 EXP SS TGFBab RO LEF ECAD VIM # Fig2E #  index_ecad_bcat
 MOD  = [sim_data[1][t48,index_ecad]+sim_data[1][t48,index_ecad_bcat],sim_data[1][t48,index_vim],sim_data[1][t48,index_lef]+sim_data[1][t48,index_lef_act],
         sim_data[4][t48,index_ecad]+sim_data[4][t48,index_ecad_bcat],sim_data[4][t48,index_vim],sim_data[4][t48,index_lef]+sim_data[4][t48,index_lef_act],
         sim_data[6][t48,index_ecad]+sim_data[6][t48,index_ecad_bcat],sim_data[6][t48,index_vim],sim_data[6][t48,index_lef]+sim_data[6][t48,index_lef_act]]
 BLOT = [Fig2E[6],Fig2E[11],Fig2E[1],
         Fig2E[7],Fig2E[12],Fig2E[2],
         Fig2E[10],Fig2E[15],Fig2E[5]]
         #@show MOD
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[4] = SSE(MOD,BLOT)
#@show BLOT-MOD #SHOWCODEOBJ
#@show MOD
 #5 EXP SS ECAD RO BCAT TCF4 # Fig3D
 # CAREFUL in setting index_bcat & index_tcf, total BCAT & TCF?
 # consistent with total TCF4 and free/nuclear bcat, but need to verify
 MOD  = [sim_data[1][t48,index_bcat_tcf],(sim_data[1][t48,index_tcf]+sim_data[1][t48,index_bcat_tcf])/6,
         sim_data[4][t48,index_bcat_tcf],(sim_data[4][t48,index_tcf]+sim_data[4][t48,index_bcat_tcf])/6,
         sim_data[5][t48,index_bcat_tcf],(sim_data[5][t48,index_tcf]+sim_data[5][t48,index_bcat_tcf])/6]
         # dividing by 6 compensates for non-active TCF4 that is in excess in the model (otherwise the normalization is skewed)
         # This implies that the relative amount of TCF4 is too much compared to BETA-CATEININ OR the blot is not measuring
         # a quantitative deifference between TCF4 and TCF4+BCAT
 BLOT = [Fig3D[1],Fig3D[6],
         Fig3D[3],Fig3D[8],
         Fig3D[5],Fig3D[10]]
         #@show MOD
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[5]  = SSE(MOD,BLOT)
#@show BLOT,MOD #SHOWCODEOBJ
#@show MOD
 #6 EXP SS ECAD RO TGFB3 ECAD VIM # Fig3F #  index_ecad_bcat
 MOD  = [sim_data[1][t48,index_tgfb3],sim_data[1][t48,index_vim],sim_data[1][t48,index_ecad]+sim_data[1][t48,index_ecad_bcat],
         sim_data[4][t48,index_tgfb3],sim_data[4][t48,index_vim],sim_data[4][t48,index_ecad]+sim_data[4][t48,index_ecad_bcat],
         sim_data[5][t48,index_tgfb3],sim_data[5][t48,index_vim],sim_data[5][t48,index_ecad]+sim_data[5][t48,index_ecad_bcat]]
 BLOT = [Fig3F[1],Fig3F[11],Fig3F[6],
         Fig3F[2],Fig3F[12],Fig3F[7],
         Fig3F[5],Fig3F[15],Fig3F[10]]
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[6]  = SSE(MOD,BLOT)
#@show MOD
#@show BLOT,MOD #SHOWCODEOBJ
 #7 EXP TGFB3 RO SS #Fig4D
 #MOD  = [sim_data[1][t48,index_ss],
#         sim_data[8][t48,index_ss]]
 #BLOT = [Fig4D[1],
#         Fig4D[3]]
 #MOD = norm_dat(MOD)
 #BLOT = norm_dat(BLOT)

 #7 EXP TGFB3 RO SS #Fig4D
 MOD  = [sim_data[1][t48,index_gsk3_p],
         sim_data[8][t48,index_gsk3_p],
         sim_data[11][t48,index_gsk3_p]]
 BLOT = [Fig4F[1],
         Fig4F[2],
         Fig4F[3]]
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[7] = SSE(MOD,BLOT)# SSE(MOD,BLOT)
 #@show MOD
 #@show BLOT,MOD #SHOWCODEOBJ
 # println(MOD,' ',BLOT,' ',SSE(MOD,BLOT))
 # if MOD == [1.0,0.0]
 #   @bp
 # end
 #8 EXP TGFB3 DNSMAD RO LEF # Fig4E
 MOD  = [sim_data[1][t48,index_lef]+sim_data[1][t48,index_lef_act],
         sim_data[8][t48,index_lef]+sim_data[8][t48,index_lef_act],
         sim_data[9][t48,index_lef]+sim_data[9][t48,index_lef_act]]
 BLOT = [Fig4E[1],
         Fig4E[2],
         Fig4E[3]]
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[8]  = SSE(MOD,BLOT)
#@show MOD
#@show BLOT,MOD #SHOWCODEOBJ
 #9 EXP TGFB3 DNSMAD DNLEF RO ECAD VIM # Fig4G #  index_ecad_bcat

 MOD  = [sim_data[1][t48,index_ecad]+sim_data[1][t48,index_ecad_bcat],sim_data[1][t48,index_vim],
         sim_data[8][t48,index_ecad]+sim_data[8][t48,index_ecad_bcat],sim_data[8][t48,index_vim],
         sim_data[9][t48,index_ecad]+sim_data[9][t48,index_ecad_bcat],sim_data[9][t48,index_vim],
         sim_data[10][t48,index_ecad]+sim_data[10][t48,index_ecad_bcat],sim_data[10][t48,index_vim]]
 BLOT = [Fig4G[1],Fig4G[5],
         Fig4G[2],Fig4G[6],
         Fig4G[3],Fig4G[7],
         Fig4G[4],Fig4G[8]]
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[9]  = SSE(MOD,BLOT)
#@show MOD
#@show BLOT,MOD #SHOWCODEOBJ
 #10 EXP TGFB2 DNMEK RO SS # Fig5A
 MOD  = [sim_data[1][t48,index_ss],
         sim_data[2][t48,index_ss],
         sim_data[3][t48,index_ss]]
 BLOT = [Fig5A[1],
         Fig5A[2],
         Fig5A[3]]
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[10] = SSE(MOD,BLOT)
#@show MOD
#@show BLOT,MOD #SHOWCODEOBJ
 #11 EXP TGFB TGFB3_Abn RO LEF1 # Fig5C
  # generate new experiment EXP TGFB + TGFB3_Ab RO: LEF1 (sim 11
 # Retain old objective using simulation 3
 MOD  = [sim_data[1][t48,index_lef]+sim_data[1][t48,index_lef_act],
         sim_data[2][t48,index_lef]+sim_data[2][t48,index_lef_act],
         sim_data[3][t48,index_lef]+sim_data[3][t48,index_lef_act]] # patched
 BLOT = [Fig5C[1],
         Fig5C[2],
         Fig5C[3]]
 MOD = norm_dat(MOD)
 BLOT = norm_dat(BLOT)
 obj_array[11] = SSE(MOD,BLOT) # OBJECTIVE DISABLED -- no effect...  ?
#@show MOD
#@show BLOT,MOD #SHOWCODEOBJ
 # Constraints are implemented as a penaltly on obj value
 lambda_value = 100.0

 # How much do we violate the constraints?
 #violation_constraint_1 = 25 - (x-5.0)^2 - y^2
 #violation_constraint_2 = (x-8.0)^2 + (y-3.0)^2 - 7.7
 penaltly_array = zeros(num_objs)

 #penaltly_array[1] = lambda_value*(min(0,violation_constraint_1))^2
 #penaltly_array[2] = lambda_value*(min(0,violation_constraint_2))^2

# For quick Cross validation, simply set obj_array[x] to 0

#obj_array[]

# pop objective from obj_array for cross val
 # return the obj_array -
 #@show obj_array
 #stop()

# apply cross validation hold out (expensive but easy way)
 if HOLD!=0
 obj_array[HOLD]=0.1 # ignore evlautation of this objective
 end
 #f@show obj_array
 #@show obj_array
 return obj_array+penaltly_array
end

#     #    #   # Provide other model specific functions for POETS
# FROM JV ...
# Generates new parameter array, given current array -
function neighbor_function(parameter_array)

  SIGMA = 0.05
  number_of_parameters = length(parameter_array)

  # calculate new parameters -
  new_parameter_array = parameter_array.*(1+SIGMA*randn(number_of_parameters))

  # Check the bound constraints -
  LOWER_BOUND = zeros(parameter_array) #+ 0.1
  UPPER_BOUND = zeros(parameter_array)
  #UPPER_BOUND = LOWER_BOUND #+ 10
  # we need a different set of bounds for the control functions
  # alpha_lower,alpha_upper = 0.1,10.0
   gain_lower,gain_upper = 0.001, 10.0
   order_lower,order_upper = 0.5,2.5
  #
  # LOWER_BOUND[1:len_rca]=LOWER_BOUND[1:len_rca]+alpha_lower
  # LOWER_BOUND[len_rca+1:A]=LOWER_BOUND[len_rca+1:A]+gain_lower
  # LOWER_BOUND[A+1:B]=LOWER_BOUND[A+1:B]+order_lower
  #
  # UPPER_BOUND[1:len_rca]=UPPER_BOUND[1:len_rca]+alpha_upper

  #
  deviation = 0.3 # 0.5
  LOWER_BOUND=original_parameter_array*(1-deviation)
  UPPER_BOUND=original_parameter_array*(1+deviation)

  # UPPER_BOUND[len_rca+1:A]=zeros(parameter_array)[len_rca+1:A]+gain_upper
  # UPPER_BOUND[A+1:B]=zeros(parameter_array)[A+1:B]+order_upper
  # LOWER_BOUND[len_rca+1:A]=zeros(parameter_array)[len_rca+1:A]+gain_lower
  # LOWER_BOUND[A+1:B]=zeros(parameter_array)[A+1:B]+order_lower

  # deviation = 10.0
  # LOWER_BOUND=original_parameter_array*(deviation)
  # UPPER_BOUND=original_parameter_array*(1./deviation)

  #println(parameter_array)
  #@show parameter_array
  # return the corrected parameter arrays -
  return parameter_bounds_function(new_parameter_array,LOWER_BOUND,UPPER_BOUND)
end

function acceptance_probability_function(rank_array,temperature)

  return (exp(-rank_array[end]/temperature))

end

function cooling_function(temperature)

  # define my new temperature -
  alpha = 0.8 #0.9   # testing 0.5
  return alpha*temperature
end

# Helper functions -
function parameter_bounds_function(parameter_array,lower_bound_array,upper_bound_array)

  # reflection_factor -
  epsilon = 0.01

  # iterate through and fix the parameters -
  new_parameter_array = copy(parameter_array)
  for (index,value) in enumerate(parameter_array)

# only allow parameters to vary slightly outside of thier original range
# See if the system remains stable

    lower_bound = lower_bound_array[index]
    upper_bound = upper_bound_array[index]
    if (value<lower_bound)
      new_parameter_array[index] = lower_bound+epsilon*upper_bound
    elseif (value>upper_bound)
      new_parameter_array[index] = upper_bound - epsilon*lower_bound
    end
  end

  # fix cpa values relating to addition terms
  # GAINS 1 2 10 19
  new_parameter_array[len_rca+1]=30.0# GAIN_tgfb3_Ab
  new_parameter_array[len_rca+2]=0.3# GAIN_vivt
  new_parameter_array[len_rca+10]=5.0# GAIN_ecad_virus
  new_parameter_array[len_rca+19]=5.0# GAIN_SS_virus
  # ORDERS
  new_parameter_array[len_rca+len_cpa_gain+1]=1.0# GAIN_tgfb3_Ab
  new_parameter_array[len_rca+len_cpa_gain+2]=1.0# GAIN_vivt
  new_parameter_array[len_rca+len_cpa_gain+10]=1.0# GAIN_ecad_virus
  new_parameter_array[len_rca+len_cpa_gain+19]=1.0# GAIN_SS_virus

   # machinery
   # these are set to zero in the balances function
  # new_parameter_array[97]=pa97
  # new_parameter_array[98]=pa98
  # new_parameter_array[99]=pa99
  # new_parameter_array[100]=pa100
  # new_parameter_array[101]=pa101
  # new_parameter_array[102]=pa102

  # addition terms
 new_parameter_array[85:100] = addition_terms_saved

#  parameter_array[117]=pa117
#  parameter_array[118]=pa118

  return new_parameter_array
end

function n_alpha(alpha,Tf=0.0001;Ti=1.0)
Ti = 1.0
(log(Tf)-log(Ti))/log(alpha)
end


text:
"""
perfect_square(number)

Check if a number is a perfect square or not.

# Example
```jula
perfect_square(9)   # returns True
perfect_square(16)  # returns True
perfect_square(1)   # returns True
perfect_square(0)   # returns True
perfect_square(10)  # returns False
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function perfect_square(number)
    return sqrt(number) * sqrt(number) == number
end
perfect_square(9)

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
module PAPI

@osx_only error("PAPI.jl currently only works on Linux")
@windows_only error("PAPI.jl currently only works on Linux")

include("events.jl")
include("retcodes.jl")

immutable PAPIError{R} <: Exception
    msg::String
end
PAPIError(R::RetCode) = PAPIError{R}(errmsg(R))

function __init__()
    # init the library and make sure that some counters are available
    if num_counters() <= 0
        error("PAPI init error: No counters are available on the current system")
    end
    atexit() do
        ccall((:PAPI_shutdown, :libpapi), Void, ())
    end
end

type EventSet
    counters::Vector{Event}
    vals::Vector{Clonglong}

    EventSet(c::Vector{Event}) = begin
        cs = new(c, zeros(Clonglong, length(c)))
        return cs
    end
end

# a nonexistent hardware event used as a placeholder
const PAPI_NULL = Cint(-1)

immutable PAPIEventSet
    val::Cint
    PAPIEventSet() = new(PAPI_NULL)
end

# TODO Low Level API
function exists(evt::Event)
    evtset = Ref(PAPIEventSet())
    ret = RetCode(ccall((:PAPI_create_eventset, :libpapi), Cint,
                        (Ref{PAPIEventSet},), evtset))
    if ret != OK
        throw(PAPIError(ret))
    end
    ret = RetCode(ccall((:PAPI_query_event, :libpapi), Cint,
                        (Cuint,), evt))
    return ret == OK
end

function start(evts::Event...)
    nevts = length(evts)
    if nevts == 0
        throw(ArgumentError("one or more PAPI.Events required"))
    end
    ncounters = num_counters()
    if nevts > ncounters
        throw(ArgumentError("number of PAPI.Events must be ≤ PAPI.num_counters(), got $nevts"))
    end
    evtcodes = Array(Cuint, nevts)
    for i = 1:length(nevts)
        evtcodes[i] = Cuint(evts[i])
    end
    ret = RetCode(ccall((:PAPI_start_counters, :libpapi), Cint,
                        (Ptr{Cuint}, Cint), evtcodes, nevts))
    if ret != OK
        throw(PAPIError(ret))
    end
    return
end

#### High Level Interface ####

@doc """
""" ->
is_initialized() = Bool(ccall((:PAPI_is_initialized, :libpapi), Cint, ()))

@doc """
Get the number of hardware counters available on the system

`PAPI.num_counters()` returns the optimal length of the values array for high-level functions.
This value corresponds to the number of hardware counters supported by the current substrate.
`PAPI.num_counters()` initializes the PAPI library using `PAPI.library_init()` if necessary.
""" ->
num_counters() = Int(ccall((:PAPI_num_counters, :libpapi), Cint, ()))

@doc """
Add current counts to array and reset counters
""" ->
function accum_counters!(values::Vector{Clonglong})
    ret = RetCode(ccall((:PAPI_accum_counters, :libpapi), Cint,
                        (Ptr{Clonglong}, Cint), values, length(values)))
    if ret != OK
        throw(PAPIError(ret))
    end
end
accum_counters!(cs::EventSet) = (accum_counters!(cs.vals); copy(cs.vals))

@doc """
Get the number of components available on the system
""" ->
num_components() = Int(ccall((:PAPI_num_components, :libpapi), Cint, ()))

@doc """
""" ->
function read_counters!(values::Vector{Clonglong})
    ret = RetCode(ccall((:PAPI_read_counters, :libpapi), Cint,
                        (Ptr{Clonglong}, Cint), values, length(values)))
    if ret != OK
        throw(PAPIError(ret))
    end
end
read_counters!(cs::EventSet) = (read_counters!(cs.vals); copy(cs.vals))

@doc """
Start counting hardware events

`PAPI.start_counters()` initializes the PAPI library (if necessary) and starts counting the events named in the events array.
This function implicitly stops and initializes any counters running as a result of a previous call to `PAPI.start_counters()`.
It is the user's responsibility to choose events that can be counted simultaneously by reading the vendor's documentation.
The number of events should be no larger than the value returned by `PAPI.num_counters()`.
""" ->
function start_counters(events::Vector{Event})
    ret = RetCode(ccall((:PAPI_start_counters, :libpapi), Cint,
                        (Ptr{Cint}, Cint), pointer(events), length(events)))
    if ret != OK
        throw(PAPIError(ret))
    end
end
start_counters(cs::EventSet) = start_counters(cs.counters)

@doc """
Stop counters and return current counts
""" ->
function stop_counters()
    ret = RetCode(ccall((:PAPI_stop_counters, :libpapi), Cint,
                        (Ptr{Cuint}, Cint), C_NULL, 0))
    if ret != OK
        throw(PAPIError(ret))
    end
end
function stop_counters(events::Vector{Event})
    ret = RetCode(ccall((:PAPI_stop_counters, :libpapi), Cint,
                        (Ptr{Cuint}, Cint), pointer(events), length(events)))
    if ret != OK
        throw(PAPIError(ret))
    end
end
stop_counters(cs::EventSet) = stop_counters(cs.counters)

const _rtime = Cfloat[0.0]
const _ptime = Cfloat[0.0]
const _flpx = Clonglong[0]
const _mflpx = Cfloat[0.0]

@doc """
Get Mflips/s (floating point instruction rate), real time and processor time
""" ->
function flips!(flpins, mflips, rtime, ptime)
    ret = RetCode(ccall((:PAPI_flips, :libpapi), Cint,
                        (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Clonglong}, Ptr{Cfloat}),
                        rtime, ptime, flpins, mflips))
    if ret != OK
        throw(PAPIError(ret))
    end
end

@doc """
Get Mflips/s (floating point instruction rate), real time and processor time
""" ->
function flips()
    flips!(_flpx, _mflpx, _rtime, _ptime)
    return (_flpx[1], _mflpx[1], _rtime[1], _ptime[1])
end

macro flips(blk)
    quote
        PAPI.flips()
        $(esc(blk))
        local res = PAPI.flips()
        PAPI.stop_counters([PAPI.FP_INS])
        ret
    end
end

@doc """
Get Mflop/s (floating point operand rate), real time and processor time
""" ->
function flops!(flpops, mflops, rtime, ptime)
    ret = RetCode(ccall((:PAPI_flops, :libpapi), Cint,
                        (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Clonglong}, Ptr{Cfloat}),
                        rtime, ptime, flpops, mflops))
    if ret != OK
        throw(PAPIError(ret))
    end
end

@doc """
Get Mflop/s (floating point operand rate), real time and processor time
""" ->
function flops()
    flops!(_flpx, _mflpx, _rtime, _ptime)
    return (_flpx[1], _mflpx[1], _rtime[1], _ptime[1])
end

macro flops(blk)
    quote
        PAPI.flops()
        $(esc(blk))
        local res = PAPI.flops()
        PAPI.stop_counters([PAPI.FP_OPS])
        res
    end
end

@doc """
Get instructions per cycle, real time and processor time
""" ->
function ipc(rtime, ptime, ins, ipc)
    ret = RetCode(ccall((:PAPI_ipc, :libpapi), Cint,
                        (Ptr{Void}, Ptr{Cfloat}, Ptr{Clonglong}, Ptr{Cfloat}),
                        rtime, ptime, ins, ipc))
    if ret != OK
        throw(PAPIError(ret))
    end
end

@doc """
Get events per cycle, real time and processor time
""" ->
function epc(event, rtime, ptime, ref, core, evt, epc)
    ret = RetCode(ccall((:PAPI_epc, :libpapi), Cint,
                        (Cint, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Clonglong},
                         Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Cfloat}),
                        event, rtime, ptime, ref, core, evt, epc))
    if ret != OK
        throw(PAPIError(ret))
    end
end

end # module


text:
"""
    sdpa_gmp_binary(arg, silent::Bool = false)
Execute the given command literal as the only argument to SDPA-GMP.

"""
function  sdpa_gmp_binary(arg::Cmd, silent = false)
    if !silent
        run(`$sdpa_gmp_path $arg`)
    else
        out = devnull
        run(pipeline(`$sdpa_gmp_path $arg`, stdout = out))
    end
end

"""
    sdpa_gmp_binary_solve!(m::Optimizer, full_input_path, full_output_path, extra_args::Cmd, redundant_entries)

Calls the binary `sdpa_gmp` to solve SDPA-formatted problem specified in a .dat-s file at `full_input_path`. Results are written into the file at `full_output_path`.
`extra_args` is passed on to the binary as additional options, allowing for e.g. custom parameter files. `redundant_entries` is a sorted vector listing indices of linearly dependent constraint which are already removed by `presolve.jl`. The corresponding decision variables are populated as zeros.

This function returns `m` with solutions already populated from results in the output file.
"""
function sdpa_gmp_binary_solve!(m::Optimizer, full_input_path::String, full_output_path::String; extra_args::Cmd = `-pt 0`, redundant_entries::Vector = [])
    # extra_args = `-p /home/jiazheng/Downloads/sdpa-gmp-7.1.3/param.sdpa`
    sdpa_gmp_binary(`-ds $full_input_path -o $full_output_path $extra_args`, m.silent);
    read_results!(m, full_output_path, redundant_entries);
    return m
end


text:
using Millboard
using Base.Test

board = []
@test """
+==+""" == table(board) |> string

board = [11 12 13; 21 22 23]
@test """
+===+====+====+====+
|   |  1 |  2 |  3 |
+===+====+====+====+
| 1 | 11 | 12 | 13 |
+---+----+----+----+
| 2 | 21 | 22 | 23 |
+---+----+----+----+""" == table(board) |> string

board = ([1 2], [5 6;7 8], [9 10; 11 12])
@test """
+===+=====+=====+=======+
|   |   1 |   2 |     3 |
+===+=====+=====+=======+
|   | 1 2 | 5 6 |  9 10 |
| 1 |     | 7 8 | 11 12 |
+---+-----+-----+-------+""" == table(board) |> string

@test """
+========+=====+=====+=======+
|        |   x |   y |     z |
+========+=====+=====+=======+
|        | 1 2 | 5 6 |  9 10 |
| result |     | 7 8 | 11 12 |
+--------+-----+-----+-------+""" ==
  table(board, :colnames=>["x" "y" "z"], :rownames=>["result"]) |> string

board = ([1 2], [5 6;7 8], [9 10; 11 12])
@test """
+===+=====+=====+=======+
|   |   1 |   2 |     3 |
+===+=====+=====+=======+
|   | 1 2 | 5 6 |  9 10 |
| 1 |     | 7 8 | 11 12 |
+---+-----+-----+-------+""" == table(board) |> string

board = ([], [5 6;7 8], [9 10; 11 12])
@test """
+===+===+=====+=======+
|   | 1 |   2 |     3 |
+===+===+=====+=======+
|   |   | 5 6 |  9 10 |
| 1 |   | 7 8 | 11 12 |
+---+---+-----+-------+""" == table(board) |> string

board = ([1 2;3 4], [5 6;7 8], [9 10; 11 12])
@test """
+===+=====+=====+=======+
|   |   1 |   2 |     3 |
+===+=====+=====+=======+
|   | 1 2 | 5 6 |  9 10 |
| 1 | 3 4 | 7 8 | 11 12 |
+---+-----+-----+-------+""" == table(board) |> string


text:
module SpringCollab2020VariableCrumbleBlock

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/variableCrumbleBlock" VariableCrumbleBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, texture::String="default",
    timer::Number=0.4, respawnTimer::Number=2.0)

const placements = Ahorn.PlacementDict(
    "Variable Crumble Blocks ($(uppercasefirst(texture))) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        VariableCrumbleBlock,
        "rectangle",
        Dict{String, Any}(
            "texture" => texture
        )
    ) for texture in Maple.crumble_block_textures
)

Ahorn.editingOptions(entity::VariableCrumbleBlock) = Dict{String, Any}(
    "texture" => Maple.crumble_block_textures
)

Ahorn.minimumSize(entity::VariableCrumbleBlock) = 8, 0
Ahorn.resizable(entity::VariableCrumbleBlock) = true, false

function Ahorn.selection(entity::VariableCrumbleBlock)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::VariableCrumbleBlock, room::Maple.Room)
    texture = get(entity.data, "texture", "default")
    texture = "objects/crumbleBlock/$texture"

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))
    tilesWidth = div(width, 8)

    Ahorn.Cairo.save(ctx)

    Ahorn.rectangle(ctx, 0, 0, width, 8)
    Ahorn.clip(ctx)

    for i in 0:ceil(Int, tilesWidth / 4)
        Ahorn.drawImage(ctx, texture, 32 * i, 0, 0, 0, 32, 8)
    end

    Ahorn.restore(ctx)
end

end


text:
using ImageQualityIndexes
using Test, ReferenceTests, TestImages
using Statistics

include("testutils.jl")

@testset "ImageQualityIndexes" begin

    include("psnr.jl")
    include("ssim.jl")
    include("colorfulness.jl")

end

nothing


text:
import JLD2, FileIO
using Serialization, Test, HDF5

for N in [1_000, 2_000, 3_000, 4_000]

A = rand(N,200_000)
sizeof(A) / 1e9

@info "N = $N"
@info "sizeof(A) = $(round(sizeof(A) / 1e9, 2))GB"
@info "Testing JLD2"
@time FileIO.save("temp.jld2", "A", A)
@time B = FileIO.load("temp.jld2", "A")
println(@test B == A)
run(`rm temp.jld2`)

@info "Testing serialize/deserialize"
function mysave(A, fname)
   f = open(fname, "w")
   serialize(f, A)
   close(f)
end
function myload(fname)
   f = open(fname, "r")
   A = deserialize(f)
   close(f)
   return A
end
@time mysave(A, "temp.dat")
@time B = myload("temp.dat")
run(`rm temp.dat`)
println(@test A == B)

@info "Testing HDF5"
function savehdf5(A, fname)
   h5open(fname, "w") do fid
      fid["A"] = A
   end
end
function loadhdf5(fname)
   fid = h5open(fname, "r")
   A = read(fid["A"])
   close(fid)
   return A
end
@time savehdf5(A, "temp.h5")
@time B = loadhdf5("temp.h5")
run(`rm temp.h5`)
println(@test A == B)

end


text:
## Julia 1.1.0

z = 1 + 2im
for T in [real, imag, conj, abs, abs2, angle, sqrt, sin, cos, sinh, cosh, exp]
    println("$(T)($(z)) -> $(T(z))")
end

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
ti = time()

# lyapunov Exponents:
include("lyapunov_exponents.jl")
# Chaos Detection:
include("chaos_detection_tests.jl")
# Orbit Diagrams:
include("orbitdiagram_tests.jl")
# Periodicity:
include("periodicity_tests.jl")
# Entropies (and attractor dimensions)
include("entropy_dimension.jl")
# Nonlinear Timeseries Analysis:
include("nlts_tests.jl")
include("R_params.jl")

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits = 33), " seconds or ",
round(ti/60, digits = 3), " minutes")


text:
function add_solver_xml!(::NEOSSolver{:Xpress, :MPS}, m::MPSModel)
	# Add solution display
	m.xmlmodel = replace(m.xmlmodel, r"(?s)<algorithm>.*</algorithm>", "<algorithm><![CDATA[SIMPLEX]]></algorithm>")

	# Add user options
	param_string = ""
	for key in keys(m.solver.params)
		param_string *= "$(key) \= $(m.solver.params[key])\n"
	end
	# Add user options
	m.xmlmodel = replace(m.xmlmodel, r"(?s)<par>.*</par>", "<par><![CDATA[$(param_string)]]></par>")
end

# Objective function value is     4.500000
#	...
#   Number   Column   At      Value      Input Cost   Reduced Cost
# C      4  VAR1      SB      1.000000      1.000000       .000000
# C      5  VAR2      SB      4.500000      1.000000       .000000
# C      6  VAR3      SB      1.000000     -1.000000       .000000

function parse_status!(::NEOSSolver{:Xpress, :MPS}, m::MPSModel)
	if contains(m.last_results, "Optimal solution found")
		m.status = OPTIMAL
	elseif contains(m.last_results, "problem is unbounded")
		m.status = UNBOUNDED
	elseif contains(m.last_results, "infeasible")
		m.status = INFEASIBLE
	end
end

function parse_objective!(::NEOSSolver{:Xpress, :MPS}, m::MPSModel)
	m.objVal = parse(Float64, match(r"Objective function value is\W+?(-?[\d.]+)", m.last_results).captures[1])
end

function parse_solution!(::NEOSSolver{:Xpress, :MPS}, m::MPSModel)
	for v in matchall(r"V(\d+).+?(-?[\d.]+)", m.last_results)
		regmatch = match(r"V(\d+).+?(-?[\d.]+)", v)
		m.solution[parse(Int, regmatch.captures[1])] = parse(Float64, regmatch.captures[2])
	end
end


text:
"""
    ProjectiveTransforms(sz; [augmentations, buffered]) <: Encoding

Encoding for spatial data that resizes blocks to a common size `sz` and applies
projective augmentations.

Encodes all spatial blocks, preserving the block type:
- `Image{N}` -> `Image{N}`
- `Mask{N}` -> `Mask{N}`
- `Keypoints{N}` -> `Keypoints{N}`

The behavior differs based on the `context` of encoding:

{.tight}
- [`DLPipelines.Training`](#):
    1. Resizes the data so the smallest side equals
    a side length in `sz` while keeping the aspect ratio.
    2. Applies `augmentations`.
    3. Crops a random `sz`-sized portion of the data
- [`DLPipelines.Validation`](#):
    1. Resizes the data so the smallest side equals
    a side length in `sz` while keeping the aspect ratio.
    2. Crops a `sz`-sized portion from the center.
- [`DLPipelines.Inference`](#):
    1. Resizes the data so the smallest side equals
    a side length in `sz` while keeping the aspect ratio.
        Note that in this context, the data does not have
        size `sz`, since no cropping happens and aspect ratio
        is preserved.

`ProjectiveTransforms` is not limited to 2D data, and works on 3D data as well.
Note, however, that some transformations in `augs_projection` (rotation, warping, flipping)
are 2D only so `augs_projection` cannot be used for 3D data.

## Keyword arguments

- `augmentations::`[`DataAugmentation.Transform`](#)` = Identity()`: Projective
    augmentation to apply during training. See [`augs_projection`](#).
- `buffered = true`: Whether to use inplace transformations. Reduces memory usage.
- `sharestate = true`: Whether to use the same random state and bounds for all blocks
    in a sample

"""
@with_kw_noshow struct ProjectiveTransforms{N} <: StatefulEncoding
    sz::NTuple{N, Int}
    buffered::Bool
    augmentations
    tfms::Dict{Context, DataAugmentation.Transform}
    sharestate::Bool
end

function ProjectiveTransforms(
        sz;
        augmentations = Identity(),
        inferencefactor = 8,
        buffered = true,
        sharestate = true)

    tfms = Dict{Context, DataAugmentation.Transform}(
        Training() => ScaleKeepAspect(sz) |> augmentations |> RandomCrop(sz) |> PinOrigin(),
        Validation() => CenterResizeCrop(sz),
        Inference() => ResizePadDivisible(sz, inferencefactor),
    )
    if buffered
        tfms[Training()] = BufferedThreadsafe(tfms[Training()])
        tfms[Validation()] = BufferedThreadsafe(tfms[Validation()])
    end

    return ProjectiveTransforms(sz, buffered, augmentations, tfms, sharestate)
end

function encodestate(enc::ProjectiveTransforms{N}, context, blocks, datas) where N
    bounds = getsamplebounds(blocks, datas, N)
    tfm = enc.tfms[context]
    randstate = DataAugmentation.getrandstate(tfm)
    return bounds, randstate
end

function encodedblock(enc::ProjectiveTransforms{N}, block::Block) where N
    return isnothing(blockitemtype(block, N)) ? nothing : block
end

function encode(
        enc::ProjectiveTransforms{N},
        context,
        block::Block,
        data;
        state=nothing) where N
    ItemType = blockitemtype(block, N)
    isnothing(ItemType) && return data
    # only init state if block is encoded
    bounds, randstate = (isnothing(state) || !enc.sharestate) ? encodestate(enc, context, block, data) : state
    # don't encode if bounds have wrong dimensionality
    bounds isa DataAugmentation.Bounds{N} || return data

    tfm = enc.tfms[context]
    item = ItemType(data, bounds)
    tdata = apply(tfm, item; randstate = randstate) |> itemdata
    return copy(tdata)
end

# ProjectiveTransforms is not invertible, hence no `decode` method!

# Conversion of `Block` to `DataAugmentation.Item`

"""
    blockitemtype(block, N)

Return a constructor for a `DataAugmentation.Item` that can be projected.
Return `nothing` by default, indicating that `block` cannot be turned into a
projectable item for bounds with dimensionality `N`.
For example, we have

    blockitemtype(Image{2}(), 2) -> DataAugmentation.Image

but

    blockitemtype(Image{3}(), 2) -> nothing
"""
blockitemtype(block::Block, n) = nothing
blockitemtype(block::Image{N}, n::Int) where N = N == n ? DataAugmentation.Image : nothing
function blockitemtype(block::Mask{N}, n::Int) where N
    return if N == n
        (data, bounds) -> DataAugmentation.MaskMulti(data, block.classes, bounds)
    else
        nothing
    end
end
blockitemtype(block::Keypoints{N}, n::Int) where N = N == n ? DataAugmentation.Keypoints : nothing
blockitemtype(block::WrapperBlock, n::Int) = blockitemtype(wrapped(block), n)

"""
    grabbounds(blocks, datas, N)

Looks through `blocks` for block data that carries `N`-dimensional
bounds information needed for projective transformations.
"""
function grabbounds(blocks::Tuple, datas::Tuple, N::Int)
    for (block, data) in zip(blocks, datas)
        bounds = grabbounds(block, data, N)
        !isnothing(bounds) && return bounds
    end
end

grabbounds(block::Image{N}, a, n) where N = N == n ? DataAugmentation.Bounds(size(a)) : nothing
grabbounds(block::Mask{N}, a, n) where N = N == n ? DataAugmentation.Bounds(size(a)) : nothing
grabbounds(block::WrapperBlock, a, n) = grabbounds(wrapped(block), a, n)

function getsamplebounds(blocks, datas, N)
    bounds = grabbounds(blocks, datas, N)
    isnothing(bounds) && error("Could not detect $N-dimensional bounds needed for projective
transformations from blocks $(blocks)! Bounds can be grabbed from arrays
like `Image{$N}` block data.")
    return bounds
end

# Augmentation helper

"""
    augs_projection([; kwargs...])

Helper to create a set of projective transformations for image, mask
and keypoint data. Similar to fastai's
[`aug_transforms`](https://github.com/fastai/fastai/blob/bdc58846753c6938c63344fcaebea7149585fd5c/fastai/vision/augment.py#L946).

## Keyword arguments

- `flipx = true`: Whether to perform a horizontal flip with probability `1/2`. See [`FlipX`](#).
- `flipy = false`: Whether to perform a vertical flip with probability `1/2`. See [`FlipY`](#).
- `max_zoom = 1.5`: Maximum factor by which to zoom. Set to `1.` to disable. See [`Zoom`](#).
- `max_rotate = 10`: Maximum absolute degree by which to rotate. Set to `0.` to disable. See [`Rotate`](#).
- `max_warp = 0.05`: Intensity of corner warp. Set to `0.` to disable. See [`WarpAffine`](#).
"""
function augs_projection(;
        flipx=true,
        flipy=false,
        max_zoom=1.5,
        max_rotate=10.,
        max_warp=0.05,
        )
    tfms = []

    flipx && push!(tfms, Maybe(FlipX()))
    flipy && push!(tfms, Maybe(FlipY()))
    max_warp > 0 && push!(tfms, WarpAffine(max_warp))
    max_rotate > 0 && push!(tfms, Rotate(max_rotate))
    push!(tfms, Zoom((1., max_zoom)))
    return DataAugmentation.compose(tfms...)
end

# Pretty-printing

function Base.show(io::IO, p::ProjectiveTransforms)
    show(io, ShowTypeOf(p))
    fields = (
        sz = ShowLimit(p.sz, limit=80),
        buffered = ShowLimit(p.buffered, limit=80),
        augmentations = ShowLimit(p.augmentations, limit=80)
    )
    show(io, ShowProps(fields, new_lines=true))
end


text:
abstract type ElectricLoad <: Injection end

include("loads/electric_loads.jl")
include("loads/controllable_loads.jl")
include("loads/shunt_elements.jl")


text:
module ADMIXTURE

using DelimitedFiles, Glob, Pkg.Artifacts

export ADMIXTURE_EXE
export admixture, read_PQ

if Sys.isapple()
    ADMIXTURE_EXE = joinpath(artifact"ADMIXTURE", "dist/admixture_macosx-1.3.0/admixture")
elseif Sys.islinux()
    ADMIXTURE_EXE = joinpath(artifact"ADMIXTURE", "dist/admixture_linux-1.3.0/admixture")
else
    error("ADMIXTURE only supports Linux and MacOS")
end

"""
    admixture(inputFile, K; kwargs...)

Calculate the admixture estimates assuming `K` populations. Input file can be
`.bed` (together with `.bim`, binary PLINK format), `geno` (STRUCTURE format),
or `ped` (12-coding PLINK format). The keyword arguments match those of
ADMIXTURE software. For details, see the ADMIXTUR manual
<http://dalexander.github.io/admixture/admixture-manual.pdf>

# Keyword arguments
- `B`: compute standard errors by bootstrap; default is `nothing`.
- `j`: number of threads to be used; default is `nothing`.
- `s`: random seed; integer or `:time`; default is `nothing`.
- `C`: major convergence criterion; float or integer; default is `nothing`.
- `c`: minor convergence criterion; float or integer; default is `nothing`.
- `qn`: secant pairs to be used in quasi-Newton acceleration; default is 3.
- `l`: tuning parameter in penalized estimation; default is `nothing`.
- `e`: tuning parameter in penalized estimation; default is `nothing`.
- `supervised`: use supervised training (requires `.pop` file); default is false
"""
function admixture(
    inputFile :: AbstractString,
    K         :: Integer;
    B         :: Union{Nothing,Integer} = nothing,
    j         :: Union{Nothing,Integer} = nothing,
    cv        :: Union{Nothing,Integer} = nothing,
    s         :: Union{Nothing,Integer,Symbol} = nothing,
    C         :: Union{Nothing,Real} = nothing,
    c         :: Union{Nothing,Real} = nothing,
    qn        :: Union{Nothing,Integer} = nothing,
    l         :: Union{Nothing,Number} = nothing,
    e         :: Union{Nothing,Number} = nothing,
    supervised::Bool = false
    )
    inputPath = inputFile |> abspath |> normpath
    args  = sizehint!(String[], 12)
    isnothing(B)  || push!(args, isinteger(B) ? "-B$B" : "-B")
    isnothing(j)  || push!(args, "-j$j")
    isnothing(cv) || push!(args, "--cv=$cv")
    isnothing(s)  || push!(args, isinteger(s) ? "-s $s" : "-s time")
    isnothing(C)  || push!(args, "-C $C")
    isnothing(c)  || push!(args, "-c $c")
    isnothing(qn) || push!(args, qn == 0 ? "-a none" : "-a qn$qn")
    isnothing(l)  || isnothing(e) || push!(args, "-l $l", "-e $e")
    !supervised   || push!(args, "--supervised")
    cmd  = `$ADMIXTURE_EXE $inputPath $K $args`
    @info "ADMIXTURE command:\n$cmd\n"
    @info "Output directory: $(pwd())\n"
    run(cmd)
    outputPath = joinpath(pwd(), splitpath(inputPath)[end][1:end-4] * ".$K")
    P, Q = read_PQ(outputPath, args=args)
end

"""
    read_PQ(outputFile)

Read the P (allele frequencies of the inferred ancestral populations) and
Q (the ancestry fractions) estimates output by the ADMIXTURE program.

# Keyword arguments
- `T`:: output data type; default is `Float64`.
- `args`:: parsed ADMIXTURE arguments as a `Vector{String}`; default is empty.
"""
function read_PQ(
    outputFile :: AbstractString;
    T          :: Type = Float64,
    args       :: Vector{String} = String[]
    )
    pfile = normpath(abspath(outputFile)) * ".P"
    qfile = normpath(abspath(outputFile)) * ".Q"
    if any(arg -> occursin("-l ", arg), args)
        # penalized estimates have different file name endings
        # assume unique file ending with P.lambda.* and Q.lambda.*
        pfile = glob("$(basename(pfile)).lambda.*", pwd())[1]
        qfile = glob("$(basename(qfile)).lambda.*", pwd())[1]
    end
    P = readdlm(pfile, ' ', T)
    Q = readdlm(qfile, ' ', T)
    P, Q
end

end # module


text:
using Documenter
using Random
using CorrelationFunctions
using CorrelationTrackers

makedocs(sitename = "CorrelationTrackers.jl documentation",
         format = Documenter.HTML(prettyurls = false))

deploydocs(repo = "github.com/shamazmazum/CorrelationTrackers.jl.git")


text:
#=
Types and methods for tensors (multidimension array).

It is cruft if it says "image" instead of "tensor"
=#

include("null.jl")

#=
A tensor (usually 2D image) and same length bit mask.

The mask is often a "selection".
That defines a non-rectangular region.
More generally, a mask is another channel of an image,
for example, defining transparency.
More generally a mask can be a byte, with degrees of masking,
e.g. for selection, defining a soft-edged or dithered selection.

Parameterized.
E.G. For usual notion of "2D image": T is Color, and N is 2.
E.G. For usual notion of "3D image": T is Color, and N is 3.
E.G. For usual notion of "text": T is Char, and N is 1.
E.G. For unusual notion of "poem": T is Char, and N is 2.

!!! Not mutable (so it is not allocated),
Only the image and mask themselves are mutable.
But the algorithm does not always mutate the image (on the target, not the corpus)
and never mutates the mask.
=#

struct MaskedImage{ValueType, DimensionCount}

    #=
    NOT AbstractArray to allow more implementations of Array.
    Here we choose a concrete type appropriate to the algorithm.
    =#
    # WAS image::AbstractArray{ValueType, DimensionCount}
    # WAS mask::AbstractArray{Bool, DimensionCount}
    image::Array{ValueType, DimensionCount}
    mask::BitArray{DimensionCount}

    #=
    inner constructor.
    Enforces invariant that mask is same dimensions as image.
    =#
    #=
    temp
    function MaskedImage(image, mask)
        @assert size(image) == size(mask) "image and mask must have same size"
        return new{typeof(image), size(image)}(image, mask)
    end
    =#
end

#=
Return CartesionIndex with each element a rand in range of indices of each dimension of array.

!!! Not unique, may repeat.

General on N i.e. DimensionCount.
We parameterize the function declaration on N.

The returned value is a CartesianIndex, or for 1D a Int64??? TODO

??? Performance.
=#
function generateRandomPoint(tensor::AbstractArray{T,N} where {T,N}) # ::CartesianIndex{N} # where {N}
    #= OLD, specific to 2D
    return CartesianIndex( rand(1:size(tensor, 1)),
                           rand(1:size(tensor, 2))  )
    =#

    #=
    eachIndex is a generator of all indices for an Array.
    See the documentation.
    "For array types that have opted into fast linear indexing
    (like Array), this is simply the range 1:length(A)"
    i.e. a Int64.

    Choose a random element of that set.
    ??? Don't know if it is fast.
    =#
    linearIndex = rand(eachindex(tensor))

    #=
    Convert to CartesianIndex.

    TODO as noted above, it might already be a CartesianIndex.

    Recommended way to convert to a CartesianIndex is
    CartesianIndices(tensor)[linearIndex]
    But requires expensive division.
    But we do that later.
    =#
    result = CartesianIndices(tensor)[linearIndex]

    return result
end

#=
Like the above, but also require the point to be masked.
(The mask true indicates is a desired or selected point.)
=#
function generateRandomMaskedPoint(maskedImage)
    point = 0
    while true
        point = generateRandomPoint(maskedImage.image)
        if maskedImage.mask[point]
            break
        end
    end
    return point
end

#=
Set elements of masked region to a null value

Re the algorithm, calling this is not necessary,
but aids in debugging (so you can see in the result what might be wrong.)

Typically typeof(value) is  RGB
=#
function nullifyMaskedRegion(maskedImage)

    nullValue = nullElementForTensor(maskedImage.image)

    # broadcasting using a logical mask as index
    maskedImage.image[maskedImage.mask] .= nullValue
end

#=
The mask defines a selected, non-rectangular region.
=#
function isPointSelectedInMaskedImage(
        maskedImage::MaskedImage{ValueType, DimensionCount},
        point::CartesianIndex{DimensionCount}
        )::Bool where {ValueType, DimensionCount}
    # is the mask value true at the point ?
    return maskedImage.mask[point]
end

#=
Is the point in bounds?
=#
#= OLD, framed
function isPointInBoundsOfImage(pointInImage)
    # use Julia method
    return checkbounds(Bool, pointInImage.maskedImage.image, pointInImage.point)
end
=#
function isPointInBoundsOfImage(
        maskedImage::MaskedImage{ValueType, DimensionCount},
        point::CartesianIndex{DimensionCount}
        )::Bool where {ValueType, DimensionCount}
    # use Julia method
    return checkbounds(Bool, maskedImage.image, point)
end

#=
FAIL: too slow and allocates memory for the excepttion?

Implementation uses an exception builtin to Julia.
Exception is expected, *not* a program error.

isInBounds = true
try
    # access the point, result is unused
    pointInImage.maskedImage.image[pointInImage.point]
catch BoundsError
    # @debug "not isPointInBoundsOfImage"
    isInBounds = false
end
return isInBounds
=#

#=
Is the given point in bounds of its MaskedImage
AND selected (mask is true)

In original C code clippedOrMaskedCorpus()
=#
function isInBoundsAndSelected(maskedImage, wildPoint)::Bool
    return (   isPointInBoundsOfImage(maskedImage, wildPoint)
            && isPointSelectedInMaskedImage(maskedImage, wildPoint)
            )
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# Random Dictionary based example
using PyPlot

dim     = 64
n_elems = 1024;
D       = randn(dim,n_elems);
s       = randn(dim)
sparse  = lca(s , D, T_soft=0.4)

println("Original  Sparsity = ", sum(s.!=0))
println("Original  RMS      = ", sqrt(mean(s.^2)))
println("Resulting Sparsity = ", sum(sparse.!=0))
println("Recovery  RMS      = ", sqrt(mean((s.-D*sparse).^2)))

figure(1)
PyPlot.clf()
title("Sparse Realization of input")
xlabel("Dictionary Element")
ylabel("Magnitude")
stem(sparse)

figure(2)
PyPlot.clf()
title("Reconstruction of random input")
plot(s)
plot(D*sparse)
legend(["Original", "Recovered"])