text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

# Helper function to format seconds part.
function formatsec(sec::Real, prec::Integer, truncate::Bool)
    sec = truncate ? trunc(sec, prec) : round(sec, prec)
    # Seconds of right ascension should be always positive (the hours part
    # holds the sign), so we don't need to take the absolute values.
    sec_frac, sec_int = modf(sec)
    # Format the integer part left padding with zeros.
    sec_int_str = lpad(string(trunc(Integer, sec_int)), 2, "0")
    # Unless precision is 0, format the fractional part with the decimal
    # separator "." followed by seconds rounded to the precision required
    # and right padded with zeros.
    sec_frac_str = prec == 0 ? "" : rpad(string(round(sec_frac, prec+1)),
                                         prec+2, "0")[2:end]
    return sec_string = string(sec_int_str, sec_frac_str)
end

"""
    adstring(ra::Real, dec::Real[, precision::Int=2, truncate::Bool=true]) -> string
    adstring([ra, dec]) -> string
    adstring(dec) -> string
    adstring([ra], [dec]) -> ["string1", "string2", ...]

### Purpose ###

Returns right ascension and declination as string(s) in sexagesimal format.

### Explanation ###

Takes right ascension and declination expressed in decimal format, converts them
to sexagesimal and return a formatted string.  The precision of right ascension
and declination can be specified.

### Arguments ###

Arguments of this function are:

* `ra`: right ascension in decimal degrees.  It is converted to hours before
  printing.
* `dec`: declination in decimal degrees.

The function can be called in different ways:

* Two numeric arguments: first is `ra`, the second is `dec`.
* A 2-tuple `(ra, dec)`.
* One 2-element numeric array: `[ra, dec]`.  A single string is returned.
* One numeric argument: it is assumed only `dec` is provided.
* Two numeric arrays of the same length: `ra` and `dec` arrays.  An array of
  strings is returned.
* An array of 2-tuples `(ra, dec)`.

Optional keywords affecting the output format are always available:

* `precision` (optional integer keyword): specifies the number of digits of
  declination seconds.  The number of digits for righ ascension seconds is
  always assumed to be one more `precision`.  If the function is called with
  only `dec` as input, `precision` default to 1, in any other case defaults to 0.
* `truncate` (optional boolean keyword): if true, then the last displayed digit
  in the output is truncated in precision rather than rounded.  This option is
  useful if `adstring` is used to form an official IAU name (see
  http://vizier.u-strasbg.fr/Dic/iau-spec.htx) with coordinate specification.

### Output ###

The function returns one string if the function was called with scalar `ra` and
`dec` (or only `dec`) or a 2-element array `[ra, dec]`.  If instead it was
feeded with arrays of `ra` and `dec`, an array of strings will be returned.  The
format of strings can be specified with `precision` and `truncate` keywords, see
above.

### Example ###

``` julia
adstring(30.4, -1.23, truncate=true)
# => " 02 01 35.9  -01 13 48"
adstring([30.4, -15.63], [-1.23, 48.41], precision=1)
# => 2-element Array{AbstractString,1}:
#     " 02 01 36.00  -01 13 48.0"
#     "-22 57 28.80  +48 24 36.0"
```
"""
function adstring{T<:AbstractFloat}(ra::T, dec::T;
                  precision::Int=0, truncate::Bool=false)
    # XXX: IDL implementation takes also real values for "precision" and
    # truncates it.  I think it's better to enforce an integer type and cure
    # only negative values.
    precision = precision < 0 ? 0 : precision
    if isnan(ra)
        # If "ra" is NaN, print only declination.
        dec_deg, dec_min, dec_sec = sixty(dec)
        ra_string = ""
    else
        ra_hr, ra_min, ra_sec, dec_deg, dec_min, dec_sec = radec(ra, dec)
        ra_sec_string = formatsec(ra_sec, precision + 1, truncate)
        ra_string = @sprintf("%03.2d %02d %s  ", ra_hr,
                             ra_min, ra_sec_string)
    end
    dec_sec_string = formatsec(dec_sec, precision, truncate)
    dec_string = @sprintf("%+03.2d %02d %s", dec_deg,
                          dec_min, dec_sec_string)
    return string(ra_string, dec_string)
end

adstring(ra::Real, dec::Real;
         precision::Int=0, truncate::Bool=false) =
             adstring(promote(float(ra), float(dec))...,
                      precision=precision, truncate=truncate)

adstring(radec::Tuple{Real, Real};
         precision::Int=0, truncate::Bool=false) =
             adstring(radec..., precision=precision, truncate=truncate)

# When printing only declination, IDL implementation defaults "precision" to 1
# instead of 0.
adstring(dec::Real; precision::Int=1, truncate::Bool=false) =
    adstring(NaN, dec, precision=precision, truncate=truncate)

function adstring{T<:Real}(radec::AbstractArray{T};
                             precision::Int=0, truncate::Bool=false)
    @assert length(radec) == 2 "provide a 2-element [ra, dec] array"
    return adstring(radec[1], radec[2], precision=precision, truncate=truncate)
end

function adstring{R<:Real, D<:Real}(ra::AbstractArray{R},
                                        dec::AbstractArray{D};
                                        precision::Int=0, truncate::Bool=false)
    @assert length(ra) == length(dec) "ra and dec arrays should be of the same length"
    result = similar(ra, AbstractString)
    for i in eachindex(ra)
        result[i] = adstring(ra[i], dec[i],
                             precision=precision, truncate=truncate)
    end
    return result
end

function adstring{N<:Real}(radec::AbstractArray{Tuple{N, N}};
                             precision::Int=0, truncate::Bool=false)
    result = similar(radec, AbstractString)
    for i in eachindex(radec)
        result[i] = adstring(radec[i]...,
                             precision=precision, truncate=truncate)
    end
    return result
end


