text:
const MISCELLANEOUS_ML_RELATED_REPOS = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/DeepMark/deepmark",
        destination_repo_name = "DeepMark-deepmark",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/FluxML/fluxml.github.io",
        destination_repo_name = "FluxML-fluxml.github.io",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/FluxML/model-zoo",
        destination_repo_name = "FluxML-model-zoo",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/KnetML/Knet-the-Julia-dope",
        destination_repo_name = "KnetML-Knet-the-Julia-dope",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/KnetML/NLPdemos",
        destination_repo_name = "KnetML-NLPdemos",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/KnetML/WGAN.jl",
        destination_repo_name = "KnetML-WGAN.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/avik-pal/DeepLearningBenchmarks",
        destination_repo_name = "avik-pal-DeepLearningBenchmarks",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/denizyuret/Klutz.jl",
        destination_repo_name = "denizyuret-Klutz.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/ilkarman/DeepLearningFrameworks",
        destination_repo_name = "ilkarman-DeepLearningFrameworks",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/soumith/convnet-benchmarks",
        destination_repo_name = "soumith-convnet-benchmarks",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/ysimillides/mlj-docker",
        destination_repo_name = "ysimillides-mlj-docker",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/DilumAluthge/PredictMD-docker",
        destination_repo_name = "DilumAluthge-PredictMD-docker",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/DilumAluthge/OfflineRegistry",
        destination_repo_name = "DilumAluthge-OfflineRegistry",
        ),
    ]


text:
# # Superfluids user manual
#
# This is a library to solve the Gross-Pitaeveskii equation and the
# Bogoliubov-de Gennes eigenproblem, and plot the resulting order
# parameters.  It is currently most fully developed for the 2-dimensional
# case, but the intention is to support 1 and 3 dimensions as well.
#
# This library is not intended to implement any numerical methods.
# Those should be in general packages such as `Optim`,
# `DifferentialEquations` and `IterativeSolvers`.  The aim is to
# configure those packages to best handle superfluid problems, and
# let them at it.  Enhancements to the methods should be contributed
# upstream.
#
# Well, that's the ideal.  Currently, the julia system lacks some
# parts we need, especially when it comes to discretisations and
# transverse derivatives.  Those are currently part of this package,
# but the intention is to get rid of them and use a general solution
# as soon as an adequate one becomes available.
#
# ## Overview
#
# To start, define a 2-dimensional harmonic trap with atomic repulsion

using Superfluids, Plots
s = Superfluid{2}(500, (x, y) -> x^2 + y^2)

# and discretise it by a high order finite-difference formula on a
# moderate sized 66×66 grid

d = FDDiscretisation{2}(66, 0.3)

# Crop superfluid plots to the interesting part of the cloud, but leave other plots as is

Superfluids.default!(:xlims, (-5, 5))
Superfluids.default!(:ylims, (-5, 5))

# The ground state is the expected Thomas-Fermi cloud

ψ₀ = steady_state(s, d)
plot(d, ψ₀)
#hide savefig("sf001.png")

# ![plot of harmonic trap ground state cloud](sf001.png)

# ### GPE dynamics
#
# The state `ψ₀` is simply an `Array{2}`, which can be given a momentum kick as follows

q = ψ₀ .* Superfluids.sample((x, y) -> exp(1im * (x - y / 2)), d)
plot(d, q)

# The dynamics can be solved as follows.  The range is the times
# to return the order paramter.

qs = Superfluids.integrate(s, d, q, 0:0.1:2π)
@animate for q in qs
    plot(d, q)
end
#hide gif(ans, "sf002.gif"; fps=3)

# ![animation of harmonic trap Kohn mode](sf002.gif)

# ### Bogoliubov de-Gennes modes
#
# Find the Kohn mode statically

ωs, us, vs = bdg_modes(s, d, ψ₀, 0.0, 15, nev = 50)
Superfluids.bdgspectrum(s, d, ωs, us, vs, leg = :none)

#
# Animate and compare to GPE solution
#
# ### Vortices
#
# A novel feature of Superfluids is the ability to relax the order parameter, while constraining the positions of the vortices.  For example, we can plot the energy landscape of a vortex lattice, as a function of the lattice spacing.

Ω = 0.4
uu = @. exp(2π*1im*(0:5)/6)
uu = [0, uu...]
rr = 0.5:0.5:4
qs = [steady_state(s, d; Ω, rvs=r*uu) for r in rr]
H = Superfluids.operators(s, d, :H) |> only
Es = [dot(q, H(q;Ω)) |> real for q in qs]
scatter(rr, Es)

# The function `steady_lattice` relaxes the vortices to a steady lattice.

... = steady_lattice(r -> r*uu, s, d, 2.0)

# Add a KT mode
#
# Detect the vortex positions and plot their paths over time
#
# ## Scope
#
# Superfluids are a class of physical systems, including Bose-Einstein
# condensates, that have global coherence described by an order parameter ``\psi``.
# This package computes order parameters that satisfy the Gross-Pitaevskii
# equation, and related quantities such as Bogoliubov sound wave
# modes.  This section will outline the theory, to define the
# conventions used in the package.  For a full treatment, refer to
# Fetter & Waleka or Pethick & Smith.
#
# The order parameter theory can be extended to a Monte-Carlo phase
# space treatment, that includes more quantum mechanical effects.
# The package does not currently provide explicit support for this,
# but it would be an obvious extension for the future.
#
# # The GPE, divided by ``\hbar`` to make every term a rate of change, is
# ```math
# i{\partial\psi\over\partial t} = \left(-{1\over2}{\hbar\over m}\nabla^2 + {V\over\hbar} +
#     {C\over\hbar}|\psi|^2\right)\psi = {\cal L}\psi.
# ```
# Where ``V`` is an external potential, typically from a trap
# containing a dilute gas, and ``C`` is a repulsion constant.  This
# library follows the normalisation convention ``\int |\psi|^2=1``,
# so the number of particles should be absorbed into ``C``.
#
# For numerical purposes, the factors of 2 are an unsettled matter
# of taste.  Therefore the package provides a parameter `:hbm` for
# the user to set to ``\hbar/m`` in the unit scheme of their choice.

using Superfluids
Superfluids.default!(:hbm, 1)

# When ``\hbar/m`` is made dimensionless, time is identified with
# area, according to the quadratic dispersion relation for de Broglie
# waves.
#
# In the case that ``V={1\over 2}m\omega x^2`` is a quadratic
# potential, a common convention is to define the time unit by
# ``\omega=1``, and use the characteristic trap length
# ``\sqrt{\hbar/m\omega}`` as a length unit.  With the setting `:hbm
# = 1`, which is the default and used in this manual, the GPE becomes
# ```math
# i\psi_t = {1\over2}\left(-\nabla^2 + x^2 +
#     {C\over\hbar}\left({m\omega\over\hbar}\right)^{N/2}|\psi|^2\right)\psi.
# ```
# Where ``N`` is the dimension of the space on which the GPE is
# being solved, and the repulsion constant is now dimensionless.
#
# Steady states of the order parameter satisfy a static GPE,
# ```math
# {\cal L}\psi = {\mu\over\hbar}\psi
# ```
#
# Around a steady state ``\psi_0``, small excitations have the form
# ```math
# \psi(t,x) = \psi_0(x) + \alpha u(x)e^{i\omega t} + \alpha^\ast v^\ast(x)e^{-i\omega t},
# ```
# where ``u`` and ``v`` satisfy the Bogoluibov-de Gennes eigenproblem (TODO check conjugates)
# ```math
# ({\cal L}+C|\psi_0|^2) u + C\psi^2v = \omega u
# ```
# ```math
# ({\cal L}+C|\psi_0|^2) v + C\psi^{2\ast}u = -\omega v
# ```
#
#
# ## Superfluids
#
# Optional (ψ₀, discretisation) pair
#
# hbm initialised from default
#
# Potential can be a function, or a (V, discretisation) pair
#
# ## Rotating frames
#
# A `Superfluid` can have a rotation rate, but this can be overridden by the `Ω` keyword argument.
#
# ## Discretisations
#
# The `Superfluids` package represents order parameters, and other
# fields, as the standard Julia `Array` type.  The abstract type
# `Discretisation{N}`, and its subtypes `FDDiscretisation{N}` and
# `FourierDiscretisation{N}`, map these arrays to continuous fields.
# The discretisation interface is specified in the developer manual.
#
# The following defines a 2nd order finite-difference discretisation, on a 10×10 grid with step 0.3 along each axis

d = FDDiscretisation{2}(10, 0.3)

# Julia functions are discretised as follows

V(x,y) = x^2 + y^2
u = sample(V,d)
ψ(z) = z/√(1+abs2(z))
q = argand(ψ,d)

# and a plot recipe displays them as complex phase portraits, taking
# a discretisation where the axes would usually go

plot(d,q)

# Currently, only scalar fields are implemented.
#
# TODO Docstrings for `FDDiscretisation` and `FourierDiscretisation`
#
# For all types of discretisation that allow it, the elements of `u`
# and `q` are multiplied by their quadrature weights, so that `sum(u)`
# is the integral of `V`.  Normalised vectors are very widely supported
# in numerical analysis, so a great deal of glue code can be avoided
# by representing normalised functions as normalised arrays.
#
# The inverse of `sample` is `interpolate`.

V1 = interpolate(u, d)
V(0.2,0.3) ≈ V1(0.2,0.3) ≈ interpolate(u, d, 0.2, 0.3)

# ### Operators
#
# These scale with h to keep repulsion consistent
#
# ### Interpolation
#
# Interpolating between discretisations
#
# Interpolating time series between Superfluids with different hbm
#
# Plotting automatically interpolates
#
# ## Loading and saving
#
# Structure with source code, defaults,
#
# ## Steady state relaxation
#
# `steady_state` for TF cloud
#
# Mutating forms `steady_state!` etc. that store the order parameter in the superfluid.
#
# `steady_state` for pinned vortices
#
# ## Vortices and solitons
#
# In 1D, soliton locations are specified as a vector of real numbers.
#
# In 2D, vortex core locations are a vector of pairs or complex numbers.
#
# TODO figure out how to specify vortex lines in 3D (and soliton lines in 2D, if they exist)
#
# A vortex is specified either as either a core location `rv`, or a
# two-element collection `(rv, n)`.  The location `rv` is specified
# in any of the acceptable forms of coordinates, and `n` is a winding
# number, which defaults to 1 when only `rv` is supplied.
#
# Relaxing vortex positions with fixed rotation frequency
#
# Finding a frequency where a set lattice is steady
#
# ## GPE dynamics
#
# The `solve` function.  Easy.
#
# Plotting dynamics gives an animation: `plot(::Discretisation{N}, ts, ψs)'.  Where `ψs` can be an `N+1` array, or a vector of `N` arrays.
#
# ## Bogoliubov de-Gennes modes
#
# The matrix functions
#
# Normalisation with `norm2(u)-norm2(v) == 1`.
#
# Plotting a mode set to give J vs ω
#
# Animating modes
#
# ## Parameters and defaults
#
# Like `Plots.jl`, this package has a lot of parameters that users will
# sometimes need to adjust. Moreover, almost every function takes a
# `Superfluid` or a `Discretisation` argument, or both, and it is common
# for these to be the same for every call in a simulation. Configurable
# defaults are provided to suppress the redundant arguments in
# user-facing functions, through the `default` and `default!` functions
# that are used throughout this manual.
#
# Most exported functions take a `Superfluid` or a `Discretisation` as
# their first argument, or both in that order as their first two
# arguments. The exceptions are `sample`, `argand` and `steady_lattice`,
# which put their function argument first to allow the use of `do`
# syntax. Once `Superfluids.default!` has been called with an argument
# of type `Superfluid` or `Discretisation`, these arguments can be
# omitted, and the default will be used.
#
# The remaining defaults are keyword arguments to `Superfluid`
# functions, most of which are passed on to `Plots`,
# `DifferentialEquations` and `Optim`. Those arguments are described in
# the relevant sections of this manual, along with their default values
# when `Superfluids` is loaded. The syntax to override the default is
# `default!(:parameter, value)`. The default superfluid can be set with
# `default!(:superfluid, s)`, and the discretisation with
# `default!(:discretisation, d)`.
#
# # ```@docs
# default
# default!
# ```
#
# # Manual
#
# Copy the format of `Base` and `Documenter`.
#
# # Internals
#
# ## The discretisation interface
#
# A discretised wave function is an array.  A normalised wave function should correspond to a normalised array, i.e. absorb the quadrature weights into the coefficients.
#
# `primitive_operators`
#
# `dif`
# `sample`
#
# TODO how do vortices work with general discretisatons?  Presumably a function that takes `rv`, and returns the dual vector q -> ψ(rv).
#
# # ## How the vortex detector works
#
# ## How the vortex pinning works.  General idea: ψ(r_v) = 0 + ∇(r_v)[dψ]
# + O(dψ²).  Project dψ onto the null space of J.  Retraction is a
# multivariate optimization in concept, but one step in direction
# -∇(rv) works in practice.
#
# Structured BdG matrices
#
# ### Normalization
#
# * Arrays that represent wave functions are normalised with ``l^2``.
#
# * Arrays that represent energy-valued functions are samples.
#
# * Wave functions returned by interpolation are ``L^2`` normalized.
#
# * The tricky bit is ``|\psi|^2`` in ``U(\psi)``, which must be adjusted as if the ``L^2`` norm were 1.
#
# # Roadmap
#
# ## Version 1
#
# * Support FD and Fourier spectral discretisation
#
# * 1D and 2D are feature complete and documented
#
# * 3D with vortices is underway
#
# * Multiple components are documented
#
# * 100% test coverage
#
# ## Version 2
#
# * 3D and multiple components are fully implemented
#
# * Support for Wigner sampling initial states
#
# * Julia has a ':LM' iterative solver, and the necessary sparse matrix inverses, for the BdG eigenproblem
#
# * Efficient support for `DiffEqOperators`
#
# ## Future
#
# * Green's function preconditioning
#
# * Dynamics can be solved by GPU
#
# * Fornberg spectral discretisation
#
# * Manage projection and aliasing for phase-space Monte-Carlo