text:
function parse_kw(ps::ParseState, ::Type{Val{Tokens.ABSTRACT}})
    startbyte = ps.t.startbyte

    # Switch for v0.6 compatability
    if ps.nt.kind == Tokens.TYPE
        # Parsing
        kw1 = INSTANCE(ps)
        format_kw(ps)
        next(ps)
        kw2 = INSTANCE(ps)
        format_kw(ps)

        @catcherror ps startbyte sig = @default ps @closer ps block parse_expression(ps)

        # Construction
        next(ps)
        ret = EXPR{Abstract}(EXPR[kw1, kw2, sig, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")
    else
        # Parsing
        kw = INSTANCE(ps)
        @catcherror ps startbyte sig = @default ps parse_expression(ps)

        # Construction
        ret = EXPR{Abstract}(EXPR[kw, sig], ps.nt.startbyte - startbyte, Variable[], "")
    end
    ret.defs = [Variable(Expr(get_id(sig)), :abstract, ret)]
    return ret
end

function parse_kw(ps::ParseState, ::Type{Val{Tokens.BITSTYPE}})
    startbyte = ps.t.startbyte

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)

    @catcherror ps startbyte arg1 = @default ps @closer ps ws @closer ps wsop parse_expression(ps)
    @catcherror ps startbyte arg2 = @default ps parse_expression(ps)

    # Construction
    ret = EXPR{Bitstype}(EXPR[kw, arg1, arg2], ps.nt.startbyte - startbyte, Variable[], "")
    ret.defs = [Variable(Expr(get_id(arg2)), :bitstype, ret)]

    return ret
end

function parse_kw(ps::ParseState, ::Type{Val{Tokens.PRIMITIVE}})
    startbyte = ps.t.startbyte

    if ps.nt.kind == Tokens.TYPE
        # Parsing
        kw1 = INSTANCE(ps)
        format_kw(ps)
        next(ps)
        kw2 = INSTANCE(ps)
        format_kw(ps)
        @catcherror ps startbyte sig = @default ps @closer ps ws @closer ps wsop parse_expression(ps)
        @catcherror ps startbyte arg = @default ps @closer ps block parse_expression(ps)

        # Construction
        next(ps)
        ret = EXPR{Primitive}(EXPR[kw1, kw2, sig, arg, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")

        ret.defs = [Variable(Expr(get_id(sig)), :bitstype, ret)]
    else
        ret = EXPR{IDENTIFIER}(EXPR[], ps.nt.startbyte - startbyte, Variable[], "primitive")
    end
    return ret
end

function parse_kw(ps::ParseState, ::Type{Val{Tokens.TYPEALIAS}})
    startbyte = ps.t.startbyte

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)

    @catcherror ps startbyte arg1 = @closer ps ws @closer ps wsop parse_expression(ps)
    @catcherror ps startbyte arg2 = parse_expression(ps)

    return EXPR{TypeAlias}(EXPR[kw, arg1, arg2], ps.nt.startbyte - startbyte, Variable[], "")
end

parse_kw(ps::ParseState, ::Type{Val{Tokens.TYPE}}) = parse_struct(ps, TRUE)
parse_kw(ps::ParseState, ::Type{Val{Tokens.IMMUTABLE}}) = parse_struct(ps, FALSE)

# new 0.6 syntax
parse_kw(ps::ParseState, ::Type{Val{Tokens.STRUCT}}) = parse_struct(ps, FALSE)

function parse_kw(ps::ParseState, ::Type{Val{Tokens.MUTABLE}})
    startbyte = ps.t.startbyte

    if ps.nt.kind == Tokens.STRUCT
        kw = INSTANCE(ps)
        format_kw(ps)
        next(ps)
        @catcherror ps startbyte ret = parse_struct(ps, TRUE)
        unshift!(ret.args, kw)
        ret.span += kw.span
    else
        ret = EXPR{IDENTIFIER}(EXPR[], ps.nt.startbyte - ps.t.startbyte, Variable[], "mutable")
    end
    return ret
end

function parse_struct(ps::ParseState, mutable)
    startbyte = ps.t.startbyte
    start_col = ps.t.startpos[2] + 4

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)
    @catcherror ps startbyte sig = @default ps @closer ps block @closer ps ws parse_expression(ps)
    block = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps parse_block(ps, block, start_col)

    # Construction
    T = mutable == TRUE ? Tokens.TYPE : Tokens.IMMUTABLE
    next(ps)
    ret = EXPR{mutable == TRUE ? Mutable : Struct}(EXPR[kw, sig, block, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")
    ret.defs = [Variable(Expr(get_id(sig)), Expr(mutable) ? :mutable : :immutable, ret)]

    return ret
end

ext: jl


text:
aaa=100
println(aaa)
vv=200
println(vv)


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

"""
    struct NestedIterator{T}
        iterators::T # Tuple of functions
        condition::Function
    end

Iterators over the tuples that are produced by a nested for loop.
For instance, if `length(iterators) == 3` , this corresponds to the tuples
`(i1, i2, i3)` produced by:
```
for i1 in iterators[1]()
    for i2 in iterator[2](i1)
        for i3 in iterator[3](i1, i2)
            if condition(i1, i2, i3)
                # produces (i1, i2, i3)
            end
        end
    end
end
```
"""
struct NestedIterator{T,C}
    iterators::T # Tuple of functions
    condition::C
end
function nested(iterators...; condition = (args...) -> true)
    return NestedIterator(iterators, condition)
end
Base.IteratorSize(::Type{<:NestedIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:NestedIterator}) = Base.EltypeUnknown()
function next_iterate(iterators, condition, elems, states, iterator, elem_state)
    while true
        elem_state === nothing && return nothing
        elem, state = elem_state
        elems_states = first_iterate(
            Base.tail(iterators),
            condition,
            (elems..., elem),
            (states..., (iterator, state, elem)),
        )
        elems_states !== nothing && return elems_states
        # This could be written as a recursive function where we call `next_iterate`
        # here with this new value of `next_iterate` instead of the `while` loop`.
        # However, if there are too many consecutive elements for which `condition`
        # is `false` for the last iterator, this will result in a `StackOverflow`.
        # See https://github.com/jump-dev/JuMP.jl/issues/2335
        elem_state = iterate(iterator, state)
    end
end
function first_iterate(::Tuple{}, condition, elems, states)
    if condition(elems...)
        return elems, states
    else
        return nothing
    end
end
function first_iterate(iterators, condition, elems, states)
    iterator = iterators[1](elems...)
    return next_iterate(
        iterators,
        condition,
        elems,
        states,
        iterator,
        iterate(iterator),
    )
end
tail_iterate(::Tuple{}, condition, elems, states, prev_states) = nothing
function tail_iterate(iterators, condition, elems, states, prev_states)
    next = tail_iterate(
        Base.tail(iterators),
        condition,
        (elems..., states[1][3]),
        Base.tail(states),
        (prev_states..., states[1]),
    )
    if next !== nothing
        return next
    end
    iterator = states[1][1]
    return next_iterate(
        iterators,
        condition,
        elems,
        prev_states,
        iterator,
        iterate(iterator, states[1][2]),
    )
end
function Base.iterate(it::NestedIterator)
    return first_iterate(it.iterators, it.condition, tuple(), tuple())
end
function Base.iterate(it::NestedIterator, states)
    return tail_iterate(it.iterators, it.condition, tuple(), states, tuple())
end
function _eltype_or_any(::NestedIterator{<:Tuple{Vararg{Any,N}}}) where {N}
    return NTuple{N,Any}
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# -*- coding: utf-8 -*-
#
# ANISODIFF - Anisotropic diffusion
#
# Usage:
#   diff = anisodiff(img, niter, kappa, lambda, option)
#
# Arguments:
# | img    - input image (2D grayscale)
# | niter  - number of iterations
# | kappa  - conduction coefficient (gradient modulus threshold)
# |          This parameter controls conduction as a function of gradient.
# |          If kappa is low, small intensity gradients are able to block
# |          conduction and hence diffusion across step edges.  A large value
# |          reduces the influence of intensity gradients on conduction.
# | lambda - integration constant for stability (0 <= lambda <= 0.25)
# |          This parameter controls the diffusion speed, and you
# |          usually want it at the maximum value of 0.25.
# |          default value: 0.25
# | option - conduction coefficient functions proposed by Perona & Malik:
# |          1: c(x, y, t) = exp(-(nablaI/kappa).^2)
# |             privileges high-contrast edges over low-contrast ones
# |          2: c(x, y, t) = 1 ./ (1 + (nablaI/kappa).^2)
# |             privileges wide regions over smaller ones
# |          default value: 1
#
# Returns:
# | diff   - anisotropic diffused image
#
# Reference:
# [1] P. Perona and J. Malik.
#     Scale-space and edge detection using ansotropic diffusion.
#     IEEE Transactions on Pattern Analysis and Machine Intelligence,
#     12(7):629-639, July 1990.
#     https://dx.doi.org/10.1109%2F34.56205
#
# Credits:
# [1] Peter Kovesi
#     pk@peterkovesi.com
#     MATLAB and Octave Functions for Computer Vision and Image Processing
#     http://www.peterkovesi.com/matlabfns/Spatial/anisodiff.m
#     --
#     June 2000  original version
#     March 2002 corrected diffusion eqn No 2.
# [2] Daniel Lopes
#     Anisotropic Diffusion (Perona & Malik)
#     http://www.mathworks.com/matlabcentral/fileexchange/14995-anisotropic-diffusion--perona---malik-
#
#
# Aaron LI <aaronly.me@gmail.com>
# 2015/07/17
#

include("calc_k_percentile.jl");

function anisodiff(img, niter, k=calc_k_percentile, lambda=0.25, option=1)
    diff = float(img)
    rows, cols = size(diff)

    for i = 1:niter
        println("anisodiff - iteration: ", i)

        # Construct diffl which is the same as diff but
        # has an extra padding of zeros around it.
        diffl = zeros(rows+2, cols+2)
        diffl[2:rows+1, 2:cols+1] = diff

        # North, South, East and West differences
        deltaN = diffl[1:rows,   2:cols+1] - diff
        deltaS = diffl[3:rows+2, 2:cols+1] - diff
        deltaE = diffl[2:rows+1, 3:cols+2] - diff
        deltaW = diffl[2:rows+1, 1:cols]   - diff

        # Calculate the kappa
        if isa(k, Function)
            kappa = k(diff)
        else
            kappa = k
        end

        println("  kappa: ", kappa)

        # Conduction
        if option == 1
            cN = exp(-(deltaN/kappa).^2)
            cS = exp(-(deltaS/kappa).^2)
            cE = exp(-(deltaE/kappa).^2)
            cW = exp(-(deltaW/kappa).^2)
        elseif option == 2
            cN = 1 ./ (1 + (deltaN/kappa).^2)
            cS = 1 ./ (1 + (deltaS/kappa).^2)
            cE = 1 ./ (1 + (deltaE/kappa).^2)
            cW = 1 ./ (1 + (deltaW/kappa).^2)
        end

        diff += lambda * (cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW)
    end

    return diff
end

max_stars_repo_name: liweitianux/atoolbox

max_issues_repo_name: liweitianux/atoolbox

max_forks_repo_name: liweitianux/atoolbox


text:
using ImplicitGraphs, Permutations

"""
    iTransposition(n::Int, adjacent::Bool=true)

Create an `ImplicitGraph` whose vertices are all permutations of `1:n`
in which two vertices are adjacent if they differ by a transpostion.

* If `adjacent` is true, then only consider transpositions of the form `(a,a+1)`
* Otherwise, consider all possible transpositions of the form `(a,b)` where `1≤a<b≤n`.
"""
function iTransposition(n::Int, adjacent::Bool = true)
    @assert n > 0 "Require n>0, got n=$n"

    vcheck(v::Permutation) = length(v) == n

    function outs(v::Permutation)
        if adjacent
            return [v * Transposition(n, i, i + 1) for i = 1:n-1]
        end
        return [v * Transposition(n, i, j) for i = 1:n for j = 1:n if i ≠ j]
    end

    return ImplicitGraph{Permutation}(vcheck, outs)
end

"""
    trans_string(p::Permutation)

We assume that `p` is a transposition.
"""
function trans_string(p::Permutation)::String
    n = length(p)
    FP = fixed_points(p)
    @assert n == length(FP) + 2 "Permutation must be a transposition"
    elts = sort(setdiff(1:n, FP))
    a, b = elts
    return "($a,$b)"
end

"""
    pscore(p::Permutation)::Int

Measures the extent to which `p` is different from the identity permutation
using this formula:

`sum(abs(p[k]-k) for k=1:length(p))`
"""
function pscore(p::Permutation)::Int
    n = length(p)
    sum(abs(i - p(i)) for i = 1:n)
end

"""
    crazy_trans_product(p::Permutation, adjacent::Bool=true)::String

Express `p` as the product of transpositions. With `adjacent=true` the transpositions
are all of the form `(a,a+1)`; otherwise, all possible transpositions are allowed.

Result is expressed as a `String`.
"""
function crazy_trans_product(p::Permutation, adjacent::Bool = true)::String
    n = length(p)
    G = iTransposition(n, adjacent)
    P = reverse(guided_path_finder(G, p, Permutation(n), score = pscore))

    nP = length(P)

    if nP == 1
        return "()"
    end

    Q = [P[k]' * P[k+1] for k = 1:nP-1]

    prod(trans_string(Q[j]) for j = 1:length(Q))
end


text:
@recipe function f(sol::AbstractNLOSolution;plot_analytic=false, evolution=false)
    ω = sol.prob.ω
    zend = sol.prob.prob.tspan[2]
    spectrum = abs2.(FT(sol, zend))
    layout := (2, 1)
    if evolution
        zmesh = linspace(0, zend, 100)
        intensity = 10*log10.(1e-30+abs2.(sol.(zmesh)))
        intensity = hcat(intensity...)
        spectrum = 10*log10.(1e-30+abs2.(FT.(sol, zmesh)))
        spectrum = hcat(spectrum...)
        @series begin
            title := "Intensity"
            subplot := 1
            #ylims := (minimum(sol.prob.tmesh)/25, maximum(sol.prob.tmesh)/25)
            clims := (maximum(intensity)-80, maximum(intensity))
            zmesh, sol.prob.tmesh, intensity
        end
        @series begin
            title := "Spectrum"
            subplot := 2
            clims := (maximum(spectrum)-80, maximum(spectrum))
            zmesh, ω, spectrum
        end
    else
        @series begin
            title := "Intensity"
            subplot := 1
            seriestype := :path
            sol.prob.tmesh, abs2.(sol(zend))
        end
        @series begin
            ylabel := ("(dB)")
            title := "Spectrum"
            subplot := 2
            seriestype := :path
            ω, spectrum/maximum(spectrum)
        end
    end
end

@recipe function f(optical_property::OpticalAttr)
    fs = frequency(optical_property)/1e12 #plot in THz
    λs = wavelength(optical_property)*1e9 #plot in nm
    @series begin
        xlabel := "Frequency (THz)"
        ylabel := get_label(optical_property)
        seriestype := :path
        xmirror := true
        fs, get_property(optical_property)
    end
end


text:
import Base: length
length(p::BifurcationPath) = length(p.p_vals)
length(p::BifurcationGridDiagram) = p.length

@recipe function f(path::BifurcationPath, var=1)
    linewidth --> 4
    color --> stab_color.(path.stability_types)
    path.p_vals, getindex.(path.vals, var)
end

@recipe function f(path::BifurcationPath, var, z_val)
    linewidth --> 4
    color --> stab_color.(path.stability_types)
    fill(z_val, length(path)), path.p_vals, getindex.(path.vals, var)
end

@recipe function f(bif::BifurcationDiagram, var, z_val)
    xlabel --> string(bif.param)
        label --> ""
    for (i, path) in enumerate(bif.paths)
        @series begin
            primary --> i==1
            path, var, z_val
        end
    end
end

@recipe function f(bif::BifurcationDiagram, var=1)
    xlabel --> string(bif.param)
        label --> ""
    for (i, path) in enumerate(bif.paths)
        @series begin
            primary --> i==1
            path, var
        end
    end
end

@recipe function f(grid::BifurcationGridDiagram, var=1)
    label --> ""
    zlabel --> "Concentration"
    ylabel --> string(grid.param2)
    xlabel --> string(grid.param1)
    for (i, bif) in enumerate(grid.bifurcation_diagrams)
        @series begin
            bif, var, grid.range1[i]
        end
    end
end

@recipe function f(u, point::BifurcationPoint, var)
    seriestype := :scatter
    color --> stab_color.(point.stability_types)
    fill(u, length(point.vals)), getindex.(point.vals,var)
end

@recipe function f(u, y, point::BifurcationPoint, var)
    seriestype := :scatter
    color --> stab_color.(point.stability_types)
    fill(u, length(point.vals)), fill(y, length(point.vals)), getindex.(point.vals,var)
end

@recipe function f(grid::BifurcationGrid, var=1)
    label --> ""
    ylabel --> "Concentration"
    xlabel --> string(grid.param)
    for (i, point) in enumerate(grid.grid_points)
        @series begin
            primary --> i==1
            grid.range[i], point, var
        end
    end
end

@recipe function f(grid::BifurcationGrid2D, var=1)
    label --> ""
    zlabel --> "Concentration"
    ylabel --> string(grid.param2)
    xlabel --> string(grid.param1)
    for i = 1:length(grid.range1), j = 1:length(grid.range2)
        @series begin
            primary --> i==1
            grid.range1[i], grid.range2[j], grid.grid_points[i,j], var
        end
    end
end

#For a specific stability types, returns a color (to be used to distinguish the types in plots).
function stab_color(stab_type::Int8)
    (stab_type==0) && (return :red)
    (stab_type==1) && (return :blue)
    (stab_type==2) && (return :orange)
    (stab_type==3) && (return :cyan)
end


text:
# create files in a temporary dir that we can modify
const tmpdir = mktempdir()
const file1 = joinpath(tmpdir, "file1")
const file2 = joinpath(tmpdir, "file2")
write(file1, ".")
write(file2, ".")

@testset "Watcher/WatchedFile struct  " begin
    wf1 = LS.WatchedFile(file1)
    wf2 = LS.WatchedFile(file2)

    # Basic struct
    @test wf1.path == file1
    @test wf2.path == file2
    @test wf1.mtime == mtime(file1)
    @test wf2.mtime == mtime(file2)

    # Apply change and check if it's detecte
    t1 = time()
    sleep(FS_WAIT)
    write(file1, "hello")
    sleep(FS_WAIT)
    @test LS.has_changed(wf1) == 1
    @test LS.has_changed(wf2) == 0

    # Set state as unchanged
    LS.set_unchanged!(wf1)
    @test LS.has_changed(wf1) == 0
    @test wf1.mtime > t1
end

@testset "Watcher/SimpleWatcher struct" begin
    sw  = LS.SimpleWatcher()

    isa(sw, LS.FileWatcher)

    sw1 = LS.SimpleWatcher(identity, sleeptime=0.0)

    # Base constructor check
    @test sw.callback === nothing
    @test sw.task === nothing
    @test sw.sleeptime == 0.1
    @test isempty(sw.watchedfiles)
    @test eltype(sw.watchedfiles) == LS.WatchedFile

    @test sw1.sleeptime == 0.05 # via the clamping
    @test sw1.callback(2) == 2 # identity function
    @test sw1.callback("blah") == "blah"
    @test isempty(sw1.watchedfiles)
    @test sw1.task === nothing
end

@testset "Watcher/watch  file routines" begin
    sw = LS.SimpleWatcher(identity)

    LS.watch_file!(sw, file1)
    LS.watch_file!(sw, file2)

    @test sw.watchedfiles[1].path == file1
    @test sw.watchedfiles[2].path == file2

    # is_watched
    @test LS.is_watched(sw, file1)
    @test LS.is_watched(sw, file2)

    # is_running?
    @test !LS.is_running(sw)

    LS.start(sw)
    sleep(0.001)

    @test LS.is_running(sw)
    @test LS.stop(sw)
    @test !LS.is_running(sw)

    #
    # modify callback to something that will eventually throw an error
    #
    LS.set_callback!(sw, log)
    @test sw.callback(exp(1.0)) ≈ 1.0

    LS.start(sw)
    sleep(0.001)

    # causing a modification will generate an error because the callback
    # function will fail on a string
    cray = Crayon(foreground=:cyan, bold=true)
    println(cray, "\n⚠ Deliberately causing an error to be displayed and handled...\n")
    write(file1, "modif")
    sleep(0.25) # needs to be sufficient to give time for propagation.
    @test sw.status == :interrupted

    #
    # deleting files
    #

    file3 = joinpath(tmpdir, "file3")
    write(file3, "hello")

    sw = LS.SimpleWatcher(identity)

    LS.watch_file!(sw, file1)
    LS.watch_file!(sw, file2)
    LS.watch_file!(sw, file3)

    @test length(sw.watchedfiles) == 3

    LS.start(sw)

    rm(file3)
    sleep(0.25) # needs to be sufficient to give time for propagation.

    # file3 was deleted
    @test length(sw.watchedfiles) == 2
    @test sw.watchedfiles[1].path == file1
    @test sw.watchedfiles[2].path == file2
end


text:
function Base.getproperty(m::KeyPoint, s::Symbol)
    if s==:pt
        return cpp_to_julia(jlopencv_KeyPoint_get_pt(m))
    end
    if s==:size
        return cpp_to_julia(jlopencv_KeyPoint_get_size(m))
    end
    if s==:angle
        return cpp_to_julia(jlopencv_KeyPoint_get_angle(m))
    end
    if s==:response
        return cpp_to_julia(jlopencv_KeyPoint_get_response(m))
    end
    if s==:octave
        return cpp_to_julia(jlopencv_KeyPoint_get_octave(m))
    end
    if s==:class_id
        return cpp_to_julia(jlopencv_KeyPoint_get_class_id(m))
    end
    return Base.getfield(m, s)
end
function Base.setproperty!(m::KeyPoint, s::Symbol, v)
    return Base.setfield!(m, s, v)
end

function KeyPoint(x::Float32, y::Float32, _size::Float32, _angle::Float32, _response::Float32, _octave::Int32, _class_id::Int32)
	return cpp_to_julia(jlopencv_cv_cv_KeyPoint_cv_KeyPoint_KeyPoint(julia_to_cpp(x),julia_to_cpp(y),julia_to_cpp(_size),julia_to_cpp(_angle),julia_to_cpp(_response),julia_to_cpp(_octave),julia_to_cpp(_class_id)))
end
KeyPoint(x::Float32, y::Float32, _size::Float32; _angle::Float32 = Float32(-1), _response::Float32 = Float32(0), _octave::Int32 = Int32(0), _class_id::Int32 = Int32(-1)) = KeyPoint(x, y, _size, _angle, _response, _octave, _class_id)

function VideoCapture(filename::String, apiPreference::Int32)
	return cpp_to_julia(jlopencv_cv_cv_VideoCapture_cv_VideoCapture_VideoCapture(julia_to_cpp(filename),julia_to_cpp(apiPreference)))
end
VideoCapture(filename::String; apiPreference::Int32 = Int32(CAP_ANY)) = VideoCapture(filename, apiPreference)

function VideoCapture(index::Int32, apiPreference::Int32)
	return cpp_to_julia(jlopencv_cv_cv_VideoCapture_cv_VideoCapture_VideoCapture(julia_to_cpp(index),julia_to_cpp(apiPreference)))
end
VideoCapture(index::Int32; apiPreference::Int32 = Int32(CAP_ANY)) = VideoCapture(index, apiPreference)

function CascadeClassifier(filename::String)
	return cpp_to_julia(jlopencv_cv_cv_CascadeClassifier_cv_CascadeClassifier_CascadeClassifier(julia_to_cpp(filename)))
end

function detect(cobj::cv_Ptr{T}, image::InputArray, mask::InputArray) where {T <: Feature2D}
	return cpp_to_julia(jlopencv_cv_cv_Feature2D_cv_Feature2D_detect(julia_to_cpp(cobj),julia_to_cpp(image),julia_to_cpp(mask)))
end
detect(cobj::cv_Ptr{T}, image::InputArray; mask::InputArray = (CxxMat())) where {T <: Feature2D} = detect(cobj, image, mask)

function detectMultiScale(cobj::CascadeClassifier, image::InputArray, scaleFactor::Float64, minNeighbors::Int32, flags::Int32, minSize::Size{Int32}, maxSize::Size{Int32})
	return cpp_to_julia(jlopencv_cv_cv_CascadeClassifier_cv_CascadeClassifier_detectMultiScale(julia_to_cpp(cobj),julia_to_cpp(image),julia_to_cpp(scaleFactor),julia_to_cpp(minNeighbors),julia_to_cpp(flags),julia_to_cpp(minSize),julia_to_cpp(maxSize)))
end
detectMultiScale(cobj::CascadeClassifier, image::InputArray; scaleFactor::Float64 = Float64(1.1), minNeighbors::Int32 = Int32(3), flags::Int32 = Int32(0), minSize::Size{Int32} = (Size{Int32}(0,0)), maxSize::Size{Int32} = (Size{Int32}(0,0))) = detectMultiScale(cobj, image, scaleFactor, minNeighbors, flags, minSize, maxSize)

function empty(cobj::CascadeClassifier)
	return cpp_to_julia(jlopencv_cv_cv_CascadeClassifier_cv_CascadeClassifier_empty(julia_to_cpp(cobj)))
end

function read(cobj::VideoCapture, image::InputArray)
	return cpp_to_julia(jlopencv_cv_cv_VideoCapture_cv_VideoCapture_read(julia_to_cpp(cobj),julia_to_cpp(image)))
end
read(cobj::VideoCapture; image::InputArray = (CxxMat())) = read(cobj, image)

function release(cobj::VideoCapture)
	return cpp_to_julia(jlopencv_cv_cv_VideoCapture_cv_VideoCapture_release(julia_to_cpp(cobj)))
end

function SimpleBlobDetector_create(parameters::SimpleBlobDetector_Params)
	return cpp_to_julia(jlopencv_cv_cv_SimpleBlobDetector_create(julia_to_cpp(parameters)))
end
SimpleBlobDetector_create(; parameters::SimpleBlobDetector_Params = (SimpleBlobDetector_Params())) = SimpleBlobDetector_create(parameters)

function imread(filename::String, flags::Int32)
	return cpp_to_julia(jlopencv_cv_cv_imread(julia_to_cpp(filename),julia_to_cpp(flags)))
end
imread(filename::String; flags::Int32 = Int32(IMREAD_COLOR)) = imread(filename, flags)

function imshow(winname::String, mat::InputArray)
	return cpp_to_julia(jlopencv_cv_cv_imshow(julia_to_cpp(winname),julia_to_cpp(mat)))
end

function namedWindow(winname::String, flags::Int32)
	return cpp_to_julia(jlopencv_cv_cv_namedWindow(julia_to_cpp(winname),julia_to_cpp(flags)))
end
namedWindow(winname::String; flags::Int32 = Int32(WINDOW_AUTOSIZE)) = namedWindow(winname, flags)

function waitKey(delay::Int32)
	return cpp_to_julia(jlopencv_cv_cv_waitKey(julia_to_cpp(delay)))
end
waitKey(; delay::Int32 = Int32(0)) = waitKey(delay)

function rectangle(img::InputArray, pt1::Point{Int32}, pt2::Point{Int32}, color::Scalar, thickness::Int32, lineType::Int32, shift::Int32)
	return cpp_to_julia(jlopencv_cv_cv_rectangle(julia_to_cpp(img),julia_to_cpp(pt1),julia_to_cpp(pt2),julia_to_cpp(color),julia_to_cpp(thickness),julia_to_cpp(lineType),julia_to_cpp(shift)))
end
rectangle(img::InputArray, pt1::Point{Int32}, pt2::Point{Int32}, color::Scalar; thickness::Int32 = Int32(1), lineType::Int32 = Int32(LINE_8), shift::Int32 = Int32(0)) = rectangle(img, pt1, pt2, color, thickness, lineType, shift)

function cvtColor(src::InputArray, code::Int32, dst::InputArray, dstCn::Int32)
	return cpp_to_julia(jlopencv_cv_cv_cvtColor(julia_to_cpp(src),julia_to_cpp(code),julia_to_cpp(dst),julia_to_cpp(dstCn)))
end
cvtColor(src::InputArray, code::Int32; dst::InputArray = (CxxMat()), dstCn::Int32 = Int32(0)) = cvtColor(src, code, dst, dstCn)

function equalizeHist(src::InputArray, dst::InputArray)
	return cpp_to_julia(jlopencv_cv_cv_equalizeHist(julia_to_cpp(src),julia_to_cpp(dst)))
end
equalizeHist(src::InputArray; dst::InputArray = (CxxMat())) = equalizeHist(src, dst)

function destroyAllWindows()
	return cpp_to_julia(jlopencv_cv_cv_destroyAllWindows())
end

function getTextSize(text::String, fontFace::Int32, fontScale::Float64, thickness::Int32)
	return cpp_to_julia(jlopencv_cv_cv_getTextSize(julia_to_cpp(text),julia_to_cpp(fontFace),julia_to_cpp(fontScale),julia_to_cpp(thickness)))
end

function putText(img::InputArray, text::String, org::Point{Int32}, fontFace::Int32, fontScale::Float64, color::Scalar, thickness::Int32, lineType::Int32, bottomLeftOrigin::Bool)
	return cpp_to_julia(jlopencv_cv_cv_putText(julia_to_cpp(img),julia_to_cpp(text),julia_to_cpp(org),julia_to_cpp(fontFace),julia_to_cpp(fontScale),julia_to_cpp(color),julia_to_cpp(thickness),julia_to_cpp(lineType),julia_to_cpp(bottomLeftOrigin)))
end
putText(img::InputArray, text::String, org::Point{Int32}, fontFace::Int32, fontScale::Float64, color::Scalar; thickness::Int32 = Int32(1), lineType::Int32 = Int32(LINE_8), bottomLeftOrigin::Bool = (false)) = putText(img, text, org, fontFace, fontScale, color, thickness, lineType, bottomLeftOrigin)

include("cv_dnn_wrap.jl")


text:
using Intervals
using IterTools: ncycle

global DEBUG = false

""" Types """

mutable struct FormulaStats
    bools # boolean vars. arraylike
    reals # real vars. arraylike
    bool_macros #arraylike
    new_bool_var_count::Int
    bool_macro_count::Int
end
FormulaStats() = FormulaStats([],[],[],0,0)

mutable struct SMTLibFormula
    formula # arraylike
    stats::FormulaStats
end
SMTLibFormula() = SMTLibFormula([], FormulaStats())

mutable struct MyError
    message::String
end

"""Printing/Writing Functions"""
function Base.show(io::IO, f::SMTLibFormula)
    s = "SMTLibFormula: "
    s *= "new_bool_var_count = " * string(f.stats.new_bool_var_count)
    s *= ", bools: " * string(f.stats.bools)
    s *= ", reals: " * string(f.stats.reals)
    println(io, s)
end

function write_to_file(f::SMTLibFormula, fname; dirname="sound/smtlibfiles/")
    # print expressions to file
    @debug println(join(f.formula, "\n"))

    # TODO: make dir before writing to file in it
    try
        mkdir(pwd() * "/" * dirname) # make dir if it doesn't exist
    catch
    end
    full_fname = pwd()* "/" * dirname * fname
    file = open(full_fname, "w")
    write(file, join(f.formula, "\n"))
    close(file)
    return full_fname
end

function write_result(fname, result; specifier="w")
    # results will be put in a txt  file of the same name but with "result" appended
    io = open(fname[1:end-5]*"_result.txt", specifier)
    write(io, result...)
    close(io)
    return nothing
end

function read_result(fname)
    # results will be put in a txt  file of the same name but with "result" appended
    io = open(fname[1:end-5]*"_result.txt", "r")
    result = read(io, String)
    close(io)
    return result
end

function define_domain(d, stats::FormulaStats)
    # d is a dictionary denoting the domain, e.g. {"x" => [-.3, 5.6], ...}
    assertions = []
    for (k,v) in d
        lb = v[1]
        ub = v[2]
        if k ∉ stats.reals
            add_real_var(k, stats)
        end
        box = assert_statement(define_box(string(k),lb, ub))
        push!(assertions, box)
    end
    return assertions
end

function handle_unary_negation(n::T where T <: Real)
    return n < 0 ? prefix_notate("-", [-n]) : string(n)
end
function handle_unary_negation(n::Symbol)
    return string(n)
end

function define_box(v::String, lb, ub)
    lb = handle_unary_negation(lb)
    ub = handle_unary_negation(ub)
    lb_e = prefix_notate("<=", [v, ub])
    ub_e = prefix_notate(">=", [v, lb])
    return prefix_notate("and", [lb_e, ub_e])
end

"""Low level functions for converting ϕ and ϕ̂ to smtlib2"""

"""
f is an array representing a conjunction.
Returns an array
"""
function assert_all(f::Array, fs::FormulaStats)
    assertions = []
    for item in f
        expr = convert_any_constraint(item, fs)::String
        push!(assertions, assert_statement(expr))
    end
    return assertions
end

function assert_conjunction(f::Array, fs::FormulaStats; conjunct_name=nothing)
    if length(f) == 1
        return [assert_literal(f[1], fs)]
    elseif length(f) > 1
        # assert conjunction
        return assert_actual_conjunction(f, fs; conjunct_name=conjunct_name)::Array
    else # empty list
        return []
    end
end

function assert_literal(l, fs::FormulaStats)
    return assert_statement(convert_any_constraint(l, fs::FormulaStats))
end

function assert_negated_literal(l, fs::FormulaStats)
    return assert_statement(negate(convert_any_constraint(l, fs::FormulaStats)))
end

function assert_negation_of_conjunction(f::Array, fs::FormulaStats; conjunct_name=nothing)
    if length(f) == 1
        return [assert_negated_literal(f[1], fs)]
    elseif length(f) >= 1
        return assert_actual_negation_of_conjunction(f, fs; conjunct_name=conjunct_name)::Array
    else # empty list
        return []
    end
end

function not(atom)
    return ~atom::Bool
end

function add_real_var(v, fs::FormulaStats)
    if not(v in fs.reals)
        push!(fs.reals, v)
    end
end

function add_real_vars(vlist::Array, fs::FormulaStats)
    for v in vlist
        add_real_var(v, fs)
    end
end

function get_new_bool(fs::FormulaStats)
    fs.new_bool_var_count += 1
    v = "b"*string(fs.new_bool_var_count) # b for boolean
    @assert not(v in fs.bools)
    push!(fs.bools, v)
    return v
end

"""
Creates prefix notation syntax of smtlib2.
Turn an op into its prefix form: (op arg1 arg2 ...).
A 'pure' style function that doesn't modify state.
"""
function prefix_notate(op, args)
    expr = "(" * op * " "
    expr *= print_args(args)
    expr *= ")"
    return expr
end

function print_args(args::Array)
    s = ""
    for a in args
        s *= string(a) * " "
    end
    s = s[1:end-1] # chop trailing whitespace
    return s
end

function declare_const(constname, consttype)
    return prefix_notate("declare-const", [constname, consttype])
end

# e.g. (define-fun _def1 () Bool (<= x 2.0))
function define_fun(name, args, return_type, body)
    arg_string = "("*print_args(args)*")"
    return prefix_notate("define-fun", [name, arg_string, return_type, body])
end

function declare_reals(fs::FormulaStats)
    real_decls = []
    for v in fs.reals
        push!(real_decls, declare_const(v, "Real"))
    end
    return real_decls
end

function define_atom(atomname, atomvalue)
    eq_expr = prefix_notate("=", [atomname, atomvalue])
    return assert_statement(eq_expr)
end

function define_bool_macro(macro_name, expr)
    return define_macro(macro_name, expr; return_type="Bool")
end

# (define-fun _def1 () Bool (<= x 2.0))
function define_macro(macro_name, expr; return_type="Bool")
    return define_fun(macro_name, [], return_type, expr)
end

function assert_statement(expr)
    return prefix_notate("assert", [expr])
end

function negate(expr)
    return prefix_notate("not", [expr])
end

function footer()
    return ["(check-sat)", "(get-model)"]
end

function header()
    h = [produce_models(), set_logic()]
    # push!(h, define_max(), define_min(), define_relu()) # min and max already defined by dreal
    return h
end

function set_logic()
    return "(set-logic QF_NRA)"
end

function produce_models()
    return "(set-option :produce-models true)"
end

# define-fun is a macro.
# if x < y return y else x
function define_max()
    return "(define-fun max ((x Real) (y Real)) Real (ite (< x y) y x))"
end

# if x < y return x else y
function define_min()
    return "(define-fun min ((x Real) (y Real)) Real (ite (< x y) x y))"
end

function define_relu()
    return "(define-fun relu ((x Real)) Real (max x 0))"
end

function assert_actual_conjunction(constraint_list, fs::FormulaStats; conjunct_name=nothing)
    formula, conjunct_name = declare_conjunction(constraint_list, fs; conjunct_name=conjunct_name) # declare conjunction
    push!(formula, assert_statement(conjunct_name)) # assert conjunction
    return formula
end

function assert_actual_negation_of_conjunction(constraint_list, fs::FormulaStats; conjunct_name=nothing)
    """
    Assert the negation of conjunction of the constraints passed in constraint_list.
    not (A and B and C and ...)
    """
    formula, conjunct_name = declare_conjunction(constraint_list, fs; conjunct_name=conjunct_name) # declare conjunction
    push!(formula, assert_statement(negate(conjunct_name))) # assert NEGATED conjunction
    return formula
end

function declare_conjunction(constraint_list, fs::FormulaStats; conjunct_name=nothing, use_macros=false)
    """
    Given a list of constraints, declare their conjunction but DO NOT
    assert their conjunction.
    can use macros e.g.
    (define-fun _def1 () Bool (<= x 2.0))
    (define-fun _def2 () Bool (>= y 3.0))

    or can use new boolean variables
    (declare-const b1 Bool)
    (assert (== b1 (<= x 2.0)))

    ...
    (declare ... phi)
    (assert (= phi (and A B)))

    But notice we are just _defining_ phi, we are not asserting that
    phi _holds_, which would be: (assert phi) [not doing that tho!]
    """
    defs, names = declare_list(constraint_list, fs; use_macros=use_macros)
    if isnothing(conjunct_name)
        conjunct_name = get_new_bool(fs)
    end
    @assert length(names) > 1
    conjunct = prefix_notate("and", names)
    conjunct_decl = [declare_const(conjunct_name, "Bool")]
    conjunct_def = [define_atom(conjunct_name, conjunct)]
    formula = vcat(defs, conjunct_decl, conjunct_def)
    return formula, conjunct_name
end

function assert_disjunction(constraint_list, fs::FormulaStats; disjunct_name=nothing)
    throw(MyError("NotImplementedError"))
end

"""
Convert julian relational operators into smtlib2 friendly ones.
"""
function convert_f(f::Symbol)
    if f == :(==)
        return "="
    elseif f == :≤ || f == :≦
        return "<="
    elseif f == :≥ || f == :≧
        return ">="
    else
        return string(f)
    end
end

function convert_any_constraint(c::Expr, fs::FormulaStats)
    # basically just prefix notate the constraint and take from expr -> string
    # base case: numerical number
    try
        constant = eval(c)
        return convert_any_constraint(constant::Real, fs)
    catch e
    end
    # recursive case
    f = convert_f(c.args[1])
    args = c.args[2:end]
    converted_args = []
    for a in args
        push!(converted_args, convert_any_constraint(a, fs))
    end
    return prefix_notate(string(f), converted_args)
end
# base cases:
function convert_any_constraint(s::Symbol, fs::FormulaStats)
    add_real_var(s, fs) # log var
    return string(s)
end
function convert_any_constraint(n::Real, fs::FormulaStats)
    if n >= 0
        return string(n)
    else #if n < 0
        return prefix_notate("-", [-n]) # smtlib2 can't handle negative numbers and needs unary minus, ugh (!)
    end
end

function declare_list(constraint_list::Array, fs::FormulaStats; use_macros=false)
    """
    turn a list of some type of AbstractConstraint <: into

    smtlib macro declarations/definitions:
    (define-fun _def1 () Bool (<= x 2.0))

    But DON'T assert either true or false for the macro (e.g. (assert _def1()) )

    OR new binary variables:

    (declare-const b1 Bool)
    (assert (= b1 (<= v 9)))

    And again DON'T assert either (assert b1) or (assert (not b1))
    """
    defs = [] # definitions + declarations
    names = [] # names
    for item in constraint_list
        expr = convert_any_constraint(item, fs)::String
        if use_macros
            macro_name = get_new_macro(fs)
            push!(names, macro_name)
            push!(defs, define_bool_macro(macro_name, expr))
        else # use new bools
            bool_name = get_new_bool(fs)
            push!(names, bool_name)
            push!(defs, declare_const(bool_name, "Bool"))
            push!(defs, define_atom(bool_name, expr))
        end
    end
    return defs, names
end

function get_new_macro(fs::FormulaStats)
    """ Get new macro name """
    fs.bool_macro_count += 1
    v = "_def"*string(fs.bool_macro_count)
    @assert not(v in fs.bool_macros)
    push!(fs.bool_macros, v)
    return v
end


text:
@testset "SGS" begin
  𝒮 = georef((z=[1.,0.,1.],), [25. 50. 75.; 25. 75. 50.])
  𝒟 = CartesianGrid((100,100), (0.5,0.5), (1.0,1.0))
  N = 3

  problem₁ = SimulationProblem(𝒮, 𝒟, :z, N)
  problem₂ = SimulationProblem(𝒟, :z=>Float64, N)

  solver = SGS(:z => (variogram=SphericalVariogram(range=35.),
                      neighborhood=NormBall(30.)))

  Random.seed!(2017)
  sol₁ = solve(problem₁, solver)
  sol₂ = solve(problem₂, solver)

  # basic checks
  reals = sol₁[:z]
  inds = LinearIndices(size(𝒟))
  @test all(reals[i][inds[25,25]] == 1. for i in 1:N)
  @test all(reals[i][inds[50,75]] == 0. for i in 1:N)
  @test all(reals[i][inds[75,50]] == 1. for i in 1:N)
end


text:
export reproject, convertmode, mappedindex, mappedbounds

# These methods are only available if ArchGDAL is loaded.
# Otherwise Projected selector crs field must be `nothing`,
# and no reprojection can occur.

DD._sel2indices(mode::Projected, dim::Dimension, sel::Contains) = begin
    selval = reproject(mappedcrs(mode), crs(mode), dim, val(sel))
    DD.contains(dim, rebuild(sel, selval))
end
DD._sel2indices(mode::Projected, dim::Dimension, sel::At) = begin
    selval = reproject(mappedcrs(mode), crs(mode), dim, val(sel))
    DD.at(dim, rebuild(sel, selval))
end
DD._sel2indices(mode::Projected, dim::Dimension, sel::Between) = begin
    selval = map(v -> reproject(mappedcrs(mode), crs(mode), dim, v), val(sel))
    DD.between(dim, rebuild(sel, selval))
end

"""
`reproject` uses ArchGDAL.reproject, but implemented for a reprojecting
a single dimension at a time.
"""
reproject(source, target, dim::Dimension, val) = val
reproject(source::GeoFormat, target::GeoFormat, dim::Lon, val::Number) =
    ArchGDAL.reproject((Float64(val), 0.0), source, target; order=:trad)[1]
reproject(source::GeoFormat, target::GeoFormat, dim::Lat, val::Number) =
    ArchGDAL.reproject((0.0, Float64(val)), source, target; order=:trad)[2]

reproject(source::GeoFormat, target::GeoFormat, ::Lon, vals::AbstractArray) =
    [r[1] for r in ArchGDAL.reproject([(Float64(v), 0.0) for v in vals], source, target; order=:trad)]
reproject(source::GeoFormat, target::GeoFormat, dim::Lat, vals::AbstractArray) =
    [r[2] for r in ArchGDAL.reproject([(0.0, Float64(v)) for v in vals], source, target; order=:trad)]

reproject(source::GeoFormat, target::GeoFormat, ::Lon, vals::Tuple) =
    Tuple(r[1] for r in ArchGDAL.reproject([(Float64(v), 0.0) for v in vals], source, target; order=:trad))
reproject(source::GeoFormat, target::GeoFormat, dim::Lat, vals::Tuple) =
    Tuple(r[2] for r in ArchGDAL.reproject([(0.0, Float64(v)) for v in vals], source, target; order=:trad))

convertmode(dstmode::Type{Mapped}, srcmode::Type{Projected}, dim::Dimension) where M = begin
    m = mode(dim)
    newindex = reproject(crs(m), mappedcrs(m), dim, val(dim))
    newbounds = reproject(crs(m), mappedcrs(m), dim, bounds(dim))
    newmode = Mapped(
        order=order(m),
        span=Irregular(newbounds),
        sampling=sampling(m),
        crs=crs(m),
        mappedcrs=mappedcrs(m),
    )
    rebuild(dim; val=newindex, mode=newmode)
end
convertmode(dstmode::Type{Projected}, srcmode::Type{Mapped}, dim::Dimension) where M = begin
    m = mode(dim)
    newindex = _projectedrange(m, dim)
    newmode = Projected(
        order=order(m),
        span=Regular(step(newindex)), sampling=sampling(m),
        crs=crs(m),
        mappedcrs=mappedcrs(m),
    )
    rebuild(dim; val=newindex, mode=newmode)
end

_projectedrange(::Projected, dim) = LinRange(first(dim), last(dim), length(dim))
_projectedrange(m::Mapped, dim) =
    _projectedrange(span(m), crs(m), m, dim)
_projectedrange(span, crs, m::Mapped, dim) = begin
    start, stop = reproject(mappedcrs(m), crs, dim, [first(dim), last(dim)])
    LinRange(start, stop, length(dim))
end
_projectedrange(::Regular, crs::Nothing, ::Mapped, dim) =
    LinRange(first(dim), last(dim), length(dim))
_projectedrange(::Irregular, crs::Nothing, ::Mapped, dim) =
    error("Cannot convert an Mapped Irregular index to Projected when projectioncrs is nothing")

projectedbounds(mode::Mapped, dim) = projectedbounds(crs(mode), mode, dim)
projectedbounds(crs::Nothing, mode::Mapped, dim) =
    error("No projection crs attached to $(name(dim)) dimension")
projectedbounds(crs::GeoFormat, mode::Mapped, dim) =
    reproject(mappedcrs(mode), crs, dim, bounds(dim))

projectedindex(mode::Mapped, dim) = projectedindex(crs(mode), mode, dim)
projectedindex(crs::Nothing, mode::Mapped, dim) =
    error("No projection crs attached to $(name(dim)) dimension")
projectedindex(crs::GeoFormat, mode::Mapped, dim) =
    reproject(mappedcrs(dim), crs, dim, index(dim))

mappedbounds(mode::Projected, dim) = mappedbounds(mappedcrs(mode), mode, dim)
mappedbounds(mappedcrs::Nothing, mode::Projected, dim) =
    error("No mappedcrs attached to $(name(dim)) dimension")
mappedbounds(mappedcrs::GeoFormat, mode::Projected, dim) =
    reproject(crs(mode), mappedcrs, dim, bounds(dim))

mappedindex(mode::Projected, dim) = mappedindex(mappedcrs(mode), mode, dim)
mappedindex(mappedcrs::Nothing, mode::Projected, dim) =
    error("No mappedcrs attached to $(name(dim)) dimension")
mappedindex(mappedcrs::GeoFormat, mode::Projected, dim) =
    reproject(crs(dim), mappedcrs, dim, index(dim))


text:
## XXX We should be hooking to the pyconvert methods XXX
## we need Py -> Sym (asSymbolic)
## we need Sym -> Py unSym
asSymbolic(x::Sym) = x
function asSymbolic(x::Py)

    cls = x.__class__

    if pyconvert(Bool, cls == pybuiltins.list)
        return asSymbolic.(x)
    elseif pyconvert(Bool, cls == pybuiltins.dict)
        return Dict(asSymbolic(k) => asSymbolic(v) for (k,v) ∈ pyconvert(PyDict, x))
    elseif pyconvert(Bool, cls == pybuiltins.tuple)
        return Tuple(asSymbolic(xᵢ) for xᵢ ∈ x)
    elseif pyconvert(Bool, cls == pybuiltins.set)
        return Set(asSymbolic(xᵢ) for xᵢ ∈ x)
    # elseif pyconvert(Bool, cls == pybuiltins.list)
    end

    Bool(pygetattr(x, "is_FiniteSet", false)) && return Set(asSymbolic(xᵢ) for xᵢ ∈ x)

    if Bool(pygetattr(x, "is_Matrix", false))
        sz = pyconvert(Tuple, x.shape)
        !isa(sz[1], Real) && return Sym(x) # MatrixSymbol special case
        return [asSymbolic(x.__getitem__((i-1, j-1))) for i ∈ 1:sz[1], j ∈ 1:sz[2]]
    end
    Sym(x)
end
asSymbolic(x::String) = x
asSymbolic(x) = (@show x; Sym(pyconvert(Py, x)))
↑(args...; kwargs...) = asSymbolic(args...; kwargs...)

# used to pass arguments down into calls
unSym(x) = x
unSym(x::Sym) = getfield(x, :py)
unSym(x::Tuple) = unSym.(x)
unSym(x::Vector) = Tuple(unSym.(x)) # <-- is this an issue?
unSym(x::Matrix) = sympy.py.Matrix(Tuple(↓(mᵢ for mᵢ ∈ eachrow(x)))) # mutable dense matrix
unSym(x::Dict) = convert(PyDict,Dict(unSym(k) => unSym(v) for (k,v) ∈ x)) # AbstractDict?
unSym(x::Set) = sympy.py.Set(SymPyCall.unSym(collect(x)))
unSym(x::Irrational{:π}) = unSym(sympy.pi)
unSymkwargs(kw) = (k=>unSym(v) for (k,v) ∈ kw)
↓(args...; kwargs...) = unSym(args...; kwargs...)

#PythonCall.getpy(x::Sym) = ↓(x)
PythonCall.getpy(x::NTuple{N,T}) where {N, T <: SymbolicObject} = ↓(x)
PythonCall.getpy(x::Vector{T}) where {T <: SymbolicObject} = unSym.(x)
PythonCall.getpy(x::Matrix{T}) where {T <: SymbolicObject} = unSym.(x)
PythonCall.getpy(x::Set{T}) where {T <: SymbolicObject} = unSym(x)

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

# generic method pattern: (pymodule, pymethod, juliamodule, juliamethod)
# these don't need to be exported
# in SymPy we use getmembers to generate this; not sure what is better
#const
generic_methods = (
    (:sympy, :cos, :Base, :cos),
    (:sympy, :sin, :Base, :sin),
    (:sympy, :tan, :Base, :tan),
    (:sympy, :sec, :Base, :sec),
    (:sympy, :csc, :Base, :csc),
    (:sympy, :cot, :Base, :cot),
    (:sympy, :acos, :Base, :acos),
    (:sympy, :asin, :Base, :asin),
#    (:sympy, :atan, :Base, :atan),
    (:sympy, :asec, :Base, :asec),
    (:sympy, :acsc, :Base, :acsc),
    (:sympy, :acot, :Base, :acot),
    #
    (:sympy, :cosh, :Base, :cosh),
    (:sympy, :sinh, :Base, :sinh),
    (:sympy, :tanh, :Base, :tanh),
    (:sympy, :sech, :Base, :sech),
    (:sympy, :csch, :Base, :csch),
    (:sympy, :coth, :Base, :coth),
    (:sympy, :acosh, :Base, :acosh),
    (:sympy, :asinh, :Base, :asinh),
    (:sympy, :atanh, :Base, :atanh),
#    (:sympy, :asech, :Base, :asech),
#    (:sympy, :acsch, :Base, :acsch),
    (:sympy, :acoth, :Base, :acoth),
    #
    (:sympy, :sqrt, :Base, :sqrt),
    (:sympy, :exp, :Base, :exp),
    (:sympy, :log, :Base, :log),
    #
    (:sympy, :Mod, :Base, :mod),
    (:sympy, :floor, :Base, :floor),
    (:sympy, :ceiling, :Base, :ceil),
    #
    (:sympy, :numer, :Base, :numerator),
    (:sympy, :denom, :Base, :denominator),
    (:sympy, :Max, :Base, :max),
    (:sympy, :Abs, :Base, :abs),
    (:sympy, :Min, :Base, :min),

    (:sympy, :re, :Base, :real),
    (:sympy, :im, :Base, :imag),

    # solve
    (:sympy, :solve, :CommonSolve, :solve),

    # Eq
    (:sympy, :Eq, :CommonEq, :Eq),
    (:sympy, :Lt, :CommonEq, :Lt),
    (:sympy, :Le, :CommonEq, :Le),
    (:sympy, :Ne, :CommonEq, :Ne),
    (:sympy, :Ge, :CommonEq, :Ge),
    (:sympy, :Gt, :CommonEq, :Gt),

    # diff
    (:sympy, :diff, :CommonEq, :diff),

    (:sympy, :collect, :Base, :collect),

    # SpecialFunctions
    (:sympy, :airyai , :SpecialFunctions, :airyai),
    (:sympy, :airyaiprime , :SpecialFunctions, :airyaiprime),
    (:sympy, :airybi , :SpecialFunctions, :airybi),
    (:sympy, :besseli , :SpecialFunctions, :besseli),
    (:sympy, :besselj , :SpecialFunctions, :besselj),
    (:sympy, :besselk , :SpecialFunctions, :besselk),
    (:sympy, :bessely , :SpecialFunctions, :bessely),
    (:sympy, :beta , :SpecialFunctions, :beta),
    (:sympy, :erf , :SpecialFunctions, :erf),
    (:sympy, :erfc , :SpecialFunctions, :erfc),
    (:sympy, :erfi , :SpecialFunctions, :erfi),
    (:sympy, :erfinv , :SpecialFunctions, :erfinv),
    (:sympy, :erfcinv , :SpecialFunctions, :erfcinv),
    (:sympy, :gamma , :SpecialFunctions, :gamma),
    (:sympy, :digamma , :SpecialFunctions, :digamma),
    (:sympy, :polygamma , :SpecialFunctions, :polygamma),
    (:sympy, :hankel1, :SpecialFunctions, :hankelh1),
    (:sympy, :hankel2, :SpecialFunctions, :hankelh2),
    (:sympy, :zeta , :SpecialFunctions, :zeta),
)

# comment, Py(...) speeds things up a bit.
for (pmod, pmeth, jmod, jmeth) ∈ generic_methods
    @eval begin
        ($(jmod).$(jmeth))(x::Sym, args...; kwargs...) =
            ↑(getpy($(pmod)).$(pmeth)(↓(x), ↓(args)...; unSymkwargs(kwargs)...))
    end
end

## --------------------------------------------------
# pmod, pmeth, meth
#const
new_exported_methods = (
    (:sympy, :simplify, :simplify),
    (:sympy, :expand_trig, :expand_trig),
    (:sympy, :expand, :expand),
    (:sympy, :together, :together),
    (:sympy, :apart, :apart),
    (:sympy, :factor, :factor),
    (:sympy, :cancel, :cancel),
    (:sympy, :degree, :degree),
    (:sympy, :integrate, :integrate),
    (:sympy, :real_roots, :real_roots),
    (:sympy, :roots, :roots),
    (:sympy, :nroots, :nroots),
    (:sympy, :dsolve, :dsolve),
    (:sympy, :nsolve, :nsolve),
    (:sympy, :linsolve, :linsolve),
    (:sympy, :nonlinsolve, :nonlinsolve),
    (:sympy, :solveset, :solveset),

    (:sympy, :series, :series),
    (:sympy, :summation, :summation),
    (:sympy, :hessian, :hessian),
)

for (pmod, pmeth, jmeth) ∈ new_exported_methods
    @eval begin
        ($(jmeth))(x::Sym, args...; kwargs...) =
            ↑(getpy($(pmod)).$(pmeth)(↓(x), ↓(args)...; unSymkwargs(kwargs)...))
        export $(jmeth)
    end
end

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

object_methods = (
    (:conjugate, :Base, :conj),
)

for (ometh, jmod, jmeth) ∈ object_methods
    @eval begin
        ($(jmod).$(jmeth))(x::Sym, args...; kwargs...) =
            ↑(getpy(x).$(ometh)(↓(args)...; unSymkwargs(kwargs)...))
    end
end

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
using ACE, Test, LinearAlgebra
using ACE: State, CylindricalBondEnvelope, EllipsoidBondEnvelope
using StaticArrays
using ACEbase.Testing: print_tf

using ACE.Random: rand_rot, rand_refl
## Use a specific example to test the Cylindrical cutoff

@info("Testing Cylindrical Bond Envelope")

r0cut = 8.0
rcut = 4.0
zcut = 2.0
env = CylindricalBondEnvelope(r0cut, rcut, zcut)

@info("Test :bond")

for i = 1:30
    local X
    rr = rand(SVector{3, Float64}) * 2 * r0cut / sqrt(3)
    r = norm(rr)
    X = State(rr = rr, rr0 = rr, be=:bond)
    print_tf(@test( filter(env, X) == (r <= r0cut) ))
end
println()

##

@info ("Test :env")

r0 = SA[r0cut/2, 0.0, 0.0]
r_centre = r0 / 2
for i = 1:30
    local X
    rr = rand(SVector{3, Float64}) * env.rcut + r_centre
    X = State(rr = rr, rr0 = r0, be=:env)
    z = rr[1] - r_centre[1]
    r = sqrt(rr[2]^2 + rr[3]^2)
    filt = (abs(z) <= env.zcut+r_centre[1]) * (r <= env.rcut)
    print_tf(@test( filter(env, X) == (filt != 0) ))

    zeff = env.zcut + norm(r_centre)
    val = ((r/env.rcut)^2 - 1)^env.pr * ( (z/zeff)^2 - 1 )^env.pz * filt
    print_tf(@test( ACE._inner_evaluate(env, X) ≈ val ))
end
println()

#%%

@info("Testing Ellipsoid Bond Envelope")

r0cut = 2.0
rcut = 1.0
zcut = 2.0
for lambda = [0,.5,.6,1]
    for floppy = [false, true]
        local env
        env = EllipsoidBondEnvelope(r0cut, rcut, zcut;floppy=floppy, λ= lambda)

        @info("Test :bond", floppy, lambda)

        for i = 1:30
            local X
            rr = rand(SVector{3, Float64}) * 2 * r0cut / sqrt(3)
            r = norm(rr)
            X = State(rr = rr, rr0 = rr, be=:bond)
            print_tf(@test( filter(env, X) == (r <= r0cut) ))
        end
        println()

        ##

        @info ("Test :env")
        local r0, r_centre, X
        r0 = @SVector [r0cut/2, 0.0, 0.0]
        r_centre = r0 * env.λ
        for i = 1:30
            rr = rand(SVector{3, Float64}) * env.rcut + r_centre
            X = State(rr = rr, rr0 = r0, be=:env)
            z = rr[1] - r_centre[1]
            r = sqrt(rr[2]^2 + rr[3]^2)
            zeff = env.zcut + env.floppy*norm(r_centre)

            filt = (((z/ zeff)^2 + (r/env.rcut)^2) <=1)

            #@show filt
            #@show filter(env, X)
            #print("------------------- \n")
            print_tf(@test( filter(env, X) == (filt != 0) ))

            val = ( (z/zeff)^2 +  (r/env.rcut)^2 - 1.0)^env.pr * filt
            #@show val
            #@show ACE._inner_evaluate(env, X)
            #print("------------------- \n")
            print_tf(@test( ACE._inner_evaluate(env, X) ≈ val ))
        end
        println()
    end
end


text:
@testset "Multi-Head Attention Layer" begin
    num_nodes = 4
    num_edges = 10
    in_dim = 16
    out_dim = 16
    num_heads = 2
    g = SimpleDiGraph(num_nodes)
    add_edge!(g, 1, 1)
    add_edge!(g, 2, 2)
    add_edge!(g, 3, 3)
    add_edge!(g, 4, 4)
    add_edge!(g, 1, 2)
    add_edge!(g, 1, 3)
    add_edge!(g, 1, 4)
    add_edge!(g, 2, 3)
    add_edge!(g, 2, 4)
    add_edge!(g, 3, 4)
    node_features = rand(in_dim, num_nodes)    # [feature_dim, num_num_nodes]
    edge_features = rand(out_dim, num_edges)   # [feature_dim, num_num_edges]
    graph = FeaturedGraph(g, nf=node_features, ef=edge_features, directed=:directed)

    layer = MultiHeadAttentionLayer(in_dim, out_dim ÷ num_heads, num_heads)
    node_feats, edge_feats = layer(graph)
    @test size(node_feats) == (out_dim, num_nodes)
    @test size(edge_feats) == (out_dim, num_edges)

    # Test the gradients
    y = rand(out_dim, num_nodes + num_edges)
    grads = gradient(params(layer)) do
        node_feats, edge_feats = layer(graph)
        logits = cat(node_feats, edge_feats, dims=2)
        return Flux.logitcrossentropy(logits, y)
    end
    for param in params(layer)
        @test size(param) == size(grads[param])
    end

end

max_stars_repo_name: DwaraknathT/Flux-Graph-Transformer

max_issues_repo_name: DwaraknathT/Flux-Graph-Transformer

max_forks_repo_name: DwaraknathT/Flux-Graph-Transformer


text:
# Evaluation methods for additively separable functions

import Base.Cartesian.lreplace

using JuMP
using Memoize

"""
Return the array of negation of expressions.

# Input
- `exprs::Array`: An array of expressions.
"""
function negate_exprs(exprs::Array{Expr, 1})::Array{Expr, 1}
    negated_exprs = Array{Expr, 1}()
    for expr in exprs
        push!(negated_exprs, Expr(:call, :-, expr))
    end
    return negated_exprs
end

"""
Return the variable indices involved in an expression.
"""
# NOTE: The extracted expressions are already flattened out by MathProgBase, i.e., multidimentional variables are represented as 1 dimensional
# NOTE: Also, all variable names are renamed as x automatically
function get_vars_expr(expr::Expr)::IntSet
    vars = IntSet()
    exprs_to_process::Array{Expr, 1} = [expr]               #stores remaining (sub)expressions in expr that need to be preocessed
    while !isempty(exprs_to_process)
        cur_expr::Expr = pop!(exprs_to_process)
        if cur_expr.head == :call           #if the current expression is a "call" type, e.g., x + y or 2*x, ...
            for i::Int in 2:length(cur_expr.args)
                if isa(cur_expr.args[i], Expr)  #if an argument is exprresion, it is added to the list to be processed
                    push!(exprs_to_process, cur_expr.args[i])
                end
            end
        elseif cur_expr.head == :ref        #if the current expression is a "reference" type, e.g., x[i]
            if cur_expr.args[1] != :x
                error("Variables in expression not supported (", expr, ")")
            end
            push!(vars, cur_expr.args[2]) # found variable, store index
        else
            error("Expression in constraint not supported (", expr, ")")
        end
    end
    return vars
end

"""
Decompose a given expression into separable terms.

# Output:
- An array that stores decomposed parts (expressions) of a given expression.
- The constant value of the expression
These parts are separated by ``+`` sign, e.g., ``3x[1] - x[1]^2 + sin(x[2]) - 4 <= 0``, will be stored in an array with 3 elements for the lhs and a constant -4
If the constraint is not additively separable, throws error.
"""

function separate_additively(expr::Expr)::Tuple{Array{Expr, 1}, Float64}
    # Examine the expression tree: first split (in)equalities and then check if
    # each term separated by a + or - uses only one type of variable
    # println(expr.head, " -- ", expr.args, " -- ", expr.typ)
    constant_val::Float64 = 0
    if expr.head == :call
        # println("Call")
        # (In)equality or non-unary -
        if expr.args[1] in [:(<=), :(==), :(>=)] || (expr.args[1] == :- && length(expr.args) >= 3)
            # If (in)equality, then move right-hand side to left-hand side by negating it
            # If subtraction, then negate second expression
            exprs1 = Array{Expr, 1}()
            if isa(expr.args[2], Float64)
                constant_val += expr.args[2]
            else
                a, b = separate_additively(expr.args[2])
                append!(exprs1, a)
                constant_val += b
            end
            if isa(expr.args[3], Float64)
                constant_val -= expr.args[3]
            else
                a, b = separate_additively(expr.args[3])
                append!(exprs1, negate_exprs(a))
                constant_val -= b
            end
            return exprs1, constant_val
        # Non-unary +
        elseif expr.args[1] in [:+] && length(expr.args) >= 3
            exprs = Array{Expr, 1}()
            for i in 2:length(expr.args)
                if isa(expr.args[i], Float64)
                    constant_val += expr.args[3]
                else
                    a, b = separate_additively(expr.args[i])
                    append!(exprs, a)
                    constant_val += b
                end
            end
            return exprs, constant_val
        # Unary +
        elseif expr.args[1] in [:+] # unary; binary operator already checked
            a, b = separate_additively(expr.args[2])
            return a, b
        # Unary -
        elseif expr.args[1] in [:-] # unary; binary operator already checked
            a, b = separate_additively(expr.args[2])
            return negate_exprs(a), -b
        # Process expression as a single component; must only contain one type of variable
        else
            vars::IntSet = get_vars_expr(expr)
            if length(vars) >= 2
                error("Non-additively separable constraint in model")
            end
            return [expr], 0
        end
    elseif expr.head == :comparison
        error("Only one-sided constraints are supported")  # to make DD construction easier
    elseif expr.head == :ref
        return [expr], 0  # isolated variable
    else
        error("Expression in constraint not supported (", expr, ")")
    end
end

"""
Return the sense of a given expression in form of an (in)equality.
Returns error if the expression is not (in)equality
"""

function get_constraint_sense(expr::Expr)::Symbol
    if expr.head == :call && expr.args[1] in [:(<=), :(==), :(>=)]
        return expr.args[1]            #Symbol for the (in)equality sense
    else
        error("Expression not in constraint form (", expr, ")")
    end
end

#************************************************************************
#*************************************************************************
#*************************************************************************

# Evaluate a univariate expression replacing the variable by the given value.
@memoize function eval_var(expr::Expr, val::Float64)::Float64
    # Replace references by the value
    if expr.head == :ref        #if the expression is composed of a single variable assigns the value
        return val
    else                        #otherwise decomposes the expression down to variables
        copy_expr = copy(expr)
        exprs_to_process = [copy_expr]
        while !isempty(exprs_to_process)
            cur_expr = pop!(exprs_to_process)
            if cur_expr.head == :call
                for i in 2:length(cur_expr.args)
                    if isa(cur_expr.args[i], Expr)
                        if cur_expr.args[i].head == :ref
                            cur_expr.args[i] = val          #replacing the variable with its value
                        else
                            push!(exprs_to_process, cur_expr.args[i])
                        end
                    end
                end
            else
                error("Expression in constraint not supported (", expr, ")")
            end
        end
        return eval(copy_expr)
    end
end

# TODO: Another method to evaluate a given expression, where we allow several variables (for nonseparable case!)

"""
Oracle that yields functional specifications of a constraint.
"""
struct SeparableFunctionEvaluation
#    separate_exprs::Array{Expr}        # each element stores the univariate expression for a variable
    single_function::Function           # evaluate the whole multivariate function in the constraint at any point (without the constant part)
    separate_functions::Array{Function, 1} # each element stores the univariate function for a variable
    constant::Float64                   # constant term of function
end

"""
Return decomposition terms of constraints of a model.

# Output
- `evals`: An array of instances that stores function specification for constraints
"""
function evaluate_separable_constraints(m::Model)::Array{SeparableFunctionEvaluation}
    d = JuMP.NLPEvaluator(m)        #obtains an NLP evaluator object from a JuMP model, which enables query of function characteristics in the model
    MathProgBase.initialize(d, [:ExprGraph])    #enables query of the functions of the model in their expression graph form

    # Gather parts for separable constraints
    nlinearcon = MathProgBase.numlinconstr(m)           # number of linear constraints (first ordered in NLPEvaluator)
    nquadcon = MathProgBase.numquadconstr(m)            # number of quadratic constraints (second ordered in NLPEvaluator)
    nlpcon = JuMP.numnlconstr(m)                        # number of nonlinear constraints, besides quadratic (third ordered in NLPEvaluator)
    nconstrs = nquadcon + nlpcon                    # total number of nonlinear constraints

    nvars = MathProgBase.numvar(m)
    evals = Array{SeparableFunctionEvaluation}(  nconstrs)

    for i in 1:nconstrs
        # NOTE: The extracted expressions are already flattened out by MathProgBase, i.e., multidimentional variables are represented as 1 dimensional
        # NOTE: Also, all variable names are renamed as x automatically
        expr = MathProgBase.constr_expr(d, i + nlinearcon)       #gets the expression of the ith constraint among the nonlinear ones
        parts, constant = separate_additively(expr)         #gets the array of decomposed part of the constraint

        # Put all constraints in f(x) <= 0 form by negating the lhs when >=.
        con_sense = get_constraint_sense(expr)
        if con_sense == :>=
            parts = negate_exprs(parts)
            constant = -constant
        # elseif con_sense == :==           #need to implement for == case: create two constraints with <=...
        end

        # Organize parts per variables
        var_parts = Array{Expr}(  nvars)
        for p in parts
            v = first(get_vars_expr(p))
            if !isassigned(var_parts, v)
                var_parts[v] = p
            else
                var_parts[v] = Expr(:call, :+, var_parts[v], p)
            end
        end

        # Create function corresponding to each univariate expression
        function_parts = Array{Function, 1}(  nvars)
        for j in 1:nvars
            if isassigned(var_parts, j)
                function_parts[j] = (x::Real) -> eval_var(var_parts[j], x)
            end
        end

        # Create a single multivariate function representing lhs of a constraint
        single_function = function (x::Array{T} where T<:Real)
            @assert(length(x) == nvars)
            sum_v = sum(function_parts[j](x[j]) for j in 1:nvars if isassigned(var_parts, j))
            return sum_v
        end

        evals[i] = SeparableFunctionEvaluation(single_function, function_parts, constant)

    end

    return evals
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
function calc_drag(buttons, drag, indrag, tracked_mousebutton)
    # only track if still the same button is pressed
    if length(buttons) == 1 && (!indrag[] || tracked_mousebutton[] == first(buttons))
        if !indrag[]
            tracked_mousebutton[] = first(buttons); indrag[] = true
            drag[] = Mouse.down # just started, so dragging is still false
            return drag[]
        else
            drag[] = Mouse.pressed
            return drag[]
        end
    end
    # already on notpressed, no need for update
    if drag[] != Mouse.notpressed
        drag[] = indrag[] ? Mouse.up : Mouse.notpressed
    end
    indrag[] = false
    return drag[]
end

function mousedrag(scene::Scene, native_window)
    indrag = RefValue(false)
    tracked_mousebutton = RefValue(Mouse.left)
    drag = RefValue(Mouse.notpressed)
    events = scene.events
    onany(events.mouseposition, events.mousebuttons) do mp, buttons
        d = calc_drag(buttons, drag, indrag, tracked_mousebutton)
        if (d == Mouse.pressed) || (d != events.mousedrag[])
            events.mousedrag[] = d
        end
        return
    end
    return
end

function disconnect!(window::AbstractScreen, signal)
    disconnect!(to_native(window), signal)
end
window_area(scene, native_window) = not_implemented_for(native_window)
window_open(scene, native_window) = not_implemented_for(native_window)
mouse_buttons(scene, native_window) = not_implemented_for(native_window)
mouse_position(scene, native_window) = not_implemented_for(native_window)
mousedrag(scene, native_window) = not_implemented_for(native_window)
scroll(scene, native_window) = not_implemented_for(native_window)
keyboard_buttons(scene, native_window) = not_implemented_for(native_window)
unicode_input(scene, native_window) = not_implemented_for(native_window)
dropped_files(scene, native_window) = not_implemented_for(native_window)
hasfocus(scene, native_window) = not_implemented_for(native_window)
entered_window(scene, native_window) = not_implemented_for(native_window)

function register_callbacks(scene::Scene, native_window)

    window_area(scene, native_window)
    window_open(scene, native_window)
    mouse_buttons(scene, native_window)
    mouse_position(scene, native_window)
    mousedrag(scene, native_window)
    scroll(scene, native_window)
    keyboard_buttons(scene, native_window)
    unicode_input(scene, native_window)
    dropped_files(scene, native_window)
    hasfocus(scene, native_window)
    entered_window(scene, native_window)

end

button_key(x::Type{T}) where {T} = error("Must be a keyboard or mouse button. Found: $T")
button_key(x::Type{Keyboard.Button}) = :keyboardbuttons
button_key(x::Type{Mouse.Button}) = :mousebuttons
button_key(x::Set{T}) where {T} = button_key(T)
button_key(x::T) where {T} = button_key(T)

"""
returns true if `button` is pressed in scene[:mousebuttons or :keyboardbuttons]
You can use nothing, to indicate it should always return true
"""
function ispressed(scene::SceneLike, button::Union{Vector, Tuple})
    all(x-> ispressed(scene, x), button)
end

# TODO this is a bit shady, but maybe a nice api!
# So you can use void whenever you don't care what is pressed
ispressed(scene::SceneLike, ::Nothing) = true

function ispressed(buttons::Set{T}, button) where T <: Union{Keyboard.Button, Mouse.Button}
    if isa(button, Set)
        return buttons == button
    else
        return length(buttons) == 1 && first(buttons) == button
    end
end
function ispressed(scene::SceneLike, button)
    buttons = getfield(events(scene), button_key(button))[]
    ispressed(buttons, button)
end

"""
Picks a mouse position
"""
function pick end

"""
    onpick(func, plot)
Calls `func` if one clicks on `plot`
"""
function onpick end


text:
function geneangles(gene::AbstractGene, chrlength::Int, ORI::Int = 1)
    loc = locus(gene)
    geneangles(loc.position.start, loc.position.stop, chrlength, ORI)
end
function geneangles(start::Number, stop::Number, chrlength::Int, ORI::Int = 1)
    a1 = mod2pi((2π * (start-ORI-1) / chrlength)) - π/2
    a2 = mod2pi((2π * (stop-ORI-1) / chrlength)) - π/2
    amid = mod2pi((2π * (mean([start, stop])-ORI-1) / chrlength)) - π/2
    a1, a2, amid
end

"""
    function drawplasmid(chr; kwargs...)

Draw a circular plasmid map. CDSs are displayed as boxes on the outside and inside of the plasmid, for genes on the + and - strands, respectively.
Supported kwargs are:
* outfile           Path to output file. Supported file types are PNG, SVG, PDF, and EPS.
* drawingsize       Size of the figure, can be a `Tuple{Int,Int}`, or a `String` such as "1000x1000".
* title             Name of the plasmid.
* ORI               The position of the ORI (defaults to 1).
* colors            A `Dict`, with `Gene`s as keys and colors as values (defaults to `defaultgenecolor`).
* defaultgenecolor  Defaults to :lightgrey.
* annotations       A `Vector` containing annotations for each CDS.
* internalcolors    A `Dict` with `UnitRange`s (representing nucleotide positions) as keys, and colors as values. Can be used to color parts of genes, or intergenic regions.
* genethickness     Size of the gene boxes.
* textoffset        Distance from the plasmid where annotations should be drawn.
* figureoffset      A `Tuple{Int,Int}` with offsets in the x and y axes, respectively. Useful if there are more annotations on one side than the other.
* genegroups        A `Dict` with `UnitRanges` as keys (representing gene numbers), and `NamedTuples` as values. The tuple must have the field `text`, which contains the annotations for the group, and can additionally take `placement` (:outer or :inner), `offset`, `textoffset`, and `thickness`.
"""
function drawplasmid(chr; outfile = "plasmid.png", offset = 0, colors = Dict(),
                     annotations = [], drawingsize = (1000, 1000),
                     internalcolors = Dict(),
                     defaultgenecolor = RGB(.95, .95, .95),
                     ORI = 1,
                     genethickness = 20, textoffset = 35, figureoffset = (0, 0),
                     genegroups = [],
                     title = "")
    xmax, ymax = drawingdimensions(drawingsize)
    chrlength = length(chr.sequence)
    d = min(xmax, ymax) / 2.5
    drawing = Drawing(xmax, ymax, outfile)
    background("white")
    origin()
    translate(figureoffset...)
    fontsize(28)
    textcentered(title, Point(0, -(ymax / 2.2)))
    circle(O, d, :stroke)
    fontsize(22)
    textcentered(string(chrlength) * " bp")
    fontsize(14)
    genes = @genes(chr, CDS)
    ## draw genes
    for (i, gene) in enumerate(genes)
        loc = locus(gene).position
        a1, a2, amid = geneangles(gene, chrlength, ORI)
        if isempty(annotations)
            genetext = gene.product
        else
            genetext = annotations[i]
        end
        if iscomplement(gene)
            gsave()
                setcolor(get(colors, get(gene, :locus_tag, ""), defaultgenecolor))
                sector(d - genethickness, d, a1, a2, :fillpreserve)
            grestore()
            strokepath()
            textradius = d - textoffset
            if -π/2 < amid < π/2
                text(genetext, Point(textradius * cos(amid), textradius * sin(amid)), halign = :right, valign = :center)
            elseif π/2 <= amid <= 3π/2
                text(genetext, Point(textradius * cos(amid), textradius * sin(amid)), halign = :left, valign = :center)
            end
            ## internal colors
            for (r, c) in get(internalcolors, gene, [])
                a1, a2, amid = geneangles(loc.stop - r.stop, loc.stop - r.start, chrlength, ORI)
                gsave()
                    setcolor(c)
                    sector(d-genethickness, d, a1, a2, :fill)
                grestore()
            end
        else
            gsave()
            setcolor(get(colors, get(gene, :locus_tag, ""), defaultgenecolor))
            sector(d, d + genethickness, a1, a2, :fillpreserve)
            grestore()
            strokepath()
            textradius = d + textoffset
            if -π/2 < amid < π/2
                text(genetext, Point(textradius * cos(amid), textradius * sin(amid)), halign = :left, valign = :center)
            elseif π/2 <= amid <= 3π/2
                text(genetext, Point(textradius * cos(amid), textradius * sin(amid)), halign = :right, valign = :center)
            end
            ## internal colors
            for (r, c) in get(internalcolors, gene, [])
                a1, a2, amid = geneangles(loc.start + r.start, loc.start + r.stop, chrlength, ORI)
                gsave()
                    setcolor(c)
                    sector(d, d+genethickness, a1, a2, :fill)
                grestore()
            end
        end
    end
    ## draw group annotations
    fontsize(18)
    for (i, kv) in enumerate(genegroups)
        k, v = kv
        geneids = k
        annotation = get(v, :text, "")
        offset = get(v, :offset, 0.1d)
        textoffset = get(v, :textoffset, 15)
        thickness = get(v, :thickness, 2)
        placement = get(v, :placement, :outer)
        a1, a2, amid = geneangles(locus(genes[geneids][1]).position.start, locus(genes[geneids][end]).position.stop, chrlength, ORI)
        if placement == :inner
            sector(d-offset, d-offset+thickness, a1, a2, :fill)
            textcurvecentered(annotation, amid, d-offset-textoffset, O)
        elseif placement == :outer
            sector(d+offset, d+offset+thickness, a1, a2, :fill)
            textcurvecentered(annotation, amid, d+offset+textoffset, O)
        end
    end
    finish()
end


text:
#InspectDR: Event handling with Gtk
#-------------------------------------------------------------------------------

import Gtk: GConstants.GdkScrollDirection

#==Constants
===============================================================================#
#Symbols not currently supported by Gtk.GdkKeySyms:
const GdkKeySyms_Plus = Gtk.GConstants.GDK_KEY_KP_Add
const GdkKeySyms_Minus = Gtk.GConstants.GDK_KEY_KP_Subtract

#==Main types
===============================================================================#
struct ISNormal <: InputState; end #Default state

#==Helper functions
===============================================================================#
modifiers_pressed(eventstate, modmask) = (modmask == (MODIFIERS_SUPPORTED & eventstate))

#Returns 0 or strip index:
#TODO: refer to graphbblist only - instead of pwidget???
function hittest(pwidget::PlotWidget, x::Float64, y::Float64)
	maplist = pwidget.rplot.strips #WANTCONST
	for i = 1:length(maplist)
		if isinside(maplist[i].bb, x, y)
			return i
		end
	end
	return 0
end

#==Wrapper functions
===============================================================================#
raiseevent_plothover(pwidget::PlotWidget, event::Gtk.GdkEventMotion) =
	raiseevent(pwidget.eh_plothover, pwidget, event.x, event.y)

#==Low-level event handlers
===============================================================================#
function handleevent_plothover(pwidget::PlotWidget, event::Gtk.GdkEventMotion, istrip::Int)
	x, y = event.x, event.y

	if istrip > 0
		#TODO: do we want "read"able coords?
		pwidget.mouseover.pos = dev2axis(pwidget.rplot.strips[istrip], Point2D(x, y))
	else
		pwidget.mouseover.pos = nothing
	end

	raiseevent_plothover(pwidget, event) #Inform owner
	return
end

#Update pwidget.mouseover & raise event to owner
function handleevent_plothover(pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	x, y = event.x, event.y
	istrip = hittest(pwidget, x, y)
	pwidget.mouseover.istrip = istrip
	handleevent_plothover(pwidget, event, istrip)
	return
end

#Last line


text:
# although you can create and manage your own *long-term* development packages
# the way Engheim suggests i.e. using your own directory, the Julia documentation
# suggests a different and better solution as a best practice:
# https://docs.julialang.org/en/v1/manual/workflow-tips/#Revise-based-workflows
# In your repl:
# pkg > add PkgTemplates
# julia> using PkgTemplates
# if your Github user is defined in git global config, you don't have to add the parameter.
# julia> t = Template(user=",your Github username>")
# julia> t("Foobar")
# As the documentation indicates this stores your project in a directory ~/.julia/dev.
# Engheim points out his reason for preferring his own directory:
# "Should your Julia installation get messed up and you want to start from scratch,
# it is tempting to erase the whole ~/.julia, but this can have catastrophic results
# if you got the habit of putting your personal projects in there.”
# However if you are in the habit of religiously pushing code to Github et al (as you should)
# there is no reason this should be catastrophic. Just check out the repository again.
# Another concern might be a disk size limit on the home directory. Whatever the case, you can always symlink
# ~/.julia/dev to another directory if you can't/don't want to store your code in your home directory.
# Or you can move your generated package directory and symlink it inside dev. Finally you can create
# a non hidden symlink to the dev directory, to easily access your packages from the Finder et al.
# NB: Engheim does discuss PkgTemplates later in the context of publishing on JuliaHub
# but it seems simple enough and useful for all projects.
# The package template does lots of heavy lifting for you, including e.g.
# creating the git repository for the package and much more (cf. the Foobarfolder in this
# repository, although for obvious reasons didnt add the .git directory)
# see PkgTemplate doco for more info:
# https://invenia.github.io/PkgTemplates.jl/stable/
# Also, it seems that adding directories to the LOAD_PATH env variable
# is not a best practice. So how do you get the Julia REPL to recognize your
# packages in the dev directory? Simple:
# (@v1.7) pkg> dev Foobar
# When you have questions you should definitely checkout Julia's friendly
# discourse. That's how I found out the above:
# https://discourse.julialang.org/t/best-practice-for-generating-long-term-packages/72853
# You can learn all about dev in the Pk doco. The whole document is worth reading:
# https://pkgdocs.julialang.org/v1/managing-packages/#developing
# One last point. Engheim seems to imply that if you add your directory to the
# LOAD_PATH, then when you do:
# (@v1.7) pkg> activate "Foobar"
# it will somehow know you are taling about *your* Foobar. My testing indicates
# that generates a Foobar package inside your current directory. Hence what
# you actually need to do is provide the full path (no bash shortcuts allowed,
# even though the response uses one!):
# (@v1.7) pkg> activate "/Users/aronet/.julia/dev/Foobar"
#  Activating project at `~/.julia/dev/Foobar`
# I mounted the ~/.julia/dev/Foobar directory into this Github repository using
# bindfs. First install macfuse as a brew cask. Then download source of bindfs
# and build and install using gmake (brew version of make). Then mkdir the mount
# directory (in this case "Foobar") and then run
# bindfs <fully qualified source directory> <fully qualified target directory>

module Foobar

# Write your package code here.
export greet

using Dates

greet() = print("Hello Earth, it's a beautiful ", dayname(today()))

end
# Revise example
# add "Revise" to current Julia version env. Fet there by running activate unmodified
# (JuliaforBeginners) pkg> activate
#   Activating project at `~/.julia/environments/v1.7`

# (@v1.7) pkg> add Revise
#     Updating registry at `~/.julia/registries/General.toml`
#    Resolving package versions...
#    Installed LoweredCodeUtils ─ v2.1.2
#    Installed CodeTracking ───── v1.0.6
#    Installed JuliaInterpreter ─ v0.8.21
#    Installed Revise ─────────── v3.2.0
#     Updating `~/.julia/environments/v1.7/Project.toml`
#   [295af30f] + Revise v3.2.0
#     Updating `~/.julia/environments/v1.7/Manifest.toml`
#   [da1fd8a2] + CodeTracking v1.0.6
#   [aa1ae85d] + JuliaInterpreter v0.8.21
#   [6f1432cf] + LoweredCodeUtils v2.1.2
#   [ae029012] + Requires v1.1.3
#   [295af30f] + Revise v3.2.0
#   [7b1f6079] + FileWatching
# Precompiling project...
#   6 dependencies successfully precompiled in 9 seconds (19 already precompiled)

# (@v1.7) pkg> activate "/Users/aronet/.julia/dev/Foobar"
#   Activating project at `~/.julia/dev/Foobar`

# julia> using Revise

# julia> using Foobar

# julia> greet()
# Hello World, it's a beautiful Friday
# modify the code and rerun function
# julia> greet()
# Hello Earth, it's a beautiful Friday

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
# Romeo LTI pretty printing.
#
# Utility functions for human-readable representations of systems.

"""
    pprint(sys [, var="s"])

Generates pretty human-readable representation of Laplace transform of the system `sys`,
with the Laplace variable denoted by `var`.
"""
pprint(sys::SisoLtiSystem) = pprint(sys, 's')

pprint(sys::UnitSys, var::Char) = "1"
pprint(sys::ZeroSys, var::Char) = "0"

function pprint(sys::Diff{<:Number}, var::Char)
    if sys.α == 0
        return "1"
    elseif sys.α == 1
        return "$var"
    else
        return "$var^$(sys.α)"
    end
end

pprint(sys::ScaledSystem, var::Char) = "$(sys.k) * ($(pprint(sys.inner, var)))"
pprint(sys::SeriesSystem, var::Char) = "($(pprint(sys.first, var))) * ($(pprint(sys.second, var)))"
pprint(sys::ParallelSystem, var::Char) = "$(pprint(sys.first, var)) + $(pprint(sys.second, var))"
pprint(sys::RationalSystem, var::Char) = "($(pprint(sys.num, var))) / ($(pprint(sys.den, var)))"
# pprint(sys::SquareRootSystem, var::Char) = "sqrt($(pprint(sys.inner, var)))"
pprint(sys::PowerSystem, var::Char) = "($(pprint(sys.inner, var)))^$(sys.α)"

macro pprint(sys)
    return esc(:(Base.println(pprint($sys))))
end

export pprint, @pprint


text:
#sbox5.jl
"""Implement a pattern for dispatching actions based on Truth table values"""

module DispatchPattern
using Match

struct Truth
  a::Bool
  b::Bool
end

actionAB()   =println("Do A and B")
actionA!B()  =println("Do A not B")
action!AB()  =println("Do B not A")
action!A!B() =println("Do nothing")

"Pattern for dispatching based on multiple truth cases"
dispatch!(item::Truth)=@match (item.a, item.b)  begin
  (true, true)          => actionAB()
  (true, false)         => actionA!B()
  (false, true)         => action!AB()
  (false, false)        => action!A!B()
  bad                   => println("Unknown dispatch: $bad")
end

truths=[Truth(true, true),
        Truth(true, false),
        Truth(false, true),
        Truth(false, false)]

map(dispatch!, truths)

end


text:
#https://juliaobserver.com/packages/TableReader
using TableReader
using DataFrames
using TimeZones
using Dates

url = "https://realized.oxford-man.ox.ac.uk/images/oxfordmanrealizedvolatilityindices.zip"
dataframe = readcsv(`unzip -p data/oxfordmanrealizedvolatilityindices.zip oxfordmanrealizedvolatilityindices.csv `)

dataframe = DataFrames.rename(dataframe, :UNNAMED_1 => :timestring)
dateformat = Dates.DateFormat("yyyy-mm-dd HH:MM:SSzzz")
#
Dates.DateTime("2019-12-13 00:00:00+00:00", dateformat)
# 2019-12-13T00:00:00
dataframe.timestring[1]
Dates.DateTime(dataframe.timestring[1], dateformat)

dataframe.timestamp =  Dates.DateTime.(dataframe.timestring, dateformat)
symbols = unique(dataframe, :Symbol).Symbol
colnames  = names(dataframe)

coldict = collect(eachcol(dataframe, true))    # build a dict of col name and values
### OLD WAY!  see https://github.com/JuliaLang/julia/issues/17886
### namesOfIndicators =  [name for (name, value ) in filter(x -> eltype(x.second) == Float64, coldict)]

# https://realized.oxford-man.ox.ac.uk/documentation/estimators
# Available Estimators
# Code	Description
# bv	Bipower Variation (5-min)
# bv_ss	Bipower Variation (5-min Sub-sampled)
# close_price	Closing (Last) Price
# close_time	Closing Time
# medrv	Median Realized Variance (5-min)
# nobs	Number of Observations
# open_price	Opening (First) Price
# open_time	Opening Time
# open_to_close	Open to Close Return
# rk_parzen	Realized Kernel Variance (Non-Flat Parzen)
# rk_th2	Realized Kernel Variance (Tukey-Hanning(2))
# rk_twoscale	Realized Kernel Variance (Two-Scale/Bartlett)
# rsv	Realized Semi-variance (5-min)
# rsv_ss	Realized Semi-variance (5-min Sub-sampled)
# rv10	Realized Variance (10-min)
# rv10_ss	Realized Variance (10-min Sub-sampled)
# rv5	Realized Variance (5-min)
# rv5_ss	Realized Variance (5-min Sub-sampled)

#indicators are Float64 types (one seem to have missing value, fot now don't include it)
#(:medrv, Union{Missing, Float64})
namesOfIndicators = [name for (name, value) in coldict if eltype(value) == Float64]

#Using shifting arrays to deal with the deltas, note the potential to use
#regular loop also (instead of specialized panda code)
#https://discourse.julialang.org/t/create-lead-and-lag-variable-in-dataframe/8658/1

# There are multiple ways to calculate the Hurst index, below is the generalized
# version that is based on calculating the moments of the differences.

# spx['sqrt']= np.sqrt(spx['SPX2.rk'])
# spx['log_sqrt'] = np.log(spx['sqrt'])
#
# def del_Raw(q, x):
#     return [np.mean(np.abs(spx['log_sqrt'] - spx['log_sqrt'].shift(lag)) ** q)
#             for lag in x]
# note data is variance, need to take sqrt to get sigma, then log, then absolute
# value of diff at
#different lags and raised to power (for the qth moment) then mean

#  m(q,Δ)∝ Δ ** ζq    This means that the log of m(q, Δ) and ζ of q
# have linear relationship..   For each q the ζ is the slope.
# Hurst index is how the ζ changees with different q
# note  q is the q-moment of the delayed time series.

log_σ(σ_squre::Float64) = log(√σ_squ)

logσt(σ_squres) = log_σ.(σ_squres)
dataframe.rv10
dataframe.rv10
Define a striucture for the fraction information of a series
    Constructor takes a series and builds the parameters

# For start just get rv10 indicator for .SPX as timeseries

Function for timeseries (ts, qvecs,range of Delta,  )
    for each q moment
        for each delta

    return h, (for each q moment, zetaq, list of (log deltas, log m))


text:
module DummyTask

using Genie, Genie.Toolbox

"""
Description of the task here
"""
function run_task()
  # Build something great
end

end

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
using Printf

function main()
    output_dir_name  = "post_result"
    input_result_dir = "result"

    #vtk_頭の文字
    front = "# vtk DataFile Version 2.0"*"\n"
    back  = "\n"*"ASCII"*"\n"*"DATASET UNSTRUCTURED_GRID"*"\n"
    write_file(output_dir_name, input_result_dir, front, back)
end

function write_file(outdir, inresult_dir, front, back)
    x = zeros(8,3)

    inf = readdir(inresult_dir)
    for i in 1:length(inf)
        if occursin(".dat", inf[i]) == true
            dname = replace(inf[i],".dat" => "")
            out_file = replace(inf[i],".dat" => ".vtk")
            print("start writing "*out_file*"\n")

            fff=[]
            open(inresult_dir*"/"*inf[i], "r") do f
                fff=read(f,String)
            end
            fff = split(fff,"\n",keepempty=false)

            for j in 2:length(fff)
                fff[j] = replace(fff[j]," \r" => "")
                temp = split(fff[j]," ")
                x[j-1,1] = parse(Float64,temp[1])
                x[j-1,2] = parse(Float64,temp[2])
                x[j-1,3] = parse(Float64,temp[3])
            end

            write_points(x, out_file, outdir, dname, front, back)
            write_cells(out_file, outdir)
            write_result(out_file, outdir)
            print("fin writing "*out_file*"\n")
        end
    end
end

function  write_points(x_arr, out_file, outdir, dname, front, back)
    a = size(x_arr)[1]
    a_st = @sprintf("%1.0f", a)

    fff = outdir*"/"*out_file
    open(fff,"w") do f
        write(f,front)
        write(f,dname)
        write(f,back)
        write(f,"POINTS "*a_st*" float\n")
        for i in 1:a
            x = @sprintf("%7.7f", x_arr[i,1])
            y = @sprintf("%7.7f", x_arr[i,2])
            z = @sprintf("%7.7f", x_arr[i,3])

            write(f, x*" "*y*" "*z*"\n")
        end
    end
    println("fin writing points")
end

function write_cells(out_file,outdir)

    fff = outdir*"/"*out_file
    open(fff,"a") do f
        write(f, "CELLS "*"1"*" "*"9"*"\n")

        d1 = @sprintf("%1.0f", 0)
        d2 = @sprintf("%1.0f", 1)
        d3 = @sprintf("%1.0f", 3)
        d4 = @sprintf("%1.0f", 2)
        d5 = @sprintf("%1.0f", 4)
        d6 = @sprintf("%1.0f", 5)
        d7 = @sprintf("%1.0f", 7)
        d8 = @sprintf("%1.0f", 6)
        write(f, "8 "*d1*" "*d2*" "*d3*" "*d4*" "*d5*" "*d6*" "*d7*" "*d8*"\n")

        write(f, "CELL_TYPES "*"1"*"\n")

        write(f, "12\n")          #四角のみ
    end
    println("fin writing cells")
end

function write_result(out_file, outdir)
    temp1 = "CELL_DATA "*"1"*"\n"
    temp2 = "SCALARS "*"null"*" float\nLOOKUP_TABLE default\n"

    fff = outdir*"/"*out_file
    open(fff,"a") do f

        write(f, temp1)

        write(f, temp2)

        data = @sprintf("%7.7f", 1.0)
        write(f, data)
        write(f, "\n")

    end
end

# main
main()