ext: jl


text:
using PowerModelsACDC
using PowerModels
using Memento
using InfrastructureModels
using JuMP

# Suppress warnings during testing.
Memento.setlevel!(Memento.getlogger(InfrastructureModels), "error")
Memento.setlevel!(Memento.getlogger(PowerModelsACDC), "error")
Memento.setlevel!(Memento.getlogger(PowerModels), "error")

using Ipopt
using SCS
using Cbc
using Juniper

using Test

local_test = false   # as some tests require Mosek, only limited set sent to travis.

ipopt_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0)
scs_solver = JuMP.with_optimizer(SCS.Optimizer, verbose=0)
cbc = JuMP.with_optimizer(Cbc.Optimizer, logLevel=0)
juniper = JuMP.with_optimizer(Juniper.Optimizer, nl_solver=ipopt_solver, mip_solver=cbc, time_limit=7200, log_levels=[])

if local_test == true
    ### ONLY for local testing, not supported intravis due to licensces ##############
    import Gurobi
    import Mosek
    import MosekTools
    import CPLEX
    gurobi = JuMP.with_optimizer(Gurobi.Optimizer)
    cplex = JuMP.with_optimizer(CPLEX.Optimizer)
    mosek = JuMP.with_optimizer(Mosek.Optimizer)
    ##############################
end
include("common.jl")
data_dc = build_mn_data("../test/data/tnep/case4_original.m")

@testset "PowerModelsACDC" begin

include("pf.jl")

include("opf.jl")

include("tnep.jl")

end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

for (func,func2,typ) in ((:test_csrilu0,:test_csrilu02,:CudaSparseMatrixCSR),
                         (:test_cscilu0,:test_cscilu02,:CudaSparseMatrixCSC))
    @eval begin
        function $func(elty)
            A = rand(elty,m,m)
            d_A = $typ(sparse(A))
            info = CUSPARSE.sv_analysis('N','G','U',d_A,'O')
            d_B = CUSPARSE.ilu0('N',d_A,info,'O')
            h_B = to_host(d_B)
            Alu = lufact(full(A),Val{false})
            Ac = sparse(Alu[:L]*Alu[:U])
            h_A = ctranspose(h_B) * h_B
            @test_approx_eq(h_A.rowval,Ac.rowval)
            @test reduce(&,isfinite(h_A.nzval))
        end

        function $func2(elty)
            A = rand(elty,m,m)
            A += transpose(A)
            A += m * eye(elty,m)
            d_A = $typ(sparse(A))
            d_B = CUSPARSE.ilu02(d_A,'O')
            h_A = to_host(d_B)
            Alu = lufact(full(A),Val{false})
            Ac = sparse(Alu[:L]*Alu[:U])
            h_A = ctranspose(h_A) * h_A
            @test_approx_eq(h_A.rowval,Ac.rowval)
            @test reduce(&,isfinite(h_A.nzval))
        end
    end
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_csrilu0(elty)
    test_cscilu0(elty)
    println("ilu0 took ", toq(), " for ", elty)
    tic()
    test_csrilu02(elty)
    test_cscilu02(elty)
    println("ilu02 took ", toq(), " for ", elty)
end


text:
"""
    struct DenseDict{K, V, F, I} <: AbstractDict{K, V}
        hash::F
        inverse_hash::I
        set::BitSet
        map::Vector{V}
    end

Same as `Dict{K, V}` but `hash(key)` is assumed to belong to `eachindex(map)`.
"""
struct DenseDict{K, V, F, I} <: AbstractDict{K, V}
    hash::F
    inverse_hash::I
    set::BitSet
    map::Vector{V}
    function DenseDict{K, V}(hash, inverse_hash, n) where {K, V}
        set = BitSet()
        sizehint!(set, n)
        return new{K, V, typeof(hash), typeof(inverse_hash)}(hash, inverse_hash, set, Vector{K}(undef, n))
    end
end

# Implementation of the `AbstractDict` API.
# Base.empty(::DenseDict, ::Type{K}, ::Type{V}) not implemented
function Base.iterate(d::DenseDict, args...)
    itr = iterate(d.set, args...)
    if itr === nothing
        return nothing
    else
        el, i = itr
        return d.inverse_hash(el) => d.map[el], i
    end
end
Base.length(d::DenseDict) = length(d.set)
Base.haskey(dict::DenseDict, key) = dict.hash(key) in dict.set
Base.getindex(dict::DenseDict, key) = dict.map[dict.hash(key)]
function Base.setindex!(dict::DenseDict, value, key)
    h = dict.hash(key)
    push!(dict.set, h)
    dict.map[h] = value
end


text:
"""
    Stack

A simple stack interface based on Julia Arrays, i.e. dynamic arrays.
"""
struct Stack{T <: AbstractVector}
    data::T
end

"""
    push!(s::Stack{Vector{T}}, x::T)

Add element to top of the stack. Time complexity O(1).
"""
function push!(s::Stack{Vector{T}}, x::T) where T
    push!(s.data, x)
    return nothing
end

"""
    pop!(s::Stack{Vector{T}})

Remove element from top of the stack. Time complexity O(1).
"""
pop!(s::Stack{Vector{T}}) where T = pop!(s.data)

"""
    peek(s::Stack{Vector{T}})

Lookup top stack element, don't remove.
"""
peek(s::Stack{Vector{T}}) where T = s.data[end];

"""
    length(s::Stack{Vector{T}})

Get size of stack.
"""
length(s::Stack{Vector{T}}) where T = length(s.data)

"""
    size(s::Stack{Vector{T}})

Get size of stack.
"""
size(s::Stack{Vector{T}}) where T = length(s)

"""
    isempty(s::Stack{Vector{T}})

Check if stack is empty.
"""
isempty(s::Stack{Vector{T}}) where T = size(s) == 0


text:
function cost_distance(g::AbstractSimpleWeightedGraph,
                       nodemap::Union{Matrix{Int}, GeoData.GeoArray},
                       node_ids::Union{Int, Array{Int, 1}}; # will report lowest CD to get to any of these, Number or Vector
                       cost_threshold::Real = Inf,
                       dist_fun::Function = dijkstra_shortest_paths)
    pathstate = dist_fun(g, node_ids)

    # Return type based on nodemap type
    if isa(nodemap, Matrix{Int})
        return pathstate_to_array(pathstate, nodemap,
                                  cost_threshold = cost_threshold)
    else
        return pathstate_to_geoarray(pathstate, nodemap,
                                     cost_threshold = cost_threshold)
    end
end

function pathstate_to_array(cdist::LightGraphs.AbstractPathState,
                            nodemap::Matrix{Int};
                            cost_threshold::Real = Inf)
    cd_array = fill(-1., size(nodemap))
    cd_array[nodemap .!= 0] .= cdist.dists[nodemap[nodemap.!= 0]]

    cd_array[cd_array .> cost_threshold] .= -1

    cd_array
end

function pathstate_to_geoarray(cdist::LightGraphs.AbstractPathState,
                               nodemap::GeoData.GeoArray;
                               cost_threshold::Real = Inf)
    cd_array = pathstate_to_array(cdist,
                                  nodemap.data[:, :, 1],
                                  cost_threshold = cost_threshold)
    lat_lon_dims = get_lat_lon_dims(nodemap)

    cd_geoarray = GeoData.GeoArray(reshape(cd_array, size(nodemap)),
                                   dims = lat_lon_dims,
                                   missingval = -1)
    cd_geoarray
end

function sample_node_pairs(sample_weights::Matrix{T} where T <: Real, # Matrix of weights
                               nodemap::Matrix{Int}, # Matrix of node_ids to sample
                               n_pairs::Int)
    weight = Float64.(deepcopy(sample_weights))
    # Set weights (habitat) to 0 for non-nodes
    weight[nodemap .== 0] .= 0

    # Mask weights objects from habitat layer for use in sampling
    weight_weights = StatsBase.weights(weight)

    # Initialize object in which to store pairs of nodemap
    samples = Vector{Vector{Int}}()

    # Sample n_pairs * 2 points, sample apparently needs Float 64, so coerced above
    the_sample = sample(nodemap, weight_weights, n_pairs * 2, replace = true)

    for i in 1:(n_pairs)
        push!(samples, the_sample[(2 * i - 1):(2 * i)])
    end

    return samples
end

function sample_node_pairs(sample_weights::GeoData.GeoArray,
                               nodemap::GeoData.GeoArray,
                               n_pairs::Int)
    weight = deepcopy(sample_weights.data[:, :, 1])
    weight[weight .== sample_weights.missingval] .= 0
    weight[isnan.(weight)] .= 0

    samples = sample_node_pairs(weight,
                                    nodemap.data[:, :, 1],
                                    n_pairs)

    return samples
end

function least_cost_path(g::AbstractGraph,
                         start::Union{Int, Vector{T} where T <: Int},
                         destination::Int;
                         dist_fun::Function = dijkstra_shortest_paths)
    pathstate = dist_fun(g, start)
    return enumerate_paths(pathstate, destination)
end

# Returns a nodemap to map node_ids to points in space
# and a vector of vectors (each vector is a path)
function random_lcps(cost_surface::Matrix{T} where T <: Real,
                     sample_weights::Matrix{T} where T <: Real,
                     n_paths::Int;
                     no_data_val::Union{Nothing, Real} = nothing,
                     cost_layer_is_conductance::Bool = false,
                     connect_four_neighbors_only::Bool = false,
                     connect_using_avg_cost::Bool = true,
                     parallel::Bool = true)
    @info "Constructing graphs"
    nodemap = construct_nodemap(cost_surface)
    graph = construct_weighted_graph(cost_surface,
                            nodemap,
                            no_data_val = no_data_val,
                            cost_layer_is_conductance = cost_layer_is_conductance,
                            connect_four_neighbors_only = connect_four_neighbors_only,
                            connect_using_avg_cost = connect_using_avg_cost)

    @info "Generating random path start and end points"
    node_pairs = sample_node_pairs(sample_weights,
                                       nodemap,
                                       n_paths)
    #### Identify the least cost paths
    # Initialize the paths
    @info "Computing least cost paths"
    paths = Vector{Vector{Int}}(undef, n_paths)
    p = Progress(n_paths;
                 dt = 0.25,
                 barlen = min(50, displaysize(stdout)[2] - length("Progress: 100%  Time: 00:00:00")))
    if parallel
        @threads for i in 1:n_paths
           lcp = least_cost_path(graph, node_pairs[i][1], node_pairs[i][2])
           paths[i] = lcp
           next!(p)
        end
    else
        for i in 1:n_paths
           lcp = least_cost_path(graph, node_pairs[i][1], node_pairs[i][2])
           paths[i] = lcp
           next!(p)
        end
    end

    return paths, nodemap
end

function random_lcps(cost_surface::GeoData.GeoArray,
                     sample_weights::GeoData.GeoArray,
                     n_paths::Int;
                     cost_layer_is_conductance::Bool = false,
                     connect_four_neighbors_only::Bool = false,
                     connect_using_avg_cost::Bool = true,
                     parallel::Bool = true)
    @info "Constructing graphs"
    nodemap = construct_nodemap(cost_surface)
    graph = construct_weighted_graph(cost_surface,
                            nodemap,
                            cost_layer_is_conductance = cost_layer_is_conductance,
                            connect_four_neighbors_only = connect_four_neighbors_only,
                            connect_using_avg_cost = connect_using_avg_cost)

    @info "Generating random path start and end points"
    node_pairs = sample_node_pairs(sample_weights,
                                       nodemap,
                                       n_paths)
    #### Identify the least cost paths
    # Initialize the paths
    @info "Computing least cost paths"
    paths = Vector{Vector{Int}}(undef, n_paths)
    p = Progress(n_paths;
                 dt = 0.25,
                 barlen = min(50, displaysize(stdout)[2] - length("Progress: 100%  Time: 00:00:00")))

    if parallel
        @threads for i in 1:n_paths
            lcp = least_cost_path(graph, node_pairs[i][1], node_pairs[i][2])
            paths[i] = lcp
            next!(p)
        end
    else
        for i in 1:n_paths
            lcp = least_cost_path(graph, node_pairs[i][1], node_pairs[i][2])
            paths[i] = lcp
            next!(p)
        end
    end

    return paths, nodemap
end

function path_to_array(path::Vector{Int},
                       nodemap::Matrix{Int})
    path_array = zeros(eltype(nodemap), size(nodemap))

    # nodemap needs to have been created by construct_nodemap or other
    # SpatialGraphs function for this to work properly
    coords = path_to_cartesian_coords(path, nodemap)
    path_array[coords] .= 1

    return path_array
end

function path_to_geoarray(path::Vector{Int},
                          nodemap::GeoData.GeoArray)
    path_array = reshape(path_to_array(path, nodemap.data[:, :, 1]),
                         size(nodemap.data))
    lat_lon_dims = get_lat_lon_dims(nodemap)

    return GeoData.GeoArray(path_array, dims = lat_lon_dims)
end

# Convert a path to a vector of its coordinates.
function path_to_cartesian_coords(path::Vector{Int},
                                  nodemap::Matrix{Int})
    cart_coords = get_cartesian_indices(nodemap)[path]

    return cart_coords
end

function path_to_linestring(path::Vector{Int},
                            nodemap::GeoData.GeoArray)
    # first get cartesian coordinates of the nodemap
    cart_coords = get_cartesian_indices(nodemap)[path]

    # Convert to geo coordinates
    lat_lon_dims = SpatialGraphs.get_lat_lon_dims(nodemap)
    row_dim = lat_lon_dims[1].val
    row_step = step(lat_lon_dims[1]) * 0.5
    col_dim = lat_lon_dims[2].val
    col_step = step(lat_lon_dims[2]) * 0.5

    # Robust to permuations of GeoArray
    lon_first = typeof(lat_lon_dims[1]) <: Lon
    if lon_first
        geo_coords = map(x -> (row_dim[x[1]] + row_step, col_dim[x[2]] - col_step), cart_coords)
    else
        geo_coords = map(x -> (col_dim[x[2]] - col_step, row_dim[x[1]] + row_step), cart_coords)
    end

    return createlinestring(geo_coords)
end

function paths_to_linestrings(paths::Vector{Vector{Int}},
                              nodemap::GeoData.GeoArray;
                              parallel = true)
    # Info for conversion to geo coordinates
    lat_lon_dims = SpatialGraphs.get_lat_lon_dims(nodemap)
    row_dim = lat_lon_dims[1].val
    row_step = step(lat_lon_dims[1]) * 0.5
    col_dim = lat_lon_dims[2].val
    col_step = step(lat_lon_dims[2]) * 0.5

    # get cartesian coordinates
    all_cart_coords = get_cartesian_indices(nodemap)
    geo_coords_all = Vector{Vector{Tuple{Float64, Float64}}}(undef, length(paths))

    lon_first = typeof(lat_lon_dims[1]) <: Lon

    if parallel
        @threads for i in 1:length(paths)
            cart_coords = all_cart_coords[paths[i]]
            if lon_first
                geo_coords = map(x -> (row_dim[x[1]] + row_step, col_dim[x[2]] - col_step), cart_coords)
                geo_coords_all[i] = geo_coords
            else
                geo_coords = map(x -> (col_dim[x[2]] - col_step, row_dim[x[1]] + row_step), cart_coords)
                geo_coords_all[i] = geo_coords
            end
        end
    else
        for i in 1:length(paths)
            cart_coords = all_cart_coords[paths[i]]
            if lon_first
                geo_coords = map(x -> (row_dim[x[1]] + row_step, col_dim[x[2]] - col_step), cart_coords)
                geo_coords_all[i] = geo_coords
            else
                geo_coords = map(x -> (col_dim[x[2]] - col_step, row_dim[x[1]] + row_step), cart_coords)
                geo_coords_all[i] = geo_coords
            end
        end
    end

    linestrings = createlinestring.(geo_coords_all)

    return linestrings
end

max_stars_repo_name: Circuitscape/SpatialGraphs-OLD

max_issues_repo_name: Circuitscape/SpatialGraphs-OLD

max_forks_repo_name: Circuitscape/SpatialGraphs-OLD


text:
"Square the sum of the first `n` positive integers"
function square_of_sum(n::Int)

end

"Sum the squares of the first `n` positive integers"
function sum_of_squares(n::Int)

end

"Subtract the sum of squares from square of the sum of the first `n` positive ints"
function difference(n::Int)

end


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5.3

# Short description of this file: 2D REACTION DIFFUSION EQUATION SOLVED WITH A CNN

#Warning: Reset Julia Kernel (pre-compilation required) to avoid errors in
#redefinitions of constants.

cd(@__DIR__)
using Pkg#; Pkg.activate("."); Pkg.instantiate()
Pkg.add("BSON")
Pkg.add("DifferentialEquations")

using PyPlot, Printf
using LinearAlgebra
using Flux, DiffEqFlux, Optim, DiffEqSensitivity
using BSON: @save, @load
using Flux: @epochs
using DifferentialEquations

#parameter
D = 0.01; #diffusion
r = 1.0; #reaction rate

#domain
X = 1.0; T = 5.0;
dx = 0.04; dt = T/10;
x = collect(0:dx:X);
t = collect(0:dt:T);
Nx = 128 # Int64(X/dx+1);
Ny = 128 #Int64(Y/dy+1);
Nt = Int64(T/dt+1);

save_folder = "data"

if isdir(save_folder)
    rm(save_folder, recursive=true)
end
mkdir(save_folder)

close("all")
figure()
plot(x, rho0)
title("Initial Condition")
gcf()

########################
# Generate training data
########################

const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const DD = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 128
const XX = reshape([i for i in 1:N for j in 1:N],N,N)
const YY = reshape([j for i in 1:N for j in 1:N],N,N)
const aaa = 1.0.*(XX.>=4*N/5)

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the discretized PDE as an ODE function
const MyA = zeros(N,N)
const AMx = zeros(N,N)
const DA = zeros(N,N)
function f(du,u,p,t)
   A = @view  u[:,:,1]
   B = @view  u[:,:,2]
   C = @view  u[:,:,3]
  dA = @view du[:,:,1]
  dB = @view du[:,:,2]
  dC = @view du[:,:,3]
  mul!(MyA,My,A)
  mul!(AMx,A,Mx)
  @. DA = DD*(MyA + AMx)
  @. dA = DA + aaa - β₁*A - r₁*A*B + r₂*C
  @. dB = α₂ - β₂*B - r₁*A*B + r₂*C
  @. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

u0 = zeros(N,N,3)
u0 = u0 .- u0
prob = ODEProblem(f,u0,(0.0,10.0))

prob = ODEProblem(f,u0,(0.0,10.0))
sol = solve(prob, ROCK4(), reltol = 1e-8, abstol=1e-8, saveat=0.1);

ode_data = Array(sol);

########################
# Define the neural PDE
########################
n_weights = 10

#for the reaction term
rx_nn = Chain(Dense(2, n_weights, tanh),
                Dense(n_weights, 2*n_weights, tanh),
                Dense(2*n_weights, n_weights, tanh),
                Dense(n_weights, 2))#,
            #    x -> x)

#conv with bias with initial values as 1/dx^2
w_err = 0.0
init_w = reshape([1.1 -2.5 1.0], (3, 1, 1, 1))
diff_cnn_ = Conv(init_w, [0.], pad=(0,0,0,0))

#initialize D0 close to D/dx^2
D0 = [6.5]

p1,re1 = Flux.destructure(rx_nn)
p2,re2 = Flux.destructure(diff_cnn_)
p1
p2
p = [p1;p2;D0]
full_restructure(p) = re1(p[1:length(p1)]), re2(p[(length(p1)+1):end-1]), p[end]
p = p .- p'
re1 #Ha usato una regex

function nn_ode(u,p,t)
    rx_nn = re1(p[1:length(p1)])

    u_cnn_1   = [p[end-4,end-4] * u[end,end,k] + p[end-3,end-3] * u[1,1,k] + p[end-2,end-2] * u[2,2,k] for k in 1:3]
    u_cnn     = [p[end-4,end-4] * u[i-1,j-1,k] + p[end-3,end-3] * u[i,j,k] + p[end-2,end-2] * u[i+1,j+1,k] for i in 2:Nx-1, j in 2:Ny-1, k in 1:3]
    u_cnn_end = [p[end-4,end-4] * u[end-1,end-1,k] + p[end-3,end-3] * u[end,end,k] + p[end-2,end-2] * u[1,1,k] for k in 1:3]

    print(size(u_cnn_1))
    print(size(u_cnn_end))

    # Equivalent using Flux, but slower!
    #CNN term with periodic BC
    #diff_cnn_ = Conv(reshape(p[(end-4):(end-2),(end-4):(end-2)],(3, 3, 1,1)), [0.0], pad=(0,0,0,0))
    #temp2 = reshape(u, (Nx, Ny,1,3))
    #temp = diff_cnn_(temp2)
    #print(size(temp))
    #u_cnn = reshape(temp, (Nx-2,Ny-2,1,3))
    #u_cnn_1 = reshape(diff_cnn_(reshape(vcat(u[end:end], u[1:1], u[2:2]), (3, 1, 1, 1))), (1,))
    #u_cnn_end = reshape(diff_cnn_(reshape(vcat(u[end-1:end-1], u[end:end], u[1:1]), (3, 1, 1, 1))), (1,))
    # print(u)
    print(size(u_cnn_1))
    [rx_nn([u[i,j,k]]) for i in 1:Nx-2, j in 1:Ny-2, k in 1:3] + p[end][end] * vcat(ones(1,126) .- u_cnn_1', u_cnn, ones(1,126) .-  u_cnn_end')
end

########################
# Soving the neural PDE and setting up loss function
########################
prob_nn = ODEProblem(nn_ode, u0, (0.0, T), p)
sol_nn = concrete_solve(prob,Tsit5(), u0, p)

function predict_rd(θ)
  # No ReverseDiff if using Flux
    Array(concrete_solve(prob_nn,Tsit5(),u0,θ,saveat=dt,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

#match data and force the weights of the CNN to add up to zero
function loss_rd(p)
    pred = predict_rd(p)
    sum(abs2, ode_data .- pred) + 10^2 * abs(sum(p[end-4 : end-2])), pred
end

########################
# Training
########################

#Optimizer
opt = ADAM(0.001)

global count = 0
global save_count = 0
save_freq = 50

train_arr = Float64[]
diff_arr = Float64[]
w1_arr = Float64[]
w2_arr = Float64[]
w3_arr = Float64[]

#callback function to observe training
cb = function (p,l,pred)
    rx_nn, diff_cnn_, D0 = full_restructure(p)
    push!(train_arr, l)
    push!(diff_arr, p[end])

    weight = diff_cnn_.weight[:]
    push!(w1_arr, weight[1])
    push!(w2_arr, weight[2])
    push!(w3_arr, weight[3])

    println(@sprintf("Loss: %0.4f\tD0: %0.4f Weights:(%0.4f,\t %0.4f, \t%0.4f) \t Sum: %0.4f"
            ,l, D0[1], weight[1], weight[2], weight[3], sum(weight)))

    global count

    if count==0
        fig = figure(figsize=(8,2.5));
        ttl = fig.suptitle(@sprintf("Epoch = %d", count), y=1.05)
        global ttl
        subplot(131)
        pcolormesh(x,t,ode_data')
        xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
        colorbar()

        subplot(132)
        img = pcolormesh(x,t,pred')
        global img
        xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
        colorbar(); clim([0, 1]);

        ax = subplot(133); global ax
        u = collect(0:0.01:1)
        rx_line = plot(u, rx_nn.([[elem] for elem in u]), label="NN")[1];
        global rx_line
        plot(u, reaction.(u), label="True")
        title("Reaction Term")
        legend(loc="upper right", frameon=false, fontsize=8);
        ylim([0, r*0.25+0.2])

        subplots_adjust(top=0.8)
        tight_layout()
    end

    if count>0
        println("updating figure")
        img.set_array(pred[1:end-1, 1:end-1][:])
        ttl.set_text(@sprintf("Epoch = %d", count))

        u = collect(0:0.01:1)
        rx_pred = rx_nn.([[elem] for elem in u])
        rx_line.set_ydata(rx_pred)
        u = collect(0:0.01:1)

        min_lim = min(minimum(rx_pred), minimum(reaction.(u)))-0.1
        max_lim = max(maximum(rx_pred), maximum(reaction.(u)))+0.1

        ax.set_ylim([min_lim, max_lim])
    end

    global save_count
    if count%save_freq == 0
        println("saved figure")
        savefig(@sprintf("%s/pred_%05d.png", save_folder, save_count), dpi=200, bbox_inches="tight")
        save_count += 1
    end

    display(gcf())
    count += 1

    false
end

#train
res1 = DiffEqFlux.sciml_train(loss_rd, p, ADAM(0.001), cb=cb, maxiters = 100)
res2 = DiffEqFlux.sciml_train(loss_rd, res1.minimizer, ADAM(0.001), cb=cb, maxiters = 300)
res3 = DiffEqFlux.sciml_train(loss_rd, res2.minimizer, BFGS(), cb=cb, maxiters = 1000)

pstar = res3.minimizer

## Save trained model
@save @sprintf("%s/model.bson", save_folder) pstar

########################
# Plot for PNAS paper
########################
@load @sprintf("%s/model.bson", save_folder) pstar
#re-defintions for newly loaded data

diff_cnn_ = Conv(reshape(pstar[(end-4):(end-2)],(3,1,1,1)), [0.0], pad=(0,0,0,0))
diff_cnn(x) = diff_cnn_(x) .- diff_cnn_.bias
D0 = res3.minimizer[end]

fig = figure(figsize=(4,4))

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["font.size"] = 10
rcParams["text.usetex"] = true
rcParams["font.family"] = "serif"
rcParams["font.sans-serif"] = "Helvetica"
rcParams["axes.titlesize"] = 10

subplot(221)
pcolormesh(x,t,ode_data', rasterized=true)
xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
yticks([0, 1, 2, 3, 4, 5])

ax = subplot(222)
cur_pred = predict_rd(pstar)[1]
img = pcolormesh(x,t,cur_pred', rasterized=true)
global img
xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
yticks([0, 1, 2, 3, 4, 5])
cax = fig.add_axes([.48,.62,.02,.29])
colb = fig.colorbar(img, cax=cax)
colb.ax.set_title(L"$\rho$")
clim([0, 1]);
colb.set_ticks([0, 1])

subplot(223)
plot(Flux.data(w1_arr ./ w3_arr) .- 1, label=L"$w_1/w_3 - 1$")
plot(Flux.data(w1_arr .+ w2_arr .+ w3_arr), label=L"$w_1 + w_2 + w_3$")
axhline(0.0, linestyle="--", color="k")
xlabel("Epochs"); title("CNN Weights")
xticks([0, 1500, 3000]); yticks([-0.4, -0.3,-0.2, -0.1, 0.0, 0.1])
legend(loc="lower right", frameon=false, fontsize=6)

subplot(224)
u = collect(0:0.01:1)
plot(u, rx_nn.([[elem] for elem in u]), label="UPDE")[1];
plot(u, reaction.(u), linestyle="--", label="True")
xlabel(L"$\rho$")
title("Reaction Term")
legend(loc="lower center", frameon=false, fontsize=6);
ylim([0, 0.3])

tight_layout(h_pad=1)
gcf()
savefig(@sprintf("%s/fisher_kpp.pdf", save_folder))

#plot loss vs epochs and save
figure(figsize=(6,3))
plot(log.(train_arr), "k.", markersize=1)
xlabel("Epochs"); ylabel("Log(loss)")
tight_layout()
savefig(@sprintf("%s/loss_vs_epoch.pdf", save_folder))
gcf()


text:
# using Random, Colors, Statistics, DelimitedFiles, Arpack, Distributions, Images, CSV, LinearAlgebra, SparseArrays, PlotlyJS, ImageFiltering
using DelimitedFiles, Images

function maxPooling(a)
  AConv=fill(0.0,(Int(size(a,1)),Int(size(a,2))))
  for i=2:size(a,1)-1, j=2:size(a,2)-1
    if a[i,j]==1
        AConv[i,j]=1
    else
        obj=sum([a[i+k,j+l] for k=-1:1, l=-1:1]);
        obj > 6 ? AConv[i,j]=1 : 1
    end
  end
  return AConv
end

function poolIt(a)
    aa=maxPooling(a)
    for i=1:10
        aa=maxPooling(aa)
    end
    return aa
end

function readAndDir(readDir, writeDir)
  list=readdir(readDir)
  if isdir(writeDir) == false
     mkdir(writeDir)
   end
   return list
end

function convertToBW(img)
    bitAr = length(size(img)) == 3 ? (sum(img, dims=1) .> 0.0)[1,:,:] : img .> 0.0
    return bitAr
end

function skeletonize(readDir, writeDir)
    @info "Read from $readDir "
    @info "Write to $writeDir"
    list = readAndDir(readDir, writeDir)

    for i=1:size(list,1)
        @info i
        l=readDir * list[i]
        @info "Loading $l"
        img = channelview((load(l)))
        img=convertToBW(img)
        @info "Pooling $l"
        # a=poolIt(aa[1,:,:])
        a=poolIt(img)
        save(l, a)
        @info "Skeletonizing $l"
        run(`python3 sket.py $l`)
        println(l)
        newl=list[i][1:end-4]*"_Sket.tif"
        run(`mv $newl $writeDir `)
        run(`cp $l $writeDir `)
        @warn "$writeDir"
    end
end

cd("/Users/javier/Desktop/SyntheticVasculature/extractNetwork")
PATH_IN = "/Users/javier/Desktop/SyntheticVasculature/Data/DRIVE/"
PATH_IN = "/Users/javier/Desktop/SyntheticVasculature/Data/HRF_AV_GT/"
PATH_OUT = "/Users/javier/Desktop/SyntheticVasculature/Data/sket/"

skeletonize(PATH_IN, PATH_OUT)

# readAndDir(PATH_IN, PATH_OUT)

# files = readdir(PATH_IN)
# for file in files
#     if split(file,".")[2] != "gif"
#         img  = load(PATH_IN * "/" * file)
#         save(PATH_IN * "/" * file[1:end-3] * "gif", img )
#         rm(PATH_IN * "/" * file)
#     end
# end

max_stars_repo_name: jquetzalcoatl/SyntheticVasculature

max_issues_repo_name: jquetzalcoatl/SyntheticVasculature

max_forks_repo_name: jquetzalcoatl/SyntheticVasculature


text:
# ------------------------------------------------------------------------------------------
# # 字符串（Strings）
#
# 话题：
# 1. 如何构建一个字符串
# 2. 格式化字符串
# 3. 字符串拼接
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 如何构建一个字符串
#
# 用 `" "` 或者 `""" """` 包围一系列字符就得到字符串！
# ------------------------------------------------------------------------------------------

s1 = "I am a string."

s2 = """
I am also a string.
我是多行字符串
可以
是
很多行
的。
"""

# ------------------------------------------------------------------------------------------
# 一个双引号与三个双引号所构造的字符串有两个功能性的区别。
#
# 一个区别是，三个双引号所构造的字符串中可以使用双引号(`"`)作为字符串的一部分。
# ------------------------------------------------------------------------------------------

"Here, we get an "error" because it's ambiguous where this string ends "

"""Look, Mom, no "errors"!!! """

# ------------------------------------------------------------------------------------------
# 注意 `' '` 定义的是字符（character），不是字符串！
# ------------------------------------------------------------------------------------------

typeof('a')

'We will get an error here'

# ------------------------------------------------------------------------------------------
# ## 格式化字符串（String interpolation）
#
# `$` 符号可以将已有的变量插入一个字符串，也可以在字符串中代入表达式的值。
#
# 下面这个例子含有高度敏感的个人信息。
# ------------------------------------------------------------------------------------------

name = "Jane"
num_fingers = 10
num_toes = 10

println("Hello, my name is $name.")
println("I have $num_fingers fingers and $num_toes toes.")

 println("That is $(num_fingers + num_toes) digits in all!!")

# ------------------------------------------------------------------------------------------
# ## 字符串拼接
#
# 下面介绍几种拼接字符串的方法。
#
# 第一种方式是使用 `string()` 函数。
# `string()` 函数将非字符串的输入转化为字符串。
# ------------------------------------------------------------------------------------------

s3 = "How many cats ";
s4 = "is too many cats?";
😺 = 10

string(s3, s4)

string("I don't know, but ", 😺, " is too few.")

# ------------------------------------------------------------------------------------------
# 也可以使用 `*` 进行拼接！
# ------------------------------------------------------------------------------------------

s3*s4

# ------------------------------------------------------------------------------------------
# ## 练习
#
# ### 2.1
# 创建一个有1000个"hi"的字符串 `hi`，先用 `repeat` 函数，在试试指数运算符，它会在后台调用 `*`。
# ------------------------------------------------------------------------------------------

# 在此作答
hi = "hi"

@assert hi == "hihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihi"

# ------------------------------------------------------------------------------------------
# ### 2.2
# 声明两个变量
#
# ```julia
# a = 3
# b = 4
# ```
# 然后用这两个变量创建两个字符串：
# ```julia
# "3 + 4"
# "7"
# ```
# 将这两个字符串分别赋值给 `c` 和 `d`
# ------------------------------------------------------------------------------------------

# 在此作答
c = ""
d = ""

@assert c == "3 + 4"
@assert d == "7"

# ------------------------------------------------------------------------------------------
# 请在完成练习运行 `@assert` 块以检测结果是否正确，若无报错即为正确。
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
using SellmeierFit
using SimpleConstants  # for micro = 1e-6
using CairoMakie

## Main section
# Load a set of measured refractive index data from https://RefractiveIndex.info.
path = "https://refractiveindex.info/tmp/data/main/SiO2/Malitson.csv"
(; λ, ε) = SellmeierFit.read(path)

# Fit the loaded data to the Sellmeier equation.
(; mdl, err) = fit_sellmeier(λ, ε)

## Analysis section
# The calculated fit parameters are the same as those shown in RefractiveIndex.info's
# Dispersion formula section.
println("B = $(mdl.str)")
println("√C (µm) = $(mdl.λres ./ micro)")
println("Error between data and fit = $err")

# Visualize the measured data and fit Sellmeier equation; note that n = √ε.
fontsize = 20
fig = Figure(; fontsize)
Axis(fig[1,1], title="Refractive Index of SiO₂", xlabel="λ (meter)", ylabel="n")
scatter!(λ, .√ε, label="Measurement", color=:black)
lines!(λ, .√mdl.(λ), label="Sellmeier equation", color=:red)
axislegend(; position=:rt)
display(fig)


text:
using YaoPlots, ZXCalculus, Graphs
using Test

@testset "helperblock" begin
    include("helperblock.jl")
end

@testset "vizcircuit" begin
    include("vizcircuit.jl")
end

@testset "zx_plot" begin
    include("zx_plot.jl")
end


text:
#Lax shock tube
using ConservationLawsDiffEq
using OrdinaryDiffEq

#Euler Equations
const CFL = 0.5
const Tend = 0.2
const γ=1.4 #gas constant
const xdiafragm = 0.0

function Jf(u::AbstractVector)
  ρ = u[1]; v = u[2]/u[1]; ϵ=u[3]
  p = (ϵ-0.5*ρ*v^2)*(γ-1)
  F =[0.0 1.0 0.0;-v^2*(1+γ)/2 v*(3-γ) (γ-1);v^3*(γ-1)+γ*ϵ*v/ρ 3/2*v^2*(1-γ)+γ*ϵ/ρ γ*v]
  F
end

function f(u::AbstractVector)
  ρ = u[1]; v = u[2]/u[1]; ϵ=u[3]
  p = (ϵ-0.5*ρ*v^2)*(γ-1)
  [u[2];u[2]^2/u[1]+p;(ϵ+p)*v]
end

#initial condition
function f0(x::Real)
  ufis_l = [1.0, 0.0,	1.0]
  ufis_r = [0.125 0.0 0.1]
  if x <= xdiafragm
    ρ,u,p = ufis_l
  else
    ρ,u,p = ufis_r
  end
  return [ρ, ρ*u, 1.0/(γ - 1.0)*p + 0.5*ρ*u*u]
end

positive(x) = max(0,x)
#define max wave speed (generic euler system)
function max_w_speed(u)
  ρ = u[:,1];v = u[:,2]./ρ;E = u[:,3]
  p = fill!(similar(ρ), zero(eltype(ρ))); wave_speed = fill!(similar(ρ), zero(eltype(ρ)))
  @. p = (γ-1)*(E-0.5*ρ*v*v)
  @. wave_speed = abs(v) + sqrt(γ*positive(p/ρ))
  return maximum(wave_speed)
end

function Nflux(ul::Vector, ur::Vector)
  ρl = ul[1]; vl = ul[2]/ul[1]; ϵl=ul[3]
  pl = (ϵl-0.5*ρl*vl^2)*(γ-1)
  ρr = ur[1]; vr = ur[2]/ur[1]; ϵr=ur[3]
  pr = (ϵr-0.5*ρr*vr^2)*(γ-1)
  zl = sqrt(ρl/pl)*[1;vl;pl]
  zr = sqrt(ρr/pr)*[1;vr;pr]
  zm = 0.5*(zl+zr)
  zln = (zr-zl)./(log(zr)-log(zl))
  F = fill!(similar(ul), zero(eltype(ul)))
  F[1] = zm[2]*zln[3]
  F[2] = zm[3]/zm[1]+zm[2]/zm[1]*F[1]
  F[2] = 0.5*zm[2]/zm[1]*((γ+1)/(γ-1)*zln[3]/zln[1]+F[2])
  F
end

function ve(u::Vector)
  ρ = u[1]; v = u[2]/u[1]; ϵ=u[3]
  p = (ϵ-0.5*ρ*v^2)*(γ-1)
  s = log(p)-γ*log(ρ)
  return [(γ-s)/(γ-1)-ρ*v^2/(2*p);ρ*v/p;-ρ/p]
end

# Now discretizate the domain
mesh = Uniform1DFVMesh(20, [-1.0, 1.0])

function get_problem(alg, mesh)
  #Compute discrete initial data
  u0 = getInitialState(mesh,f0,use_threads = true)

  # Now get a explicit semidiscretization (discrete in space) du_h(t)/dt = f_h(u_h(t))
  f_h = getSemiDiscretization(f,alg,mesh,[ZeroFlux()]; Df = Jf, use_threads = false,numvars = 3)

  #Setup ODE problem for a time interval = [0.0,1.0]
  ode_prob = ODEProblem(f_h,u0,(0.0,Tend))

  #Setup callback in order to fix CFL constant value
  cb = getCFLCallback(f_h, CFL)

  #Estimate an initial dt
  dt = update_dt!(u0, f_h, CFL)
  return ode_prob, cb, dt
end

ode_prob, cb, dt = get_problem(FVSKTScheme(), mesh)
sol = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

ode_prob, cb, dt = get_problem(FVCompWENOScheme(), mesh)
sol2 = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

u1_h = fv_solution(sol, mesh; vars = 3)
u2_h = fv_solution(sol2, mesh; vars = 3)

#Plot
using Plots
plot(u1_h, tidx=1, vars=1, lab="ρ_0",line=(:dot,2))
plot!(u1_h, vars=1, lab="KT ρ",line = (:dot,2))
plot!(u2_h, vars=1, lab="CompWENO5 ρ",line = (:dot,2))


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: importexport
using AWS.Compat
using AWS.UUIDs

"""
    CancelJob()

This operation cancels a specified job. Only the job owner can cancel it. The operation fails if the job has already started or is complete.

# Required Parameters
- `JobId`:

# Optional Parameters
- `APIVersion`:
"""
cancel_job(JobId; aws_config::AWSConfig=global_aws_config()) = importexport("CancelJob", Dict{String, Any}("JobId"=>JobId); aws_config=aws_config)
cancel_job(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = importexport("CancelJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobId"=>JobId), args)); aws_config=aws_config)

"""
    CreateJob()

This operation initiates the process of scheduling an upload or download of your data. You include in the request a manifest that describes the data transfer specifics. The response to the request includes a job ID, which you can use in other operations, a signature that you use to identify your storage device, and the address where you should ship your storage device.

# Required Parameters
- `JobType`:
- `Manifest`:
- `ValidateOnly`:

# Optional Parameters
- `APIVersion`:
- `ManifestAddendum`:
"""
create_job(JobType, Manifest, ValidateOnly; aws_config::AWSConfig=global_aws_config()) = importexport("CreateJob", Dict{String, Any}("JobType"=>JobType, "Manifest"=>Manifest, "ValidateOnly"=>ValidateOnly); aws_config=aws_config)
create_job(JobType, Manifest, ValidateOnly, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = importexport("CreateJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobType"=>JobType, "Manifest"=>Manifest, "ValidateOnly"=>ValidateOnly), args)); aws_config=aws_config)

"""
    GetShippingLabel()

This operation generates a pre-paid UPS shipping label that you will use to ship your device to AWS for processing.

# Required Parameters
- `jobIds`:

# Optional Parameters
- `APIVersion`:
- `city`:
- `company`:
- `country`:
- `name`:
- `phoneNumber`:
- `postalCode`:
- `stateOrProvince`:
- `street1`:
- `street2`:
- `street3`:
"""
get_shipping_label(jobIds; aws_config::AWSConfig=global_aws_config()) = importexport("GetShippingLabel", Dict{String, Any}("jobIds"=>jobIds); aws_config=aws_config)
get_shipping_label(jobIds, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = importexport("GetShippingLabel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("jobIds"=>jobIds), args)); aws_config=aws_config)

"""
    GetStatus()

This operation returns information about a job, including where the job is in the processing pipeline, the status of the results, and the signature value associated with the job. You can only return information about jobs you own.

# Required Parameters
- `JobId`:

# Optional Parameters
- `APIVersion`:
"""
get_status(JobId; aws_config::AWSConfig=global_aws_config()) = importexport("GetStatus", Dict{String, Any}("JobId"=>JobId); aws_config=aws_config)
get_status(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = importexport("GetStatus", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobId"=>JobId), args)); aws_config=aws_config)

"""
    ListJobs()

This operation returns the jobs associated with the requester. AWS Import/Export lists the jobs in reverse chronological order based on the date of creation. For example if Job Test1 was created 2009Dec30 and Test2 was created 2010Feb05, the ListJobs operation would return Test2 followed by Test1.

# Optional Parameters
- `APIVersion`:
- `Marker`:
- `MaxJobs`:
"""
list_jobs(; aws_config::AWSConfig=global_aws_config()) = importexport("ListJobs"; aws_config=aws_config)
list_jobs(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = importexport("ListJobs", args; aws_config=aws_config)

"""
    UpdateJob()

You use this operation to change the parameters specified in the original manifest file by supplying a new manifest file. The manifest file attached to this request replaces the original manifest file. You can only use the operation after a CreateJob request but before the data transfer starts and you can only use it on jobs you own.

# Required Parameters
- `JobId`:
- `JobType`:
- `Manifest`:
- `ValidateOnly`:

# Optional Parameters
- `APIVersion`:
"""
update_job(JobId, JobType, Manifest, ValidateOnly; aws_config::AWSConfig=global_aws_config()) = importexport("UpdateJob", Dict{String, Any}("JobId"=>JobId, "JobType"=>JobType, "Manifest"=>Manifest, "ValidateOnly"=>ValidateOnly); aws_config=aws_config)
update_job(JobId, JobType, Manifest, ValidateOnly, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = importexport("UpdateJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobId"=>JobId, "JobType"=>JobType, "Manifest"=>Manifest, "ValidateOnly"=>ValidateOnly), args)); aws_config=aws_config)


text:
#production and acquisition

function DDCost!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    key = (options.postprocess ? "postprocessing" : "evaluating")
    println(key*" decommissioning cost")
    decommissioning!(options,sizes,data,output)
    println(key*" wasteManagement cost")
    wasteManagement!(options,sizes,data,output)
    println(key*" siteClearance cost")
    siteClearance!(options,sizes,data,output)
    println(key*" postMonitoring cost")
    postMonitoring!(options,sizes,data,output)

    nothing
end

function decommissioning!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.decommissioning == 0
        if options.postprocess
            nothing
        else
            CDDport = (data.Cddport + data.Nddport .* data.Lddport).*transpose(data.turbineQuantity)
            Cremov = (data.Nremove .* data.Vremove).*transpose(data.turbineQuantity)
            output.DECOM .+= CDDport + Cremov
            output.decommissioning .+= CDDport + Cremov
        end
    else
        error("decommissioning mode not found.")
    end
    nothing

end

function wasteManagement!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.wasteManagement == 0
        if options.postprocess
            nothing
        else
            Cwproc = data.Wjproc .* data.Cjproc
            Cwtrans = data.Wjtrans .* data.Cjproc ./ data.Wttrans
            Cland = data.Wnr .* data.Cnr
            Csv = data.Wr .* data.SV
            output.DECOM .+= transpose(Cwproc + Cwtrans + Cland - Csv).*transpose(data.turbineQuantity)
            output.wasteManagement .+= transpose(Cwproc + Cwtrans + Cland - Csv).*transpose(data.turbineQuantity)
        end
    else
        error("wasteManagement mode not found.")
    end
    nothing

end

function siteClearance!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.siteClearance == 0
        if options.postprocess
            nothing
        else
            output.DECOM .+= data.Asc.*data.Cscunit.*transpose(data.turbineQuantity)
            output.siteClearance .+= data.Asc.*data.Cscunit.*transpose(data.turbineQuantity)
        end
    else
        error("siteClearance mode not found.")
    end
    nothing

end

function postMonitoring!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.postMonitoring == 0
        if options.postprocess
            nothing
        else
            output.DECOM .+= data.postCost
            output.postMonitoring .+= data.postCost
        end
    else
        error("postMonitoring mode not found.")
    end
    nothing

end

max_stars_repo_name: Thiago-NovaesB/zefiro

max_issues_repo_name: Thiago-NovaesB/zefiro

max_forks_repo_name: Thiago-NovaesB/zefiro


text:
module CompressedSensing

# package code goes here
using LinearAlgebra
include("IRLS.jl")
include("UIRLS.jl")
#include("GiniIndex.jl")
#include("nGMCA.jl")
#include("ZAP.jl")
#include("Coherence.jl")

#export IRLS,UIRLS,GI,nGMCA,ZAP,Coherence,CSExamples
export IRLS, UIRLS

#function CSExamples()
#	include("../examples/Example 1.jl")
#
#end

end # module


text:
#!/usr/bin/env julia
#
# usage:
#  script......... plaint.jl /path/to/system.yml
#  interactively.. julia -qi plaint.jl

using YAML              #read data file
using DelimitedFiles    #output results
using LinearAlgebra

"""
    ksnr(e, M, eps, max_iterations)

Solve Kepler's equation with NR method.
"""
function ksnr(e::Float64, M::Float64,
              eps::Float64 = 1e-16, nm::Int64 = 16)
    x = M + e * sin(M)
    d = (x - e * sin(x) - M) / (1 - e * cos(x))
    while abs(d) > eps && nm > 0
        x -= d
        d = (x - e * sin(x) - M) / (1 - e * cos(x))
        nm -= 1
    end
    return x
end

"""
    euler_rotd(Ω, i, ω)

Return Eulerian rotation.
"""
function euler_rotd(Ω::Float64, i::Float64, ω::Float64)
    return [ cosd(Ω) -sind(Ω) 0; sind(Ω) cosd(Ω) 0; 0 0 1 ] *
           [ 1 0 0; 0 cosd(i) -sind(i); 0 sind(i) cosd(i) ] *
           [ cosd(ω) -sind(ω) 0; sind(ω) cosd(ω) 0; 0 0 1 ]
end

"""
    oe_to_xv!(mu, j, orbital_elements, positions, velocities)

Convert orbital elements to cartesian vectors.
"""
function oe_to_xv!(mu::Float64, j::Int64,
                   oe::Array{Float64,2},
                   vr::Array{Float64,2}, vv::Array{Float64,2})
    a, e, i, Ω, ω, M = oe[:, j]
    E = ksnr(e, M * π / 180)
    r = a * (1 - e * cos(E))
    pof = [cos(E) - e, sqrt(1 - e^2) * sin(E), 0] * a
    vof = [  - sin(E), sqrt(1 - e^2) * cos(E), 0] * (sqrt(mu * a) / r)
    rot = euler_rotd(Ω, i, ω)
    vr[:, j] = rot * pof
    vv[:, j] = rot * vof
end

"""
    xv_to_oe!(mu, j, positions, velocities, orbital_elements)

Convert cartesian vectors to orbital elements.
"""
function xv_to_oe!(mu::Float64, j::Int64,
                   rt::Array{Float64,2}, vt::Array{Float64,2},
                   oe::Array{Float64,2})
    tiny = 1e-10
    vr, vv = rt[:, j], vt[:, j]
    r = norm(vr)
    v = norm(vv)
    vh = cross(vr, vv)
    h = norm(vh)
    vn = [-vh[2], vh[1], 0]
    n = norm(vn)
    ve = cross(vv, vh) / mu - vr / r
    a = 1 / (2 / r - v^2 / mu)
    e = norm(ve)
    i = acos(vh[3] / h)
    if e > tiny
        if n > tiny
            Ω = acos(vn[1] / n)
            if vn[2] < 0
                Ω = 2π - Ω
            end
            ω = acos(dot(ve, vn) / (e * n))
            if ve[3] < 0
                ω = 2π - ω
            end
        else
            i = Ω = .0
            ω = atan2(ve[2], ve[3])
            if cross(ve, vr)[3] < 0
                ω = 2π - ω
            end
        end
        nu = acos(dot(ve, vr)/(e * r))
        if dot(vr, vv) < 0
            nu = 2π - nu
        end
    else
        #not implemented
    end
    E = 2 * atan(sqrt((1 - e) / (1 + e)) * tan(nu / 2))
    M = E - e * sin(E)
    if M < 0
        M += 2π
    end
    i, Ω, ω, M = [i, Ω, ω, M] * 180 / π
    oe[:, j] = [a, e, i, Ω % 360, ω % 360, M % 360]
end

"""
    grav_a!(no_of_bodies, masses, positions, acceleration)

Compute gravitational acceleration.
!!! note
    The `a` should be multiplied with the appropriate gravitational constant.
"""
function grav_a!(n::Int64, m::Array{Float64,1},
                 x::Array{Float64,2}, a::Array{Float64,2})
    fill!(a, 0)
    @inbounds for i = 1:(n - 1), j = (i + 1):n
        R = x[:, i] - x[:, j]
        Rd3 = R / norm(R)^3
        a[:, i] -= m[j] * Rd3
        a[:, j] += m[i] * Rd3
    end
end

"""
    mksnr(dt, α, n, σ0, r0, eps, max_iterations)

Solve modified Kepler's equation with NR method.
!!! note
    Argument `α` is the inverse of semimajor axis.
"""
function mksnr(dt::Float64, α::Float64, n::Float64, s::Float64, r0::Float64,
               eps::Float64 = 1e-16, nm::Int64 = 16)
    r0a, s0a, ndt = 1 - r0 * α, s * sqrt(α), n * dt
    x = π
    d = (x - r0a * sin(x) - s0a * (cos(x) - 1) - ndt) /
            (1 - r0a * cos(x) + s0a * sin(x))
    while abs(d) > eps && nm > 0
        x -= d
        d = (x - r0a * sin(x) - s0a * (cos(x) - 1) - ndt) /
                (1 - r0a * cos(x) + s0a * sin(x))
        nm -= 1
    end
    return x
end

"""
    gauss_step!(dt, mu, no_of_bodies, positions, velocities)

Keplerian drift using Gauss' f, g functions.
"""
function gauss_step!(dt::Float64, mu::Float64, np::Int64,
                     rt::Array{Float64,2}, vt::Array{Float64,2})
    @inbounds @simd for i = 1:np
        vr0, vv0 = rt[:, i], vt[:, i]
        r0, v0 = norm(vr0), norm(vv0)
        α = 2 / r0 - v0^2 / mu
        a, n = 1 / α, sqrt(mu * α^3)
        s = dot(vr0, vv0) / sqrt(mu)
        E = mksnr(dt, α, n, s, r0)
        cose, sine = cos(E), sin(E)
        r1 = a + (r0 - a) * cose + sqrt(a) * s * sine
        a1ce = a * (cose - 1)
        ivr0, ivr1 = 1 / r0, 1 / r1
        ft = 1 + a1ce * ivr0
        gt = dt + (sine - E) / n
        fd = - a^2 * n * sine * ivr1 * ivr0
        #gd = 1 + a1ce * ivr1
        gd = (1 + gt * fd) / ft
        rt[:, i] = ft * vr0 + gt * vv0
        vt[:, i] = fd * vr0 + gd * vv0
    end
end

"""
    eam_dh(no_of_bodies, masses, positions, velocities)

Compute energy and angular momentum, with democratic heliocentric vectors.
"""
function eam_dh(np::Int64, m::Array{Float64,1},
                x::Array{Float64,2}, v::Array{Float64,2})
    et, am = norm(sum(m' .* v, dims = 2))^2 / 2, zeros(3)
    @inbounds @simd for i = 1:np
        et += m[i] * norm(v[:, i])^2 / 2
        et -= m[i] / norm(x[:, i])
        am += m[i] * cross(x[:, i], v[:, i])
    end
    @inbounds for i = 1:(np - 1), j = (i + 1):np
        et -= m[i] * m[j] / norm(x[:, i] - x[:, j])
    end
    return et, norm(am)
end

"""
    plaint_dhi(iostream, total_iterations, dt, step_iterations,
               mu, masses, orbital_elements)

Perform planetary integration using democratic heliocentric integrator.
"""
function plaint_dhi(of::Array{IOStream, 1},
                    kmax::Float64, dt::Float64, kpri::Float64, mu::Float64,
                    mp::Array{Float64,1}, oe::Array{Float64,2})
    mt = 1 / (sum(mp) + 1)
    np = size(mp, 1)
    an = Array{Float64}(undef, 3, np)
    x = similar(an)
    v = similar(an)
    @inbounds @simd for j = 1:np
        oe_to_xv!(mu, j, oe, x, v)
    end
    v .-= mt * sum(mp' .* v, dims = 2)                  #to_dh
    el0 = eam_dh(np, mp, x, v)
    println("Initial energy:           ", el0[1], "\n",
            "Initial angular momentum: ", el0[2], "\n")
    d2 = dt / 2
    @inbounds for k = 1:kmax
        try
            x .+= d2 * sum(mp' .* v, dims = 2)
            gauss_step!(d2, mu, np, x, v)
            grav_a!(np, mp, x, an)
            v .+= dt * an
            gauss_step!(d2, mu, np, x, v)
            x .+= d2 * sum(mp' .* v, dims = 2)
            if k % kpri == 0
                eln = eam_dh(np, mp, x, v)
                va = v .+ mt * sum(mp' .* v, dims = 2)  #to_ac
                writedlm(of[end], [collect(eln .- el0); collect(@. 1 - eln / el0)]')
                @simd for j = 1:np
                    xv_to_oe!(mu, j, x, va, oe)
                    writedlm(of[j], oe[:, j]')
                end
            end
        catch
            println("Error: Exiting on ", k, "\n")
            break
        end
    end
    close.(of)
end

"""
    data_get(file)

Get planetary system's data and integration options.
"""
function data_get(file::String)
    data = YAML.load_file(file)
    tmax::Float64 = data["options"]["tmax"]
    step::Float64 = data["options"]["step"]
    ever::Float64 = data["options"]["ever"]
    unit = "solar"
    try
        unit = data["system"]["units"]
    catch
    end
    ms = data["system"]["star_mass"]
    datp = data["system"]["planets"]
    sort!(datp, by = x -> x["orbital_elements"][1])
    np = size(datp, 1)
    name = Array{String}(undef, np)
    mp = Array{Float64}(undef, np)
    oe = Array{Float64}(undef, 6, np)
    @inbounds @simd for i = 1:np
        name[i] = data["system"]["star_name"] * " " * datp[i]["id"]
        mp[i] = datp[i]["mass"]
        oe[:, i] = datp[i]["orbital_elements"]
    end
    return tmax, step, ever, ms, name, mp, oe, unit
end

"""
    data_prep(tmax, step, ever,
              star_mass, file_names, masses, orbital_elements, units)

Prepare data and options to feed to the integrator.
"""
function data_prep(tmax::Float64, step::Float64, ever::Float64,
                   ms::Float64, name::Array{String, 1},
                   mp::Array{Float64, 1}, oe::Array{Float64,2},
                   unit::String = "solar")
    if unit == "solar"
        mp = mp * 9.547919e-4 / ms
    else #natural
        mp = mp / ms
    end
    oe[1, :] /= oe[1, 1]
    tp = 2π
    dt = step * tp
    kmax = div(tmax * tp, dt)
    kpri = div(ever * tp, dt)
    of = map(f -> open("out/$f.txt", "w"), name)
    push!(of, open("out/eamerr.txt", "w"))
    println("Total simulation time is: ", tmax, " yr\n",
            "Time of every step to be: ", step, " yr\n",
            "Total number of steps is: ", kmax, " steps\n",
            "Output will be every:     ", ever, " yr\n",
            "or differently every:     ", kpri, " steps\n")
    return of, kmax, dt, kpri, 1.0, mp, oe
end

"""
    usage()

Return an exemplary session.
"""
function usage()
    println("## Example session:\n",
            "#set file containing data\n",
            "file = \"data/TRAPPIST-1.yml\";\n",
            "#get data from file\n",
            "tmax, step, ever, ms, name, mp, oe = data_get(file);\n",
            "#or enter manually\n",
            "dgr = (tmax, step, ever, ms, name, mp, oe);\n",
            "#prepare data\n",
            "of, kmax, dt, kpri, mu, mp, oe = data_prep(dgr...);\n",
            "#integrate, available: plaint_dhi, plaint_s6i\n",
            "plaint_dhi(of, kmax, dt, kpri, mu, mp, oe);\n",
            "#benchmark integration of a particular system\n",
            "@time plaint_dhi(data_prep(data_get(file)...)...);")
end

## Session
# Run as script or, interactively in a Julia session.

if ! isinteractive()
    @time plaint_dhi(data_prep(data_get(ARGS[1])...)...)
else
    usage()
end


text:
###############################
# @dfn, @req, @prp, @mth, @sprp
###############################
#=
(Example)

@prp tp = x -> Base.typeof(x)
tp = @prp x -> Base.typeof(x)

[1,2,3].tp == Array{Int64,1}

@mth sz = (x,i) -> size(x,i)
sz = @mth (x,i) -> size(x,i)
[1,2,3].sz(1) == 3

g = NS()
g.a = 3
@prp g.b = g -> 3 * g.a
g.c = @prp g -> 3 * g.b
@prp g.d(g) = 3 * g.c

@mth g.f = (g,x) ->  g.a + x
g.g = @mth (g,x) ->  g.b + x
@mth g.h(g,x) =  g.c + x
=#

macro prp(ex)
    # f = @prp x -> 2x
    if ex.head != :(=)
        return esc(:(Wild.NSTagFunc{:prp}($(ex))))
    end

    #=
    # @prp f(x) = 2x
    if hasfield(typeof(ex.args[1]), :head) && ex.args[1].head == :call
        name = ex.args[1].args[1]
        ex.args[1].args[1] = gensym()
        ex2 = :($(name) = Wild.NSTagFunc{:prp}($(ex.args[1].args[1])))
        ex = Meta.parse(string(ex) * ";" * string(ex2))
        return esc(ex)
    end
    =#

    # @prp f = x -> 2x
    return esc(:($(ex.args[1]) = Wild.NSTagFunc{:prp}($(ex.args[2]))))
end

macro mth(ex)
    # f = @mth (x,y) -> x+y
    if ex.head != :(=)
        return esc(:(Wild.NSTagFunc{:mth}($(ex))))
    end

    #=
    # @mth f = (x,y) -> x+y
    if hasfield(typeof(ex.args[1]), :head) && ex.args[1].head == :call
        name = ex.args[1].args[1]
        ex.args[1].args[1] = gensym()
        ex2 = :($(name) = Wild.NSTagFunc{:mth}($(ex.args[1].args[1])))
        ex = Meta.parse(string(ex) * ";" * string(ex2))
        return esc(ex)
    end
    =#

    # @mth f(x,y) = x+y
        return esc(:($(ex.args[1]) = Wild.NSTagFunc{:mth}($(ex.args[2]))))
end

###############################
# @prpfnc, @mthfnc
###############################
abstract type AbstTagFunc <: Function end
abstract type AbstPrpFunc <: AbstTagFunc end
abstract type AbstMthFunc <: AbstTagFunc end

"""
    @prpfnc fncname

# Examples
@prpfnc dtype
(::dtype.type)(itr) = eltype(itr)
(::dtype.type)(itr::Base.Generator) = Base.return_types(itr.f, (dtype(itr.iter),))[1]
```
"""
macro prpfnc(name)
    typename = Base.gensym()
    ex = quote
        struct $(typename)  <: AbstPrpFunc
            type::Type
        end
        $name = ($typename)(:($$typename))
    end
    esc(ex)
end

"""
    @mthfnc fncname

# Examples
@mthfnc astype
(::astype.type)(m::AbstractArray, t::Type) =
    (n = Base.similar(m, t, m.size()...); n .= m; n)
(::astype.type)(x, t::Type) = convert(t, x)
```
"""
macro mthfnc(name)
    typename = Base.gensym()
    ex = quote
        struct $(typename)  <: AbstMthFunc
            type::Type
        end
        $name = ($typename)(:($$typename))
        (f::$typename)(x) = (a...; ka...) -> f(x, a...; ka...)
    end
    esc(ex)
end


text:
module TestCore

using Test
using MLJIteration
using IterationControl
using MLJBase
using MLJModelInterface
using ..DummyModel

X, y = make_dummy(N=20)
model = DummyIterativeModel(n=0)

@testset "integration: resampling=nothing" begin

    model = DummyIterativeModel(n=0)
    controls=[Step(2), Threshold(0.01), TimeLimit(0.005)]

    # using IterationControl.jl directly:
    mach = machine(deepcopy(model), X, y);
    function IterationControl.train!(mach::Machine{<:DummyIterativeModel},
                                     n::Int)
        mlj_model = mach.model
        mlj_model.n = mlj_model.n + n
        fit!(mach, verbosity=0)
    end
    IterationControl.loss(mach::Machine{<:DummyIterativeModel}) =
        last(training_losses(mach))
    IterationControl.train!(mach, controls..., verbosity=0)
    losses1 = report(mach).training_losses
    yhat1 = predict(mach, X)

    # using IteratedModel wrapper:
    imodel = IteratedModel(model=model,
                           controls=controls,
                           resampling=nothing)
    mach = machine(imodel, X, y)
    @test_logs((:info, r"Training"),
               (:info, MLJIteration.info_unspecified_iteration_parameter(:n)),
               (:info, r"final loss"),
               (:info, r"final training loss"),
               (:info, r"Stop"),
               (:info, r"Total"),
               fit!(mach, verbosity=1))
    imodel.resampling = nothing
    @test_logs fit!(mach, verbosity=0)

    losses2 = report(mach).model_report.training_losses
    yhat2 = predict(mach, X);
    @test imodel.model == DummyIterativeModel(n=0) # hygeine check

    # compare:
    @test losses1 ≈ losses2
    @test yhat1 ≈ yhat2

    # report:
    r = report(mach)
    @test r.model_report.training_losses ≈ losses1
    @test :controls in keys(r)
    i = r.n_iterations
    @test i == 38

    # warm restart when changing controls:
    noise(n) = fill((:info, r"43"), n)
    imodel.controls = [Step(1),
                       WithIterationsDo(i-> i>39,
                                        stop_if_true=true),
                       Info(x->"43")]
    @test_logs((:info, r"Updating"),
               noise(2)...,
               (:info, r""),
               (:info, r""),
               (:info, r""),
               (:info, r""),
               fit!(mach))
    @test report(mach).n_iterations == i + 2

    # warm restart when changing iteration parameter (trains one more
    # iteration because stopping control comes after `Step(...)`):
    imodel.model.n = 1
    @test_logs((:info, r"Updating"),
               noise(1)...,
               (:info, r""),
               (:info, r""),
               (:info, r""),
               (:info, r""),
               fit!(mach))
    @test report(mach).n_iterations == i + 3

    # cold restart when changing anything else:
    imodel.check_measure=false
    imodel.controls = [Step(1),
                       WithIterationsDo(i-> i>4,
                                        stop_if_true=true),
                       Info(x->"43")]
    @test_logs((:info, r"Updating"),
               (:info, MLJIteration.info_unspecified_iteration_parameter(:n)),
               noise(5)...,
               (:info, r""),
               (:info, r""),
               (:info, r""),
               (:info, r""),
               fit!(mach))
    @test report(mach).n_iterations == 5
end

@testset "integration: resampling=Holdout()" begin

    controls=[Step(2), Patience(4), TimeLimit(0.001)]

    # using IterationControl.jl directly:
    mach = machine(deepcopy(model), X, y);
    train, test = partition(eachindex(y), 0.7)
    function IterationControl.train!(mach::Machine{<:DummyIterativeModel},
                                     n::Int)
        mlj_model = mach.model
        mlj_model.n = mlj_model.n + n
        fit!(mach, rows=train, verbosity=0)
    end
    function IterationControl.loss(mach::Machine{<:DummyIterativeModel})
        yhat = predict(mach, rows=test)
        return mae(yhat, y[test]) |> mean
    end
    IterationControl.train!(mach, controls...; verbosity=0)
    losses1 = report(mach).training_losses;
    yhat1 = predict(mach, X[test]);
    niters = mach.model.n
    @test niters == length(losses1)

    # using IteratedModel wrapper:
    imodel = IteratedModel(model=model,
                           resampling=Holdout(fraction_train=0.7),
                           controls=controls)
    mach = machine(imodel, X, y)
    @test_logs((:info, r"Training"),
               (:info, MLJIteration.info_unspecified_iteration_parameter(:n)),
               (:info, MLJIteration.info_unspecified_measure(rms)),
               (:info, r"final loss"),
               (:info, r"final train"),
               (:info, r"Stop"),
               (:info, r"Total"),
               fit!(mach, verbosity=1))
    imodel.measure = mae
    @test_logs fit!(mach, verbosity=0)
    losses2 = report(mach).model_report.training_losses;
    yhat2 = predict(mach, X[test]);

    # compare:
    @test losses1 ≈ losses2
    @test yhat1 ≈ yhat2

    # now repeat wrapper train with retrain=true:
    imodel = IteratedModel(model=model,
                           resampling=Holdout(fraction_train=0.7),
                           controls=controls,
                           retrain=true,
                           measure=mae)
    mach = machine(imodel, X, y)
    fit!(mach, verbosity=0)
    losses3 = report(mach).model_report.training_losses
    yhat3 = predict(mach, X[test])
    @test fitted_params(mach).machine.model.n == niters

    # to be compared with atomic model trained on all data for same
    # number of iterations:
    model2 = DummyIterativeModel(n=niters)
    mach2 = machine(model2, X, y)
    fit!(mach2, verbosity=0)
    losses4 = report(mach2).training_losses
    yhat4 = predict(mach2, X[test])

    # check report:
    r = report(mach2)

    # compare:
    @test losses3 ≈ losses4
    @test yhat3 ≈ yhat4
end

@testset "integration: cyclic learning rates" begin

    X, y = make_dummy(N=100)
    ss = 3
    model = DummyIterativeModel()

    # constant learning rate:
    controls = [CycleLearningRate(stepsize=ss,
                                  lower=1,
                                  upper=1,
                                  learning_rate_parameter=:learning_rate),
                Step(1),
                NumberLimit(5*ss)]

    imodel = IteratedModel(model=model,
                           resampling=Holdout(fraction_train=0.7),
                           controls=controls,
                           measure=mae)

    mach = machine(imodel, X, y)
    fit!(mach, verbosity=0)
    yhat1 = predict(mach, X)

    # no cycling:
    imodel.controls = [Step(1), NumberLimit(5*ss)]
    fit!(mach, force=true, verbosity=0)
    yhat2 = predict(mach, X)

    # with real cycling:
    controls = [CycleLearningRate(stepsize=ss,
                                  lower=0.5,
                                  upper=1.5,
                                  learning_rate_parameter=:learning_rate),
                Step(1),
                NumberLimit(5*ss)]
    imodel.controls = controls
    fit!(mach, force=true, verbosity=0)
    yhat3 = predict(mach, X)

    # compare:
    @test yhat1 ≈ yhat2
    @test !(yhat3 ≈ yhat2)

end

@testset "issue #40" begin
    # change trait to an incorrect value:
    MLJModelInterface.iteration_parameter(::Type{<:DummyIterativeModel}) = :junk
    @test iteration_parameter(model) == :junk

    # user specifies correct value:
    iterated_model = IteratedModel(model, measure=mae, iteration_parameter=:n)
    mach = machine(iterated_model, X, y)

    # and model still runs:
    fit!(mach, verbosity=0)

    # change trait back:
    MLJModelInterface.iteration_parameter(::Type{<:DummyIterativeModel}) = :n
    @test iteration_parameter(model) == :n
end

end
true


text:
function show()
    myfunc = (x, y) -> x + y
    println(myfunc(1, 2))
end

if abspath(PROGRAM_FILE) == @__FILE__
    show()
end


text:
using Random

#=
length(::WeirdType{T,N}) where {T,N} = N
Dual(1, svec(1,2,3)) # rocks
up = SVector
Dual(1, partials) --> Dual(1, SVector(partials...))
Dual(1, partials::StaticArrays)
Dual(1, [1], y, z) --> Dual(1, [[1], y, z])
Dual(1, [1]) --> Dual(1, 1)
Dual(1, [1]) --> Dual(1, Partials([1,])) # fail Any

Dual(::Number1, [1]) --> Dual(1, Partials([1,])) # fail Any
=#

# Modified from `ForwardDiff.jl`

###########
# Prelude #
###########

const NANSAFE_MODE_ENABLED = false

const UNARY_PREDICATES = Symbol[:isinf, :isnan, :isfinite, :iseven, :isodd, :isreal, :isinteger]

const DEFAULT_CHUNK_THRESHOLD = 12

struct Chunk{N} end

const CHUNKS = [Chunk{i}() for i in 1:DEFAULT_CHUNK_THRESHOLD]

function Chunk(input_length::Integer, threshold::Integer = DEFAULT_CHUNK_THRESHOLD)
    N = pickchunksize(input_length, threshold)
    N <= DEFAULT_CHUNK_THRESHOLD && return CHUNKS[N]
    return Chunk{N}()
end

function Chunk(x::AbstractArray, threshold::Integer = DEFAULT_CHUNK_THRESHOLD)
    return Chunk(length(x), threshold)
end

# Constrained to `N <= threshold`, minimize (in order of priority):
#   1. the number of chunks that need to be computed
#   2. the number of "left over" perturbations in the final chunk
function pickchunksize(input_length, threshold = DEFAULT_CHUNK_THRESHOLD)
    if input_length <= threshold
        return input_length
    else
        nchunks = round(Int, input_length / DEFAULT_CHUNK_THRESHOLD, RoundUp)
        return round(Int, input_length / nchunks, RoundUp)
    end
end

chunksize(::Chunk{N}) where {N} = N

########
# Dual #
########

struct Dual{T,V,P} <: Real
    value::V
    partials::P
end

################
# Constructors #
################

# we intend for the right cassette context to
# intercept calls to dualtag
dualtag() = nothing

@inline Dual{T}(value::V, partials::P) where {T,V,P} = Dual{T,V,P}(value, partials)

#@inline Dual{T}(value::V, ::Chunk{N}, p::Val{i}) where {T,V,P,i} = Dual{T}(value, single_seed(Partials{N,V}, p))
@inline Dual(value, partials) = Dual{typeof(dualtag())}(value, partials)

##############################
# Utility/Accessor Functions #
##############################

@inline value(x) = x
@inline value(d::Dual) = d.value

@inline valtype(::V) where {V} = V
@inline valtype(::Type{V}) where {V} = V
@inline valtype(::Dual{T,V}) where {T,V} = V
@inline valtype(::Type{Dual{T,V,P}}) where {T,V,P} = V

@inline partials(d::Dual) = d.partials
@inline partials(d::Number) = Zero()
@inline npartials(d::Dual) = (ps = partials(d)) isa ChainRulesCore.AbstractDifferential ? 1 : length(d.partials)

#####################
# Generic Functions #
#####################

Base.copy(d::Dual) = d

Base.eps(d::Dual) = eps(value(d))
Base.eps(::Type{D}) where {D<:Dual} = eps(valtype(D))

Base.rtoldefault(::Type{D}) where {D<:Dual} = Base.rtoldefault(valtype(D))

Base.floor(::Type{R}, d::Dual) where {R<:Number} = floor(R, value(d))
Base.floor(d::Dual) = floor(value(d))

Base.ceil(::Type{R}, d::Dual) where {R<:Number} = ceil(R, value(d))
Base.ceil(d::Dual) = ceil(value(d))

Base.trunc(::Type{R}, d::Dual) where {R<:Number} = trunc(R, value(d))
Base.trunc(d::Dual) = trunc(value(d))

Base.round(::Type{R}, d::Dual) where {R<:Number} = round(R, value(d))
Base.round(d::Dual) = round(value(d))

Base.hash(d::Dual) = hash(value(d))
Base.hash(d::Dual, hsh::UInt) = hash(value(d), hsh)

function Base.read(io::IO, ::Type{Dual{T,V,P}}) where {T,V,P}
    value = read(io, V)
    partials = read(io, P)
    return Dual{T,V,P}(value, partials)
end

function Base.write(io::IO, d::Dual)
    write(io, value(d))
    write(io, partials(d))
end

@inline Base.zero(d::Dual) = zero(typeof(d))
@inline Base.zero(::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(zero(V), zero(P))

@inline Base.one(d::Dual) = one(typeof(d))
@inline Base.one(::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(one(V), zero(P))

@inline Random.rand(rng::AbstractRNG, d::Dual) = rand(rng, value(d))
@inline Random.rand(::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(rand(V), zero(P))
@inline Random.rand(rng::AbstractRNG, ::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(rand(rng, V), zero(P))
@inline Random.randn(::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(randn(V), zero(P))
@inline Random.randn(rng::AbstractRNG, ::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(randn(rng, V), zero(P))
@inline Random.randexp(::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(randexp(V), zero(P))
@inline Random.randexp(rng::AbstractRNG, ::Type{Dual{T,V,P}}) where {T,V,P} = Dual{T}(randexp(rng, V), zero(P))

# Predicates #
#------------#

isconstant(d::Dual) = iszero(partials(d))

for pred in UNARY_PREDICATES
    @eval Base.$(pred)(d::Dual) = $(pred)(value(d))
end

##############
# Conversion #
##############

Base.convert(::Type{Dual{T,V,P}}, d::Dual{T}) where {T,V,P} = Dual{T}(convert(V, value(d)), convert(P, partials(d)))
Base.convert(::Type{Dual{T,V,P}}, x) where {T,V,P} = Dual{T}(convert(V, x), zero(P))
Base.convert(::Type{Dual{T,V,P}}, x::Number) where {T,V,P} = Dual{T}(convert(V, x), zero(P))
Base.convert(::Type{D}, d::D) where {D<:Dual} = d
Base.float(d::Dual{T,V,P}) where {T,V,P} = convert(Dual{T,float(V),P}, d)
Base.AbstractFloat(d::Dual) = float(d)

###################
# Pretty Printing #
###################

subscript_num(n) = join(reverse(map(d -> '₀' + d, digits(n))))

function tag_show(t, n=0)
    if t isa Nothing
        return subscript_num(n)
    elseif t isa Tag
        tag_show(oldertag(t), n+1)
    else
        return "{" * repr(t) * "}" * subscript_num(n)
    end
end

function Base.show(io::IO, d::Dual{T,V,N}) where {T,V,N}
    print(io, "(")
    print(io, value(d))
    print(io, " + ")
    color = isbits(partials(d)) ? 2 : 3
    if npartials(d) > 0
        printstyled(io, partials(d), color=color)
    else
        printstyled(io, "0", color=color)
    end
    print(io, "ϵ")
    print(io, tag_show(T()))
    print(io, ")")
    return nothing
end


text:
#THIS CODE BASED ON N. REDDY'S THICWEED.JL DISPLAY SCRIPT

#CONSTANTS
# each tuple specifies the x, y, fontsize, xscale to print the character
# in a 100x100 box at position 100,100.
charvals_dict = Dict{Char,Tuple}('A'=> (88.5,199.,135.,1.09),
                 'C'=> (79.,196.,129.,1.19),
                 'G'=> (83.,196.,129.,1.14),
                 'T'=> (80.,199.,135.,1.24))

colour_dict = Dict{Char,String}('A'=>"(0,128,0)", #green
                 'C'=>"(0,0,128)", #blue
                 'G'=>"(255,69,0)", #yellow-brown
                 'T'=>"(150,0,0)")  #red

#char output params for string display of PWMs in nested sampler instrumentation
lwcs_vec=['a','c','g','t']
upcs_vec=['A','C','G','T']
cs_vec=[:green,:blue,:yellow,:red]
thresh_dict=Dict([(0.,'_'),(.25,'.'),(.5,"lwcs"),(1,"upcs")])

source_left_x = 10
xlen = 20
yscale = 250
scale_factor = 0.9
ypixels_per_source = 250
ypad = 10
xpad = 20
xpixels_per_position = 20
fontsize1=60
fontsize2=40

#function to convert ICA_PWM_model sources lists to PWM sequence logo diagrams
function logo_from_model(model::ICA_PWM_model,svg_output::String;freq_sort::Bool=false)
    source_tups = Vector{Tuple{AbstractFloat, Integer, Matrix{AbstractFloat}}}() #(%of sequences w/ source, prior index, weight matrix)
    mix = model.mix_matrix #o x s
    for (prior, source) in enumerate(model.sources)
        push!(source_tups, (sum(model.mix_matrix[:,prior])/size(model.mix_matrix)[1], prior, exp.(source[1])))
    end

    freq_sort && sort(source_tups)

    file = open(svg_output, "w")
    write(file, svg_header(xpad+xpixels_per_position*maximum([length(source[1]) for source in model.sources]),ypixels_per_source*length(model.sources)+ypad))

    curry = 0
    for (frequency, index, source) in source_tups
        curry += yscale
        font1y = curry-190
        ndig = ndigits(index+1)
        write(file, pwm_to_logo(source,source_left_x,curry,xlen,yscale*scale_factor))
        write(file, "<text x=\"10\" y=\"$font1y\" font-family=\"Helvetica\" font-size=\"$fontsize1\" font-weight=\"bold\" >$(index)</text>")
        write(file, "<text x=\"$(20+fontsize1*0.6*ndig)\" y=\"$(font1y-fontsize1+1.5*fontsize2)\" font-family=\"Helvetica\" font-size=\"$fontsize2\" font-weight=\"bold\" >$(frequency*100) % of sequences</text>\n")
    end

    write(file, svg_footer())
    @info "Logo written."
    close(file)
end

function svg_header(canvas_size_x, canvas_size_y)
    return """
    <?xml version="1.0" standalone="no"?>
    <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
    "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
    <svg width="$canvas_size_x" height="$canvas_size_y"
     version="1.1" xmlns="http://www.w3.org/2000/svg">
    """
end

function pwm_to_logo(source,xpos,ypos,xlen,yscale)
    outstr = ""
    for n in 1:size(source,1)
        outstr *= print_weightvec_at_x_y(source[n,:],xpos+xlen*n,ypos,xlen,yscale)
    end
    outstr *= "\n"
    return outstr
end

function pwmstr_to_io(io,source;log=true)
    log && (source=exp.(source))
    for pos in 1:size(source,1)
        char,color=uni_wvec_params(source[pos,:])
        printstyled(io, char; color=color)
    end
end

function uni_wvec_params(wvec) #assign a single unicode char and color symbol for the most informational position in a position weight vector
    wvec.+=10^-99
    wvec = [x/sum(wvec) for x in wvec]
    infscore=(2.0 + sum([x*log(2,x) for x in wvec]))
    infvec = [x*infscore for x in wvec]
    val,idx=findmax(infvec)
    return char,color=get_char_by_thresh(idx,val)
end
                function get_char_by_thresh(idx,val)
                    char = '?'; seen_thresh=0.
                    for (thresh,threshchar) in thresh_dict
                        val >= thresh && thresh >= seen_thresh && (char=threshchar; seen_thresh=thresh)
                    end
                    char=='?' && println(val)
                    char=="lwcs" && (char=lwcs_vec[idx])
                    char=="upcs" && (char=upcs_vec[idx])
                    color=cs_vec[idx]

                    return char,color
                end

function print_weightvec_at_x_y(wvec, xpos, ypos, xlen, yscale)
    # xlen is the length occupied by that column
    # yscale is the total height for 2 information bits i.e. the
    # maximum height available for a "perfect" nucleotide
    outstr = ""
    basestr = "ACGT"
    wvec.+=10^-99  #prevent log(0) = -Inf
    wvec = [x/sum(wvec) for x in wvec] #renorm
    infscore = (2.0 + sum([x*log(2,x) for x in wvec]))
    if infscore==0.0
        return ""
    end
    wvec = [x*infscore*yscale/2.0 for x in wvec]
    # at this point, the sum of all wvec is a maximum of yscale
    wveclist = [(wvec[n],basestr[n]) for n in 1:4]
    wveclist = sort(wveclist)
    curr_ypos = ypos
    for n in 1:4
        curr_ypos -= wveclist[n][1]
        outstr *= print_char_in_rect(wveclist[n][2],xpos,curr_ypos,xlen,wveclist[n][1])
    end
    return outstr
end

blo = 1

function print_char_in_rect(c,x,y,width,height)
    raw_x, raw_y, raw_fontsize, raw_xscale = charvals_dict[c]
    raw_x = (raw_x*raw_xscale-100)/raw_xscale
    raw_y = raw_y-100

    xscale = width/100.0 * raw_xscale
    yscale = height/100.0

    scaled_x = x/xscale + raw_x
    scaled_y = y/yscale + raw_y

    return "<text x=\"$scaled_x\" y=\"$scaled_y\" font-size=\"$raw_fontsize\" font-weight=\"bold\" font-family=\"Helvetica\" fill=\"rgb"*colour_dict[c]*"\" transform=\"scale($xscale,$yscale)\">$c</text>\n"
end

function svg_footer()
    return "</svg>"
end


text:
path = joinpath(dirname(@__FILE__), "..", "..")
push!(Base.LOAD_PATH, path)
datapath = joinpath(dirname(@__FILE__), "data")

using ProxSDP

sympackedlen(n) = div(n*(n+1), 2)
sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    n1, n2 = size(out)
    @assert n == n1 == n2
    c = 0
    for j in 1:n, i in 1:j
        c += 1
        out[i,j] = v[c]
    end
    return out
end
function ivech(v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    out = zeros(n, n)
    ivech!(out, v)

    return out
end

include("moitest.jl")