text:
# This file contains SDESystem prototypes

import DifferentialEquations: LambaEM, SDEProblem
import UUIDs: uuid4

"""
    @def_sde_system ex

where `ex` is the expression to define to define a new AbstractSDESystem component type. The usage is as follows:
```julia
@def_sde_system mutable struct MySDESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractSDESystem
    param1::T1 = param1_default                 # optional field
    param2::T2 = param2_default                 # optional field
    param3::T3 = param3_default                 # optional field
        ⋮
    paramN::TN = paramN_default                 # optional field
    drift::DR = drift_function                  # mandatory field
    diffusion::DF = diffusion_function          # mandatory field
    readout::RO = readout_functtion             # mandatory field
    state::ST = state_default                   # mandatory field
    input::IP = input_default                   # mandatory field
    output::OP = output_default                 # mandatory field
end
```
Here, `MySDESystem` has `N` parameters. `MySDESystem` is represented by the `drift`, `diffusion` and `readout` function. `state`, `input` and `output` is the initial state, input port and output port of `MySDESystem`.

!!! warning
    `drift` must have the signature
    ```julia
    function drift((dx, x, u, t, args...; kwargs...)
        dx .= .... # update dx
    end
    ```
    and `diffusion` must have the signature
    ```julia
    function diffusion((dx, x, u, t, args...; kwargs...)
        dx .= .... # update dx
    end
    ```
    and `readout` must have the signature
    ```julia
    function readout(x, u, t)
        y = ...
        return y
    end
    ```

!!! warning
    New SDE system must be a subtype of `AbstractSDESystem` to function properly.

# Example
```julia
julia> @def_sde_system mutable struct MySDESystem{DR, DF, RO, IP, OP} <: AbstractSDESystem
           η::Float64 = 1.
           drift::DR = (dx, x, u, t) -> (dx .= x)
           diffusion::DF = (dx, x, u, t, η=η) -> (dx .= η)
           readout::RO = (x, u, t) -> x
           state::Vector{Float64} = rand(2)
           input::IP = nothing
           output::OP = Outport(2)
       end

julia> ds = MySDESystem();
```
"""
macro def_sde_system(ex)
    checksyntax(ex, :AbstractSDESystem)
    appendcommonex!(ex)
    foreach(nex -> appendex!(ex, nex), [
        :( alg::$ALG_TYPE_SYMBOL = Causal.LambaEM{true}() ),
        :( integrator::$INTEGRATOR_TYPE_SYMBOL = Causal.construct_integrator(Causal.SDEProblem, input,
            (drift, diffusion), state, t, modelargs, solverargs; alg=alg, modelkwargs=modelkwargs, solverkwargs=solverkwargs, numtaps=3) )
        ])
    quote
        Base.@kwdef $ex
    end |> esc
end

##### Define SDE system library

"""
    SDESystem(; drift, diffusion, readout, state, input, output)

Constructs a SDE system.
"""
@def_sde_system mutable struct SDESystem{DR, DF, RO, ST, IP, OP} <: AbstractSDESystem
    drift::DR
    diffusion::DF
    readout::RO
    state::ST
    input::IP
    output::OP
end

@doc raw"""
    NoisyLorenzSystem()

Constructs a noisy Lorenz system
"""
@def_sde_system mutable struct NoisyLorenzSystem{ET, DR, DF, RO, IP, OP} <: AbstractSDESystem
    σ::Float64 = 10.
    β::Float64 = 8 / 3
    ρ::Float64 = 28.
    η::ET = 1.
    γ::Float64 = 1.
    drift::DR = function lorenzdrift(dx, x, u, t, σ=σ, β=β, ρ=ρ, γ=γ)
        dx[1] = σ * (x[2] - x[1])
        dx[2] = x[1] * (ρ - x[3]) - x[2]
        dx[3] = x[1] * x[2] - β * x[3]
        dx .*= γ
    end
    diffusion::DF = (dx, x, u, t, η=η) -> (dx .= η)
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = nothing
    output::OP = Outport(3)
end

@doc raw"""
    NoisyLorenzSystem()

Constructs a noisy Lorenz system
"""
@def_sde_system mutable struct ForcedNoisyLorenzSystem{ET, CM, DR, DF, RO, IP, OP} <: AbstractSDESystem
    σ::Float64 = 10.
    β::Float64 = 8 / 3
    ρ::Float64 = 28.
    η::ET = 1.
    cplmat::CM = I(3)
    γ::Float64 = 1.
    drift::DR = function forcedlorenzdrift(dx, x, u, t, σ=σ, β=β, ρ=ρ, γ=γ, cplmat=cplmat)
        dx[1] = σ * (x[2] - x[1])
        dx[2] = x[1] * (ρ - x[3]) - x[2]
        dx[3] = x[1] * x[2] - β * x[3]
        dx .*= γ
        dx .+= cplmat * map(ui -> ui(t), u.itp)   # Couple inputs
    end
    diffusion::DF = (dx, x, u, t, η=η) -> (dx .= η)
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = Inport(3)
    output::OP = Outport(3)
end

##### Pretty printing
show(io::IO, ds::SDESystem) = print(io,
    "SDESystem(drift:$(ds.drift), diffusion:$(ds.diffusion), readout:$(ds.readout), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::NoisyLorenzSystem) = print(io,
    "NoisyLorenzSystem(σ:$(ds.σ), β:$(ds.β), ρ:$(ds.ρ), η:$(ds.η), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ForcedNoisyLorenzSystem) = print(io,
    "ForcedNoisyLorenzSystem(σ:$(ds.σ), β:$(ds.β), ρ:$(ds.ρ), η:$(ds.η), γ:$(ds.γ), cplmat:$(ds.cplmat), ",
    "state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")


text:
"""
Import SO6 fileformat as a DataFrame

See EUROCONTROL NEST Manual for original SO6 fileformat description

SO6 DataFrame Column Names

#                   Field                   Type    Size    Comment
:SEGMENT_ID          segment identifier      String          first point name "_" last point name (see note)
:ADEP                origin of flight        String  4       ICAO code
:ADES                destination of flight   String  4       ICAO code
:ACTYPE              aircraft type           String  4
:DATETIMEBEGINSEGMENT                        DateTime
:DATETIMEENDSEGMENT                          DateTime
:FLBEGINSEGMENT      FL begin segment        Int64   1-3
:FLENDSEGMENT        FL begin segment        Int64   1-3
:STATUS              flight status           Int64   1       0=climb, 1=descent, 2=cruise
:CALLSIGN                                    String
:LATBEGINSEGMENT_DEG Latitude of begin of segment in degrees     Float64
:LONBEGINSEGMENT_DEG Longitude of begin of segment in degrees    Float64
:LATENDSEGMENT_DEG   Latitude of begin of segment in degrees     Float64
:LONENDSEGMENT_DEG   Longitude of begin of segment in degrees Float64
:FLIGHT_ID           Flight identifier       Int64           Unique ID
:SEQUENCE            Sequence number         Int64           Start at 1 for new flight
:SEGMENT_LENGTH_M    Segment length in m     Float64         In meters
:SEGMENT_PARITY      Segment parity          Int64           * see below

* Segment parity
0=NO
1=ODD
2=EVEN
3=ODD_LOW
4=EVEN_LOW
5=ODD_HIGH
6=EVEN_HIGH
7 =General Purpose
8=General Purpose
9=General Purpose

Example SO6 file format
EHAM_*AM80 EHAM UUEE A321 214900 214926 0 10 0 AFL2193 171214 171214 3138.483333 285.850000 3137.483333 284.183333 213765625 1 1.427867 0
*AM80_!AAEW EHAM UUEE A321 214926 215004 10 25 0 AFL2193 171214 171214 3137.483333 284.183333 3135.833333 284.033333 213765625 2 1.652550 0
!AAEW_!AAEX EHAM UUEE A321 215004 215026 25 35 0 AFL2193 171214 171214 3135.833333 284.033333 3134.166667 283.883333 213765625 3 1.669195 0
!AAEX_*AM46 EHAM UUEE A321 215026 215047 35 45 0 AFL2193 171214 171214 3134.166667 283.883333 3133.066667 283.783333 213765625 4 1.101704 0
"""
module SO6

export read

# include("utility.jl")
using ..util
using CSV
using Dates
using DataFrames

const fileformat = Dict(1=>String, 2=>String, 3=>String, 4=>String, 5=>String,
6=>String, 7=>Int64, 8=>Int64, 9=>String, 10=>String, 11=>String, 12=>String,
13=>Float64, 14=>Float64, 15=>Float64, 16=>Float64, 17=>Int64, 18=>Int64,
19=>Float64, 20=>Int64)

const header_format = ["SEGMENT_ID", "ADEP", "ADES", "ACTYPE",
"TIMEBEGINSEGMENT", "TIMEENDSEGMENT", "FLBEGINSEGMENT", "FLENDSEGMENT",
"STATUS", "CALLSIGN", "DATEBEGINSEGMENT", "DATEENDSEGMENT", "LATBEGINSEGMENT_DEG",
"LONBEGINSEGMENT_DEG", "LATENDSEGMENT_DEG", "LONENDSEGMENT_DEG", "FLIGHT_ID",
"SEQUENCE", "SEGMENT_LENGTH_M", "SEGMENT_PARITY"]

const yymmdd = DateFormat("YYmmdd")
const hhmmss = DateFormat("HHMMSS")
const yymmddhhmmss = DateFormat("YYmmddHHMMSS")
const year2000 = Year(2000)

function read(file)
    df = CSV.read(file, delim=" ", types=fileformat, header=header_format,
    copycols=true)
    reformat!(df);
    return df
end

function format_status(str::Union{String, Missing})
    if str === missing
        return missing
    end
    maybeint = tryparse(Int64, str)
    if maybeint == nothing
        return missing
    else
        return parse(Int64, str)
    end
end

function reformat!(df)
    #Date/Time conversion
    df[:,:DATETIMEBEGINSEGMENT] = format_datetime.(df[:,:DATEBEGINSEGMENT] .*
    df[:,:TIMEBEGINSEGMENT], yymmddhhmmss, addyear=year2000)

    df[:,:DATETIMEENDSEGMENT] = format_datetime.(df[:,:DATEENDSEGMENT] .*
    df[:,:TIMEENDSEGMENT], yymmddhhmmss, addyear=year2000)

    select!(df, Not(:TIMEBEGINSEGMENT))
    select!(df, Not(:TIMEENDSEGMENT))
    select!(df, Not(:DATEBEGINSEGMENT))
    select!(df, Not(:DATEENDSEGMENT))

    df[:,:TEMP] = format_status.(df[:,:STATUS])
    select!(df, Not(:STATUS))
    df[:,:STATUS] = df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:LATBEGINSEGMENT_DEG] = df[:,:LATBEGINSEGMENT_DEG] / 60.0
    df[:,:LONBEGINSEGMENT_DEG] = df[:,:LONBEGINSEGMENT_DEG] / 60.0
    df[:,:LATENDSEGMENT_DEG] = df[:,:LATENDSEGMENT_DEG] / 60.0
    df[:,:LONENDSEGMENT_DEG] = df[:,:LONENDSEGMENT_DEG] / 60.0
    df[:,:SEGMENT_LENGTH_M] = df[:,:SEGMENT_LENGTH_M] * 1852.0

    select!(df,[:SEGMENT_ID, :ADEP, :ADES, :ACTYPE, :DATETIMEBEGINSEGMENT,
    :DATETIMEENDSEGMENT, :FLBEGINSEGMENT, :FLENDSEGMENT, :STATUS, :CALLSIGN,
    :LATBEGINSEGMENT_DEG, :LONBEGINSEGMENT_DEG, :LATENDSEGMENT_DEG,
    :LONENDSEGMENT_DEG, :FLIGHT_ID, :SEQUENCE, :SEGMENT_LENGTH_M,
    :SEGMENT_PARITY])
end

end # module

ext: jl


text:
#===============
    STRUCTS
===============#

abstract type RTLight end

struct point_light <: RTLight
    position::RayTracer.rt_tuple
    intensity::RayTracer.rt_color
end

#===============
    FUNCTIONS
===============#


text:
################################################################################

using MindReader
using Test
using DelimitedFiles

################################################################################

tests = [

  # private functions

  # export functions
  "screening.jl",
  "writeCSV.jl",

]

################################################################################

@testset verbose = true "MindReader" begin

  for τ ∈ tests
    include(τ)
  end
end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DelimitedFiles

include("../../src/pecuzal_method.jl")
include("../../src/data_analysis_functions.jl")

Y_GA = readdlm("./scripts/Fooling systems/noise results/Y_GA.csv")
Y_mdop = readdlm("./scripts/Fooling systems/noise results/Y_mdop.csv")
Y_pec = readdlm("./scripts/Fooling systems/noise results/Y_pec.csv")
τ_vals_GA = readdlm("./scripts/Fooling systems/noise results/taus_GA.csv")
τ_vals_mdop = readdlm("./scripts/Fooling systems/noise results/taus_mdop.csv")
τ_vals_pec = readdlm("./scripts/Fooling systems/noise results/taus_pec.csv")
ts_vals_GA = readdlm("./scripts/Fooling systems/noise results/ts_GA.csv")
ts_vals_mdop = readdlm("./scripts/Fooling systems/noise results/ts_mdop.csv")
ts_vals_pec = readdlm("./scripts/Fooling systems/noise results/ts_pec.csv")


text:
#Vectors of uncertain values
uvals_x = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100]
uvals_y = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100];

# UncertainIndexDataset
UVX = UncertainValueDataset(uvals_x)
UVY = UncertainValueDataset(uvals_y)

constraints = [[TruncateStd(0.3) for x in 1:50]; [TruncateQuantiles(0.3, 0.7) for x in 1:50]]

@test resample(UVX, constraints) isa Vector{<:Real}
@test resample(UVY, constraints) isa Vector{<:Real}

# Test for some more complicated uncertain values
#----------------------------------------------
a = [UncertainValue(Normal, 1, 0.5) for i = 1:10]
b = UncertainValue(rand(1000))
c = UncertainValue(Uniform, rand(10000))
s = rand(-20:20, 100)
e = UncertainValue(s, rand(length(s)))
f = UncertainValue(2)
g = UncertainValue(3.1)

uvals = [a; b; c; e; f; g]
udata = UncertainValueDataset(uvals)

@test resample(udata, NoConstraint(), 10) isa Vector{Vector{Real}}
@test resample(udata, TruncateLowerQuantile(0.2), 10) isa Vector{Vector{Real}}
@test resample(udata, TruncateUpperQuantile(0.2), 10) isa Vector{Vector{Real}}
@test resample(udata, TruncateQuantiles(0.2, 0.8), 10) isa Vector{Vector{Real}}
@test resample(udata, TruncateMaximum(0.2), 10) isa Vector{Vector{Real}}
@test resample(udata, TruncateMinimum(0.2), 10) isa Vector{Vector{Real}}
@test resample(udata, TruncateRange(-40, 10), 10) isa Vector{Vector{Real}}
@test resample(udata, TruncateStd(1), 10) isa Vector{Vector{Real}}


text:
function Base.show(io::IO, diagram::Diagram)
    println(io, "already fixed: ", diagram.partial_sol)
    println(io, "root: ", only(diagram.layers[1]))
    for (l, variable) in enumerate(diagram.variables)
        println(io, "variable: ", variable)
        for pair in diagram.layers[l + 1]
            println(io, " ", pair)
        end
    end
end

function fixed_variables(diagram)
    return vcat(diagram.partial_sol, diagram.variables)
end

function add_transition(layer, new_state, new_node)
    if haskey(layer, new_state)
        if new_node.distance > layer[new_state].distance
            layer[new_state] = new_node
        end
    else
        layer[new_state] = new_node
    end
end

function build_layer(instance, diagram::Diagram{S,D,V}, variable) where {S,D,V}
    layer = Layer{S,D,V}()

    # Collect new states, keep only "best" arcs.
    for (state, node) in diagram.layers[end]
        for (arc, new_state) in transitions(instance, state, variable)
            new_node = Node{S,D,V}(node.distance + arc.value, arc)
            add_transition(layer, new_state, new_node)
        end
    end

    return layer
end

struct InOrder end

function next_variable(instance, diagram, variable_order::InOrder)
    num_variables = length(instance)
    fixed = fixed_variables(diagram)
    for i in 1:num_variables
        if !(i in fixed)  # TODO: efficient!
            return i
        end
    end
    return nothing
end

struct KeepAllNodes end

process(::KeepAllNodes, layer) = layer

"""
    top_down!(diagram::Diagram{S,D,V}, instance; variable_order, processing)

Build a decision diagram for given instance. It is developed layer by layer,
from the top, by following the state transitions. The first layer with a single
root node should already be present in the diagram.

Without any layer processing, the diagram will be exact, that is will contain a
path representing an optimal solution.

Restrictions or relaxations are achieved by dropping or merging the nodes and
states in each layer.
"""
function top_down!(
    diagram::Diagram{S,D,V}, instance;
    variable_order=InOrder(), processing=KeepAllNodes()
) where {S,D,V}
    @assert length(diagram.layers) == 1   # root layer

    # Intermediate layers
    while true
        variable = next_variable(instance, diagram, variable_order)
        if variable === nothing
            break
        end

        layer = build_layer(instance, diagram, variable)
        layer = process(processing, layer)   # restrict/relax

        push!(diagram.layers, layer)
        push!(diagram.variables, variable)
    end

    # Terminal node (last layer reduced to best)
    max_state, max_node = nothing, Node{S,D,V}(typemin(V))
    for (state, node) in diagram.layers[end]
        if node.distance > max_node.distance
            max_state = state
            max_node = node
        end
    end
    diagram.layers[end] = Layer{S,D,V}(Dict(max_state => max_node))
end

"""
    longest_path(diagram::Diagram{S,D,V})

Extract optimal solution by following a longest path from root to terminal.
"""
function longest_path(diagram::Diagram{S,D,V}) where {S,D,V}
    # Collect path in reverse, from terminal to root.
    terminal = only(values(diagram.layers[end].nodes))
    num_variables = length(diagram.partial_sol) + length(diagram.variables)
    decisions = Vector{D}(undef, num_variables)
    node, depth = terminal, length(diagram.layers) - 1
    while depth != 0
        decisions[diagram.variables[depth]] = node.inarc.decision
        state = node.inarc.tail
        node = diagram.layers[depth][state]
        depth -= 1
    end

    return Solution(decisions, terminal.distance)
end

function last_exact_layer(diagram)
    for (l, layer) in enumerate(diagram.layers)
        if !layer.exact
            # Current layer has at least one relaxed node.
            @assert l > 1

            # Return previous layer (all exact)
            return l - 1
        end
    end
    # If we reached the end then even the terminal layer is exact.
    return len(diagram.layers)
end

# TODO: pass solver object to store options, statistics, results?
"""
    branch_and_bound(instance; restrict, relax, variable_order)

Solve given instance using branch-and-bound. For each subproblem, a restriction
is computed to update the incumbent. Then a relaxation is computed. If the
problem can not be pruned from the relaxation bound, new subproblems are created
from an exact cutset in the relaxation's diagram.
"""
function branch_and_bound(instance; restrict, relax, variable_order=InOrder())
    state = initial_state(instance)
    S = typeof(state)
    D = domain_type(instance)
    V = value_type(instance)
    original_problem = Subproblem(Int[], D[], zero(V), state)

    # Assume maximization problems, so relaxation provides upper bound.
    # Following the book, we will use a node's distance (from root) as priority
    # and pick the smallest. This corresponds to picking nodes that are on early
    # layers?!
    #
    # TODO: turn this into a configurable strategy.
    problems = PriorityQueue(original_problem => zero(V))
    incumbent = Solution(D[], typemin(V))

    # Solve subproblems, one at a time.
    while !isempty(problems)
        current = dequeue!(problems)

        root_layer = Layer{S,D,V}(
            Dict(current.state => Node{S,D,V}(current.distance)))

        # solve restriction
        diagram = Diagram{S,D,V}(current.variables, [root_layer], [])
        top_down!(diagram, instance,
                  variable_order=variable_order, processing=restrict)
        solution = longest_path(diagram)

        # update incumbent
        if solution.objective > incumbent.objective
            for (variable, decision) in zip(current.variables, current.decisions)
                solution.decisions[variable] = decision
            end
            incumbent = solution
        end

        # have we solved the subproblem already?
        if all(l -> l.exact, diagram.layers)
            continue
        end

        # solve relaxation
        diagram = Diagram{S,D,V}(current.variables, [root_layer], [])
        top_down!(diagram, instance;
                  variable_order=variable_order, processing=relax)
        solution = longest_path(diagram)

        # create subproblems if not pruned
        if solution.objective > incumbent.objective
            cutset = last_exact_layer(diagram)
            @assert length(diagram.layers[cutset]) > 1
            for (sub_state, sub_node) in diagram.layers[cutset]
                depth = cutset - 1
                new_decisions = Vector{D}(undef, depth)
                node = sub_node
                while depth != 0
                    new_decisions[depth] = node.inarc.decision
                    state = node.inarc.tail
                    node = diagram.layers[depth][state]
                    depth -= 1
                end

                variables = vcat(current.variables,
                                 diagram.variables[1:cutset - 1])
                decisions = vcat(current.decisions, new_decisions)

                subproblem = Subproblem(variables, decisions,
                                        sub_node.distance, sub_state)
                problems[subproblem] = subproblem.distance
            end
        end
    end

    return incumbent
end


text:
import Base.isempty
mutable struct ElectrodermalData{T₁ <: Number} <: AbstractModel
    timestamp::Float64
    sampling_frequency::Int
    unprocessed_eda::Vector{T₁}
    eda::Vector{T₁}
end

function get_timestamp(e::ElectrodermalData)
    e.timestamp
end

function set_timestamp!(e::ElectrodermalData, timestamp::Float64)
    e.timestamp = timestamp
end

function get_sampling_frequency(e::ElectrodermalData)
    e.sampling_frequency
end

function set_sampling_frequency!(e::ElectrodermalData, sampling_frequency::Int)
    e.sampling_frequency = sampling_frequency
end

function get_unprocessed_eda(e::ElectrodermalData)
    e.unprocessed_eda
end

function set_unprocessed_eda!(e::ElectrodermalData, eda::AbstractVector)
    e.unprocessed_eda = eda
end

function get_eda(e::ElectrodermalData)
    e.eda
end

function set_eda!(e::ElectrodermalData, eda::AbstractVector)
    e.eda = eda
end

function isempty(e::ElectrodermalData)
    isempty(get_unprocessed_eda(e))
end


text:
module SubgridScaleParameters

using CLIMA.ParametersType

  @exportparameter C_smag         0.15                    "Standard Smagorinsky Coefficient"
  @exportparameter inv_Pr_turb    3                       "Turbulent Prandtl Number"
  @exportparameter Prandtl_air    71//100                 "Molecular Prandtl Number, dry air"
  @exportparameter c_a_KASM       0.10                    "cₐ KASM (2006)"
  @exportparameter c_e1_KASM      0.19                    "cₑ₁ KASM (2006)"
  @exportparameter c_e2_KASM      0.51                    "cₑ₂ KASM (2006)"
  @exportparameter c_1_KASM       c_a_KASM*0.76^2         "c₁  KASM (2006)"
  @exportparameter c_2_KASM       c_e2_KASM+2*c_1_KASM    "c₂ KASM (2006)"
  @exportparameter c_3_KASM       c_a_KASM^(3/2)          "c₃ KASM (2006)"

  end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
function load_bound_states(directory::AbstractString, k::Integer)
    dbl_pat = r"([ -])([0-9]).([0-9]+)D([+-])([0-9]{2})" # Double-precision in Fortran: 0.0000D+00
    flt_pat = r"([ -])([0-9]).([0-9]+)E([+-])([0-9]{2})"
    open(joinpath(directory, @sprintf("bound.%03d", k))) do file
        ns,kch,kcp,nhm,nbound,lpar,ispar,ipar =
            map(s -> parse(Float64, s),
                split(split(readline(file), "=>")[1]))
        map(1:nbound) do i
            label = split(readline(file))[2]
            E_au,E_bind,n_eff = map(s -> parse(Float64, s),
                                    split(split(readline(file), "=>")[1]))
            c = Vector{Float64}()
            while length(c) < nhm
                line = readline(file)
                ms = matchall(dbl_pat, line)
                t = map(ms) do m
                    parse(Float64, replace(m, "D", "E"))
                end
                append!(c, t)
            end
            Dict(:label => label,
                 :E_au => E_au, :E_bind => E_bind, :n_eff => n_eff,
                 :c => c)
        end
    end
end

function asmatrix(v)
    V = Matrix{Float64}(length(v[1][:c]), length(v))
    for j = 1:length(v)
        V[:,j] = v[j][:c]
    end
    V
end

function eigenbasis(directory::AbstractString, k::Integer, Ip::Real)
    v = load_bound_states(directory, k)
    V = asmatrix(v)
    E = Float64[v[j][:E_bind]/27.211-Ip for j = eachindex(v)]
    E,V
end

function eigenbasis(directory::AbstractString,
                    krange::UnitRange{Int}, Ip::Real)
    @time E,V = zip(map(k -> eigenbasis(directory, k, Ip), krange)...)
    E = Vector{Float64}[EE for EE in E]
    V = Matrix{Float64}[V[k] for k=1:length(V)]
    E,V
end

export load_bound_states, eigenbasis


text:
# Single Layer
function (SL::SingleLayerKernel{T,S})(x::SVector,y::SVector)::T  where {T,S<:Helmholtz}
    N = ambient_dimension(S)
    k = SL.op.k
    r = x - y
    d = norm(r)
    d == 0 && return zero(T)
    if N==2
        return im/4 * hankelh1(0,k*d)
    elseif N==3
        return 1/(4π)/d * exp(im*k*d)
    end
end

# Double Layer Kernel
function (DL::DoubleLayerKernel{T,S})(x::SVector,y::SVector,ny::SVector)::T where {T,S<:Helmholtz}
    N = ambient_dimension(S)
    k = DL.op.k
    r = x - y
    d = norm(r)
    d == 0 && return zero(T)
    if N==2
        return im*k/4/d * hankelh1(1,k*d) .* dot(r,ny)
    elseif N==3
        return 1/(4π)/d^2 * exp(im*k*d) * ( -im*k + 1/d ) * dot(r,ny)
    end
end

# Adjoint double Layer Kernel
function (ADL::AdjointDoubleLayerKernel{T,S})(x::SVector,y::SVector,nx::SVector)::T where {T,S<:Helmholtz}
    N = ambient_dimension(S)
    k = ADL.op.k
    r = x - y
    d = norm(r)
    d == 0 && return zero(T)
    if N==2
        return -im*k/4/d * hankelh1(1,k*d) .* dot(r,nx)
    elseif N==3
        return -1/(4π)/d^2 * exp(im*k*d) * ( -im*k + 1/d ) * dot(r,nx)
    end
end

# Hypersingular kernel
function (HS::HyperSingularKernel{T,S})(x::SVector,y::SVector,nx::SVector,ny::SVector)::T where {T,S<:Helmholtz}
    N = ambient_dimension(S)
    k = HS.op.k
    r = x - y
    d = norm(r)
    d == 0 && return zero(T)
    if N==2
        RRT = r*transpose(r) # r ⊗ rᵗ
        # TODO: rewrite the operation below in a more clear/efficient way
        return transpose(nx)*((-im*k^2/4/d^2*hankelh1(2,k*d)*RRT + im*k/4/d*hankelh1(1,k*d)*I)*ny)
    elseif N==3
        RRT   = r*transpose(r) # r ⊗ rᵗ
        term1 = 1/(4π)/d^2 * exp(im*k*d) * ( -im*k + 1/d ) * I
        term2 = RRT/d * exp(im*k*d)/(4*π*d^4) * (3*(d*im*k-1) + d^2*k^2)
        return  transpose(nx)*(term1 + term2) * ny
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# Common facilities

typealias AttributeDict Dict{String, Any}

#################################################
#
#  vertex types
#
################################################

immutable KeyVertex{K}
    index::Int
    key::K
end

KeyVertex{K}(idx::Int, key::K) = KeyVertex{K}(idx, key)
make_vertex{V<:KeyVertex}(g::AbstractGraph{V}, key) = V(num_vertices(g) + 1, key)
vertex_index(v::KeyVertex) = v.index

type ExVertex
    index::Int
    label::Compat.UTF8String
    attributes::AttributeDict

    ExVertex(i::Int, label::AbstractString) = new(i, label, AttributeDict())
end

make_vertex(g::AbstractGraph{ExVertex}, label::AbstractString) = ExVertex(num_vertices(g) + 1, String(label))
vertex_index(v::ExVertex) = v.index
attributes(v::ExVertex, g::AbstractGraph) = v.attributes

typealias ProvidedVertexType @compat(Union{KeyVertex, ExVertex})

vertex_index{V<:ProvidedVertexType}(v::V, g::AbstractGraph{V}) = vertex_index(v)

function vertex_index{V}(v::V, g::AbstractGraph{V})
    @graph_requires g vertex_list
    if applicable(vertex_index, v)
        return vertex_index(v)
    end
    return vertex_index(v, vertices(g)) # slow linear search
end

vertex_index(v, vs::AbstractArray) = findfirst(vs, v)

#################################################
#
#  edge types
#
################################################

immutable Edge{V}
    index::Int
    source::V
    target::V
end
typealias IEdge Edge{Int}

Edge{V}(i::Int, s::V, t::V) = Edge{V}(i, s, t)
make_edge{V,E<:Edge}(g::AbstractGraph{V,E}, s::V, t::V) = Edge(num_edges(g) + 1, s, t)

revedge{V}(e::Edge{V}) = Edge(e.index, e.target, e.source)

edge_index(e::Edge) = e.index
source(e::Edge) = e.source
target(e::Edge) = e.target
source{V}(e::Edge{V}, g::AbstractGraph{V}) = e.source
target{V}(e::Edge{V}, g::AbstractGraph{V}) = e.target

type ExEdge{V}
    index::Int
    source::V
    target::V
    attributes::AttributeDict
end

=={V}(e1::ExEdge{V}, e2::ExEdge{V}) = (e1.index == e2.index &&
                                       e1.source == e2.source &&
                                       e1.target == e2.target)

ExEdge{V}(i::Int, s::V, t::V) = ExEdge{V}(i, s, t, AttributeDict())
ExEdge{V}(i::Int, s::V, t::V, attrs::AttributeDict) = ExEdge{V}(i, s, t, attrs)
make_edge{V}(g::AbstractGraph{V}, s::V, t::V) = ExEdge(num_edges(g) + 1, s, t)

revedge{V}(e::ExEdge{V}) = ExEdge{V}(e.index, e.target, e.source, e.attributes)

edge_index(e::ExEdge) = e.index
source(e::ExEdge) = e.source
target(e::ExEdge) = e.target
source{V}(e::ExEdge{V}, g::AbstractGraph{V}) = e.source
target{V}(e::ExEdge{V}, g::AbstractGraph{V}) = e.target
attributes(e::ExEdge, g::AbstractGraph) = e.attributes

#################################################
#
#  iteration
#
################################################

# general reindexed vector

immutable ReindexedVec{T, Vec<:AbstractVector, I<:AbstractVector{Int}}
    src::Vec
    inds::I
end

ReindexedVec{T}(a::AbstractVector{T}, inds::AbstractVector{Int}) =
    ReindexedVec{T,typeof(a),typeof(inds)}(a, inds)

length(a::ReindexedVec) = length(a.inds)
isempty(a::ReindexedVec) = isempty(a.inds)
getindex(a::ReindexedVec, i::Integer) = a.src[a.inds[i]]

start(a::ReindexedVec) = start(a.inds)
done(a::ReindexedVec, s) = done(a.inds, s)
next(a::ReindexedVec, s) = ((i, s) = next(a.inds); (a.src[i], s))

# iterating over targets

immutable TargetIterator{G<:AbstractGraph,EList}
    g::G
    lst::EList
end

TargetIterator{G<:AbstractGraph,EList}(g::G, lst::EList) =
    TargetIterator{G,EList}(g, lst)

length(a::TargetIterator) = length(a.lst)
isempty(a::TargetIterator) = isempty(a.lst)
getindex(a::TargetIterator, i::Integer) = target(a.lst[i], a.g)

start(a::TargetIterator) = start(a.lst)
done(a::TargetIterator, s) = done(a.lst, s)
next(a::TargetIterator, s::Int) = ((e, s) = next(a.lst, s); (target(e, a.g), s))

# iterating over sources

immutable SourceIterator{G<:AbstractGraph,EList}
    g::G
    lst::EList
end

SourceIterator{G<:AbstractGraph,EList}(g::G, lst::EList) =
    SourceIterator{G,EList}(g, lst)

length(a::SourceIterator) = length(a.lst)
isempty(a::SourceIterator) = isempty(a.lst)
getindex(a::SourceIterator, i::Integer) = source(a.lst[i], a.g)

start(a::SourceIterator) = start(a.lst)
done(a::SourceIterator, s) = done(a.lst, s)
next(a::SourceIterator, s::Int) = ((e, s) = next(a.lst, s); (source(e, a.g), s))

#################################################
#
#  Edge Length Visitors
#
################################################

abstract AbstractEdgePropertyInspector{T}

edge_property_requirement{T, V}(visitor::AbstractEdgePropertyInspector{T}, g::AbstractGraph{V}) = nothing

type ConstantEdgePropertyInspector{T} <: AbstractEdgePropertyInspector{T}
  value::T
end

edge_property{T}(visitor::ConstantEdgePropertyInspector{T}, e, g) = visitor.value

type VectorEdgePropertyInspector{T} <: AbstractEdgePropertyInspector{T}
  values::Vector{T}
end

edge_property{T,V}(visitor::VectorEdgePropertyInspector{T}, e, g::AbstractGraph{V}) = visitor.values[edge_index(e, g)]

edge_property_requirement{T, V}(visitor::AbstractEdgePropertyInspector{T}, g::AbstractGraph{V}) = @graph_requires g edge_map

type AttributeEdgePropertyInspector{T} <: AbstractEdgePropertyInspector{T}
  attribute::Compat.UTF8String
end

function edge_property{T}(visitor::AttributeEdgePropertyInspector{T},edge::ExEdge, g)
    convert(T,edge.attributes[visitor.attribute])
end
#################################################
#
#  convenient functions
#
################################################

isnz(x::Bool) = x
isnz(x::Number) = x != zero(x)

intrange(n::Integer) = 1:convert(Int,n)

multivecs{T}(::Type{T}, n::Int) = [T[] for _ =1:n]

function collect_edges{V,E}(graph::AbstractGraph{V,E})
    local edge_list
    if implements_edge_list(graph)
        collect(edges(graph))

    elseif implements_vertex_list(graph) && implements_incidence_list(graph)
        edge_list = Array(E, 0)
        sizehint!(edge_list, num_edges(graph))

        for v in vertices(graph)
            for e in out_edges(v, graph)
                push!(edge_list, e)
            end
        end
        edge_list
    else
        throw(ArgumentError("graph must implement either edge_list or incidence_list."))
    end
end

immutable WeightedEdge{E,W}
    edge::E
    weight::W
end

isless{E,W}(a::WeightedEdge{E,W}, b::WeightedEdge{E,W}) = a.weight < b.weight

function collect_weighted_edges{V,E,W}(graph::AbstractGraph{V,E}, weights::AbstractEdgePropertyInspector{W})

    edge_property_requirement(weights, graph)

    wedges = Array(WeightedEdge{E,W}, 0)
    sizehint!(wedges, num_edges(graph))

    if implements_edge_list(graph)
        for e in edges(graph)
            w = edge_property(weights, e, graph)
            push!(wedges, WeightedEdge(e, w))
        end

    elseif implements_vertex_list(graph) && implements_incidence_list(graph)
        for v in vertices(graph)
            for e in out_edges(v, graph)
                w = edge_property(weights, e, graph)
                push!(wedges, WeightedEdge(e, w))
            end
        end
    else
        throw(ArgumentError("graph must implement either edge_list or incidence_list."))
    end

    return wedges
end

function collect_weighted_edges{V,E,W}(graph::AbstractGraph{V,E}, weights::AbstractVector{W})
    visitor::AbstractEdgePropertyInspector{D} = VectorEdgePropertyInspector(edge_dists)
    collect_weighted_edges(graph, visitor)
end


text:
include("lobpcg_sep.jl")

mutable struct Ham
    Ns::Int64
    Ls::Float64
    kmul::Array{Float64,2} # wait until everything is clearer
    dx::Float64
    gridpos::Array{Float64,2}
    posstart::Int64
    posidx::Int64
    H                  # no idea what is this
    rhoa::Array{Float64,2}               # pseudo-charge for atoms (negative)
    rho::Array{Float64,2}                # electron density
    Vhar::Array{Float64,2}               # Hartree potential for both electron and nuclei
    Vtot::Array{Float64,2}               # total potential
    drhoa::Array{Float64,2}              # derivative of the pseudo-charge (on each atom)
    ev::Array{Float64,1}
    psi::Array{Float64,2}
    fermi::Float64
    occ
    nspin::Int64
    Neigs::Int64    # QUESTION: Number of eigenvalues?
    atoms::Atoms
    Eband::Float64              # Total band energy
    Fband::Float64               # Helmholtz band energy
    Ftot::Float64                # Total Helmholtz energy
    YukawaK::Float64             # shift for the potential
    epsil0::Float64
    Tbeta::Float64               # temperature 1beta = 1/T

    function Ham(Lat, Nunit, n_extra, dx, atoms,YukawaK, epsil0,Tbeta)
        Ls = Nunit * Lat;
        Ns = round(Integer, Ls / dx);

        #
        dx = Ls / Ns;
        # defining the grid
        gridpos = zeros(Ns,1) # allocating as a 2D Array
        gridpos[:,1] = transpose(collect(0:Ns-1)).*dx;
        posstart = 0;
        posidx   = 0;

        # Initialize the atom positions
        atoms   = atoms;
        Neigs = sum(atoms.nocc) + n_extra;

        Ls_glb = Ls;
        Ns_glb = Ns;

        # we define the Fourier multipliers as an 2D array
        kx = zeros(Ns,1);
        kx[:,1] = vcat( collect(0:Ns/2-1), collect( -Ns/2:-1) )* 2 * pi / Ls;
        kmul = kx.^2/2;

        rhoa, drhoa = pseudocharge(gridpos, Ls_glb, atoms,YukawaK,epsil0);

        # TODO: we need to figure out the type of each of the fields to properlu
        # initialize them
        H  = []
        rho = zeros(1,1);
        Vhar = zeros(1,1);
        Vtot = zeros(1,1);
        # drhoa = []
        ev = []
        psi = zeros(1,1);
        fermi = 0.0;
        occ = []
        Eband = 0.0
        Fband = 0.0
        Ftot = 0.0
        nspin = 1;

        new(Ns, Ls, kmul, dx, gridpos, posstart, posidx, H, rhoa,
            rho, Vhar, Vtot, drhoa, ev, psi, fermi, occ,nspin, Neigs, atoms,
            Eband, Fband, Ftot, YukawaK, epsil0, Tbeta)
    end
end

# importing the necessary functions to comply with Julia duck typing
using LinearAlgebra
import Base.*
import LinearAlgebra.mul!
import Base.eltype
import Base.size
import LinearAlgebra.issymmetric

function *(H::Ham, x::Array{Float64,1})
    # matvec overloading
    # Function  to  overload the application of the Hamiltonian times avector
    y_lap  = lap(H,x);
    y_vtot = Vtot(H,x);

    return y_lap + y_vtot;
end

# TODO: this should be optimized in order to take advantage of BLAS 3 operations
function *(H::Ham, X::Array{Float64,2})
    # Function  to  overload the application of the Hamiltonian times a matrix
    Y = zeros(size(X));
    mul!(Y, H, X)
    return Y
    # # new version that acts on the full matrices
    # Y_lap  = lap(H,X);
    # Y_vtot = Vtot(H,X);
    #
    # return Y_lap + Y_vtot;
end

function LinearAlgebra.mul!(Y::Array{Float64,2}, H::Ham, V::Array{Float64,2})
    # in place matrix matrix multiplication
    @assert(size(Y) == size(V))
    Y[:,:] = lap(H,V)
    Y += Vtot(H,V)

    return Y
end

# optimized version for eigs (it uses sub arrays to bypass the inference step)
function LinearAlgebra.mul!(Y::SubArray{Float64,1,Array{Float64,1}},
                  H::Ham,
                  V::SubArray{Float64,1,Array{Float64,1}})
    # in place matrix matrix multiplication
    @assert(size(Y) == size(V))
    for ii = 1:size(V,2)
        Y[:,ii] = H*V[:,ii]
    end

end

function size(H::Ham)
    return (H.Ns, H.Ns)
end

function eltype(H::Ham)
    # we work always in real space, the Fourier operations are on ly meant as
    # a pseudo spectral discretization
    return typeof(1.0)
end

function LinearAlgebra.issymmetric(H::Ham)
    return true
end

function update_psi!(H::Ham, eigOpts::eigOptions)
    # we need to add some options to the update
    # functio to solve the eigenvalue problem for a given rho and Vtot

    if eigOpts.eigmethod == "eigs"
        # TODO: make sure that eigs works with overloaded operators
        # TODO: take a look a the interface of eigs in Julia
        (ev,psi,nconv,niter,nmult,resid) = eigs(H,   # Hamiltonian
                                                nev=H.Neigs+4, # number of eigs
                                                which=:SR, # small real part
                                                ritzvec=true, # provide Ritz v
                                                tol=eigOpts.eigstol, # tolerance
                                                maxiter=eigOpts.eigsiter) #maxiter
        #assert(flag == 0);
        # TODO: this has a bug somewhere !!!! fix me!!!!
    elseif  eigOpts.eigmethod == "lobpcg_sep"
        # not working... to be fixed
        X0 = qr(rand(H.Ns, H.Neigs)).Q[:,1:H.Neigs]
        precond(x) = inv_lap(H,x)

        (ev,psi, iter) = lobpcg_sep(H, X0, precond, H.Neigs,
                                    tol=eigOpts.eigstol,
                                    maxiter=eigOpts.eigsiter)

    elseif  eigOpts.eigmethod == "eig"
        # we use a dense diagonalization
        A = create_Hamiltonian(H)
        # checkign that A is symetric
        @assert issymmetric(A)
        (ev, psi) = (eigen(A)...,)

    end

    # sorting the eigenvalues, eigs already providesd them within a vector
    ind = sortperm(ev)[1:H.Neigs];
    # updating the eigenvalues
    H.ev = ev[ind]
    # updating the eigenvectors
    H.psi = psi[:, ind];
end

function update_rho!(H::Ham, nocc::Int64)

    ev = H.ev;
    (occ, fermi) = get_occ(ev, nocc, H.Tbeta);
    occ = occ * H.nspin;
    rho = sum(H.psi.^2*Matrix(Diagonal(occ)),dims=2)/H.dx;

    # Total energy
    E = sum(ev.*occ);

    # Helmholtz free energy
    intg = H.Tbeta*(fermi .- ev);
    ff = zeros(H.Neigs,1);
    for i = 1 : H.Neigs
      if( intg[i] > 30 )  # Avoid numerical problem.
        ff[i] = ev[i]-fermi;
      else
        ff[i] = -1/H.Tbeta * log(1+exp(intg[i]));
      end
    end
    F = sum(ff.*occ) + fermi * nocc * H.nspin;

    H.occ = occ;
    H.fermi = fermi;
    H.Eband = E;
    H.Fband = F;
    H.rho = rho;
end

function update_rhoa!(H::Ham)
    H.rhoa, H.drhoa = pseudocharge(H.gridpos, H.Ls, H.atoms,H.YukawaK,H.epsil0);
end

function create_Hamiltonian(H::Ham)
    # create the matrix version of the Hmailtonian
      A = real(ifft(diagm(0 => H.kmul[:])*fft(Matrix{Float64}(I, length(H.kmul), length(H.kmul)),1),1));
      A += diagm(0 => H.Vtot[:])
      # we symmetrize A
    return 0.5*(A + A')
end

function lap(H::Ham,x::Array{Float64,1})
    # we ask for a vector, given that we will consider the vector to be
    # a nx1 matrix
    # TODO: this can be optimized using rfft
    # TODO: we can surther accelerate this using a in-place multiplication

    ytemp = H.kmul.*fft(x);
    return real(ifft(ytemp))
end

function lap(H::Ham, x::Array{Float64,2})
    # application of the laplacian part to a matrix.
    # TODO: this can be optimized using rfft
    # TODO: we can surther accelerate this using a in-place multiplication

    ytemp = broadcast(*,H.kmul,fft(x,1));
    return real(ifft(ytemp,1))
end

function inv_lap(H::Ham, x::Array{Float64,1})
    # inverse laplacian, to be used as a preconditioner for the
    # lobpcg algorithm

    inv_kmul = zeros(size(H.kmul))
    inv_kmul[1] = 0;
    inv_kmul[2:end] = 1 ./H.kmul[2:end];

    ytemp = inv_kmul.*fft(x);
    return real(ifft(ytemp))
end

function prec(H::Ham, x::Array{Float64,1})
    # preconditioner for lobpcg
## matlab code to translate
#     for ik = 1:nkpts
#     X  = gkincell{ik};
#     Y  = 27.0 + X.*(18.0 + X.*(12.0 + 8.0*X));
#     p{ik}  = Y./(Y + 16.0*X.^4);
# end;

    # inv_kmul = zeros(size(H.kmul))
    # inv_kmul[1] = 0;
    # inv_kmul[2:end] = 1./H.kmul[2:end];

    inv_kmul = zeros(size(H.kmul))
    inv_kmul[1] = 0;
    inv_kmul[2:end] = H.kmul[2:end];

    Y = 27.0 .+ inv_kmul.*(18.0 .+ inv_kmul.*(12.0 .+ 8.0 .*inv_kmul));
    yfft = Y./(Y + 16.0 .*inv_kmul.^4);
    ytemp = yfft.*fft(x);
    return real(ifft(ytemp))
end

function Vtot(H::Ham,x::Array{Float64,1})
    # application of the potential part of the Hamiltonian to a vector
    return (H.Vtot).*x
end

function Vtot(H::Ham, X::Array{Float64,2})
    # application of the potential part of the Hamiltonian to a matric
    return broadcast(*,H.Vtot, X)
end

function hartree_pot_bc(rho, H::Ham)
# computes the Hartree potential and energy in periodic boundary potential
# by solving the 1-d Yakuwa equation.
    # we call the (hopefully compiler optmized version)
    return hartree_pot_bc(rho, H.Ls, H.YukawaK, H.epsil0)
end

function update_vtot!(H::Ham, mixOpts)
    # TODO here we need to implement the andersson mix
    # I added the signature

    (Vtotnew,Verr) = update_pot!(H)

    # TODO: add a swtich to use different kinds of mixing here
    betamix = mixOpts.betamix;
    mixdim = mixOpts.mixdim;
    ymat = mixOpts.ymat;
    smat = mixOpts.smat;
    iter = mixOpts.iter;

    (Vtotmix,ymat,smat) = anderson_mix(H.Vtot,Vtotnew,
        betamix, ymat, smat, iter, mixdim);

    # they are already modified inside the function
    # mixOpts.ymat = ymat;
    # mixOpts.smat = smat;
    mixOpts.iter += 1;

    # updating total potential
    H.Vtot = Vtotmix;
    return Verr
end

function update_vtot!(H::Ham, mixOpts::kerkerMixOptions)
    # TODO here we need to implement the andersson mix
    # I added the signature

    (Vtotnew,Verr) = update_pot!(H)

    # println("using kerker mixing")
    # computign the residual
    res     = H.Vtot - Vtotnew;

    # println("appliying the preconditioner")
    resprec = kerker_mix(res, mixOpts.KerkerB, mixOpts.kx,
                              mixOpts.YukawaK, mixOpts.epsil0);

    # println("performing the linear mixing ")
    Vtotmix =  H.Vtot - mixOpts.betamix * resprec;

    # updating total potential
    # println("saving the total potentail ")
    H.Vtot = Vtotmix;

    # println("returning the error")
    return Verr
end

function update_vtot!(H::Ham, mixOpts::andersonPrecMixOptions)
    # TODO here we need to implement the andersson mix
    # I added the signature

    (Vtotnew,Verr) = update_pot!(H)

    betamix = mixOpts.betamix;
    mixdim = mixOpts.mixdim;
    ymat = mixOpts.ymat;
    smat = mixOpts.smat;
    iter = mixOpts.iter;

    # TODO change this so this is properly modified inside the function
    (Vtotmix,ymat,smat) = prec_anderson_mix(H.Vtot,Vtotnew,
        betamix, ymat, smat, iter, mixdim, mixOpts.prec, mixOpts.precargs)

    # they are already modified inside the function
    # just to be sure we can leave them here
    mixOpts.ymat = ymat;
    mixOpts.smat = smat;
    mixOpts.iter += 1;

    # updating total potential
    H.Vtot = Vtotmix;
    return Verr
end

import Statistics.mean
using FFTW

function init_pot!(H::Ham, nocc::Int64)
    # nocc number of occupied states
    #function to initialize the potential in the Hamiltonian class
    rho  = -H.rhoa;
    rho  = rho / ( sum(rho)*H.dx) * (nocc*H.nspin);
    H.rho = rho;
    H.Vhar = hartree_pot_bc(H.rho+H.rhoa, H);
    H.Vtot = H.Vhar;   # No exchange-correlation
    H.Vtot = H.Vtot .- mean(H.Vtot);# IMPORTANT (zero mean?)
end

function update_pot!(H::Ham)
     # TODO: I dont' know in how many different ways this si updateded
     # computing the hartree potenatial
    H.Vhar = hartree_pot_bc(H.rho+H.rhoa,H);
    # here Vtotnew only considers the
    Vtotnew  = H.Vhar;  # no exchange-correlation so far
    Verr = norm(Vtotnew-H.Vtot)./norm(H.Vtot); # computing the relative error

    # NOTE: H.Fband is only the band energy here.  The real total energy
    # is calculated using the formula below:
    H.Ftot = H.Fband + 1/2 * sum((H.rhoa-H.rho).*H.Vhar)*H.dx;
    return (Vtotnew,Verr) # returns the differnece betwen two consecutive iterations
end

using Arpack

# TODO: add a default setting for the scfOpts
function scf!(H::Ham, scfOpts::scfOptions)

    # vector containing the hostorical of the results
    VtoterrHist = zeros(scfOpts.scfiter)
    eigOpts = scfOpts.eigOpts;
    mixOpts = scfOpts.mixOpts;

    # number of occupied states
    Nocc = round(Integer, sum(H.atoms.nocc) / H.nspin);

    # we test first updating the psi

    for ii = 1:scfOpts.scfiter
        # solving the linear eigenvalues problem
        update_psi!(H, eigOpts);

        # update the electron density
        update_rho!(H,Nocc);

        # update the total potential, and compute the
        # differnce between the potentials
        Verr = update_vtot!(H, mixOpts);

        # save the error
        VtoterrHist[ii] = Verr ;
        # test if the problem had already satiesfied the tolerance
        if scfOpts.SCFtol > Verr
            break
        end
    end

    return VtoterrHist[VtoterrHist.>0]
end

function lap_opt(H::Ham,x::Array{Float64,1})
    # we ask for a 2 vector, given that we will consider the vector to be
    xFourier = rfft(x)
    laplacian_fourier_mult!(xFourier, H.Ls)
    return irfft(xFourier, H.Ns )
end

function lap_opt!(H::Ham,x::Array{Float64,1})
    # we ask for a 2 vector, given that we will consider the vector to be
    xFourier = rfft(x)
    laplacian_fourier_mult!(xFourier, H.Ls)
    x[:] = irfft(xFourier, H.Ns )
end

function laplacian_fourier_mult!(R::Vector{ComplexF64}, Ls::Float64 )
    c = (2 * pi / Ls)^2/2
    @simd for ii = 1:length(R)
        @inbounds R[ii] = (ii-1)^2*c*R[ii]
    end
end

# Function to compute the force and update ham.atoms.force
function get_force!(H::Ham)
    atoms  = H.atoms
    Natoms = atoms.Natoms
    rhotot = H.rho + H.rhoa
    for i=1:Natoms
        # IMPORTANT: derivative is taken w.r.t atom positions, which introduces the minus sign
        dV = -hartree_pot_bc_opt_vec(H.drhoa[:,i], H.Ls, H.YukawaK, H.epsil0)
        # Force is negative gradient
        atoms.force[i] = - sum(dV.*rhotot)*H.dx
    end
end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
export ZeroPad
using Knet.Ops21: zeropad

mutable struct ZeroPad; padding; xsize; y; end

ZeroPad(p) = ZeroPad(p, nothing, nothing)

function (z::ZeroPad)(x)
    if size(x) !== z.xsize || typeof(x) !== typeof(z.y)
        z.y = zeropad(x, z.padding)
        z.xsize = size(x)
    else
        z.y = zeropad(x, z.padding; y=z.y, fillzero=false)
    end
    return z.y
end


text:
function gaussianpulse(t::Real,p::Vector)
    # Gaussian pulse centered on t=0, area normalized, and starting and finishing at y=0
    σ   = p[1] # gaussian standard deviation (FWHM = 2√(2ln2)σ ≈ 2.355σ)
    tg  = p[2] # total gate time
    ω10 = p[3] # qubit 0-1 transition angular freq
    ϕ   = p[4] # phase
    A   = p[5] # amplitude, (π for π-pulse, etc...)
    B = inv(√(2π)*σ*erf(tg/(√(8)*σ))-tg*gaussian(0.5tg,σ))
    Ɛπ = A*B*(gaussian(t,σ)-gaussian(0.5tg,σ))
    return Ɛπ*cos(ω10*t+ϕ)
end

function cosinepulse(t::Real,p::Vector)
    # Cosine pulse centered on t=0
    σ   = p[1] # not used
    tg  = p[2] # total gate time
    ω10 = p[3] # qubit 0-1 transition angular freq
    ϕ   = p[4] # phase
    A   = p[5] # amplitude, (π for π-pulse, etc...)
    Ɛπ = A/tg*(cos(2π*t/tg)+1)
    return Ɛπ*cos(ω10*t+ϕ)
end

function dragpulse(t::Real,p::Vector)
    # 2nd order DRAG-enhanced Gaussian pulse (Y-only correction) centered on t=0
    # 10.1103/PhysRevA.83.012308
    σ   = p[1] # gaussian standard deviation (FWHM = 2√(2ln2)σ ≈ 2.355σ)
    tg  = p[2] # total gate time
    ω10 = p[3] # qubit 0-1 transition angular freq
    Δ   = p[4] # ω21 - ω10
    λ   = p[5] # drive 1-2 transition amplitude ~√(2)
    ϕ   = p[6] # phase
    A   = p[7] # amplitude, (π for π-pulse, etc...)
    B = inv(√(2π)*σ*erf(tg/(√(8)*σ))-tg*gaussian(0.5tg,σ))
    Ɛπ = A*B*(gaussian(t,σ)-gaussian(0.5tg,σ))
    Ɛπ′ = A*B*gaussianprime(t,σ)
    ωt = ω10*t
    return (Ɛπ - λ^2*(λ^2-4)*Ɛπ^3/32Δ^2)*cos(ωt+ϕ) + (-λ^2*Ɛπ′/4Δ)*sin(ωt+ϕ)
end

function dragpulse_det(t::Real,p::Vector)
    # 2nd order DRAG-enhanced Gaussian pulse (Y-only) centered on t=0
    # with detuning
    # PhysRevLett.116.020501
    σ   = p[1] # gaussian standard deviation (FWHM = 2√(2ln2)σ ≈ 2.355σ)
    tg  = p[2] # total gate time
    ω10 = p[3] # qubit 0-1 transition angular freq
    Δ   = p[4] # ω21 - ω10
    λ   = p[5] # drive 1-2 transition amplitude ~√(2)
    δω  = p[6]
    ϕ   = p[7] # phase
    A   = p[8] # amplitude, (π for π-pulse, etc...)
    Δ = Δ - δω
    B = inv(√(2π)*σ*erf(tg/(√(8)*σ))-tg*gaussian(0.5tg,σ))
    Ɛπ = A*B*(gaussian(t,σ)-gaussian(0.5tg,σ))
    Ɛπ′ = A*B*gaussianprime(t,σ)
    C = (Ɛπ - λ^2*(λ^2-4)*Ɛπ^3/32Δ^2)
    S = (-λ^2*Ɛπ′/4Δ)
    I = C*cos(δω*t) - S*sin(δω*t)
    Q = S*cos(δω*t) + C*sin(δω*t)
    ωt = ω10*t
    return I*cos(ωt+ϕ) + Q*sin(ωt+ϕ)
end

function dragpulse_cos(t::Real,p::Vector)
    # 2nd order DRAG-enhanced cosine pulse (Y-only) centered on t=0
    # 10.1103/PhysRevA.83.012308
    σ   = p[1] # gaussian standard deviation (FWHM = 2√(2ln2)σ ≈ 2.355σ)
    tg  = p[2] # total gate time
    ω10 = p[3] # qubit 0-1 transition angular freq
    Δ   = p[4] # ω21 - ω10
    λ   = p[5] # drive 1-2 transition amplitude ~√(2)
    ϕ   = p[6] # phase
    A   = p[7] # amplitude, (π for π-pulse, etc...)
    B = 1
    Ɛπ = A*(cos(2π*t/tg)+1)/tg
    Ɛπ′ = -A*sin(2π*t/tg)*2π/tg^2
    ωt = ω10*t
    return (Ɛπ - λ^2*(λ^2-4)*Ɛπ^3/32Δ^2)*cos(ωt+ϕ) + (-λ^2*Ɛπ′/4Δ)*sin(ωt+ϕ)
end

function dragpulse5(t::Real,p::Vector)
    # 5th order DRAG-enhanced gaussian pulse centered on t=0
    σ   = p[1] # gaussian standard deviation (FWHM = 2√(2ln2)σ ≈ 2.355σ)
    tg  = p[2] # total gate time
    ω10 = p[3] # qubit 0-1 transition angular freq
    Δ   = p[4] # ω21 - ω10
    λ   = p[5] # drive 1-2 transition amplitude ~√(2)
    ϕ   = p[6] # phase
    A   = p[7] # amplitude, (π for π-pulse, etc...)
    B = inv(√(2π)*σ*erf(tg/(√(8)*σ))-tg*gaussian(0.5tg,σ))
    Ɛπ = A*B*(gaussian(t,σ)-gaussian(0.5tg,σ))
    Ɛπ′ = A*B*gaussianprime(t,σ)
    Ɛx = Ɛπ + (λ^2-4)*Ɛπ^3/8Δ^2 - (13λ^4-76λ^2+112)*Ɛπ^5/128Δ^4
    Ɛy = -Ɛπ′/Δ + 33(λ^2-2)*Ɛπ^2*Ɛπ′/24Δ^3
    δ = (λ^2-4)*Ɛπ^2/4Δ - (λ^4-7λ^2+12)*Ɛπ^4/16Δ^3
    ωt = (ω10-δ)*t
    return Ɛx*cos(ωt+ϕ) + Ɛy*sin(ωt+ϕ)
end

function simplegaussian(t::Real,p::Vector)
    # Gaussian pulse centered on t=0
    FWHM = p[1] # full width at half-maximum
    ω10  = p[2] # qubit 0-1 transition angular freq
    ϕ    = p[3] # phase
    A    = p[4] # amplitude
    #B = sqrtln16/(√(π)*FWHM)
    return A*gaussian_fwhm(t,FWHM)*cos(ω10*t+ϕ)
end

function simpledragpulse(t::Real,p::Vector)
    # 1st order DRAG-enhanced gaussian pulse centered on t=0
    FWHM = p[1] # full width at half-maximum
    ω10  = p[2] # qubit 0-1 transition angular freq
    Δ    = p[3] # ω21 - ω10
    λ    = p[4] # drive 1-2 transition amplitude ~√(2)
    ϕ    = p[5] # phase
    A    = p[6] # amplitude, (π for π-pulse, etc...)
    B = sqrtln16/(√(π)*FWHM)
    Ɛx = A*gaussian_fwhm(t,FWHM)*B
    Ɛy = -A*gaussianprime_fwhm(t,FWHM)*B/Δ
    δ = (λ^2-4)*Ɛx^2/4Δ
    ωt = (ω10-δ)*t
    return Ɛx*cos(ωt+ϕ) + Ɛy*sin(ωt+ϕ)
end

function simpledragpulse5(t::Real,p::Vector)
    # 5th order DRAG-enhanced gaussian pulse centered on t=0
    FWHM = p[1] # full width at half-maximum
    ω10  = p[2] # qubit 0-1 transition angular freq
    ω21  = p[3] # qubit 1-2 transition angular freq
    λ    = p[4] # drive 1-2 transition amplitude ~√(2)
    ϕ    = p[5] # phase
    A    = p[6] # amplitude, (π for π-pulse, etc...)
    Δ = ω21-ω10
    B = sqrtln16/(√(π)*FWHM)
    Ɛπ = A*gaussian_fwhm(t,FWHM)*B
    Ɛπ′ = A*gaussianprime_fwhm(t,FWHM)*B
    Ɛx = Ɛπ + (λ^2-4)*Ɛπ^3/8Δ^2 - (13λ^4-76λ^2+112)*Ɛπ^5/128Δ^4
    Ɛy = -Ɛπ′/Δ + 33(λ^2-2)*Ɛπ^2*Ɛπ′/24Δ^3
    δ = (λ^2-4)*Ɛx^2/4Δ - (λ^4-7λ^2+12)*Ɛπ^4/16Δ^3
    ωt = (ω10-δ)*t
    return Ɛx*cos(ωt+ϕ) + Ɛy*sin(ωt+ϕ)
end

function stepfun(t::Real,p::Vector,u::Vector)
    tg = p[1] # total gate time
    to = p[2] # time offset (to make the time start at 0, e.g. if tspan=(-2,3), to=2)
    return u[round(Int,1+(length(u)-1)*(t+to)/tg)]
end


text:
function serial_add()
s=0.0
for i = 1:1000000
s=s+randn()
end
return s
end
function parallel_add()
return @parallel (+) for i=1:1000000
randn()
end
end

max_stars_repo_path: Module 2/Chapter07/parallel_for.jl

max_issues_repo_path: Module 2/Chapter07/parallel_for.jl

max_forks_repo_path: Module 2/Chapter07/parallel_for.jl


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function GitAnnotated(repo::GitRepo, commit_id::GitHash)
    ann_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_annotated_commit_lookup, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{GitHash}),
                   ann_ptr_ptr, repo.ptr, Ref(commit_id))
    return GitAnnotated(repo, ann_ptr_ptr[])
end

function GitAnnotated(repo::GitRepo, ref::GitReference)
    ann_ref_ref = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_annotated_commit_from_ref, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}),
                   ann_ref_ref, repo.ptr, ref.ptr)
    return GitAnnotated(repo, ann_ref_ref[])
end

function GitAnnotated(repo::GitRepo, fh::FetchHead)
    ann_ref_ref = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_annotated_commit_from_fetchhead, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring, Ptr{GitHash}),
                   ann_ref_ref, repo.ptr, fh.name, fh.url, Ref(fh.oid))
    return GitAnnotated(repo, ann_ref_ref[])
end

function GitAnnotated(repo::GitRepo, comittish::AbstractString)
    obj = GitObject(repo, comittish)
    cmt = peel(GitCommit, obj)
    return GitAnnotated(repo, GitHash(cmt))
end

function GitHash(ann::GitAnnotated)
    unsafe_load(ccall((:git_annotated_commit_id, :libgit2), Ptr{GitHash}, (Ptr{Void},), ann.ptr))
end

function merge_analysis(repo::GitRepo, anns::Vector{GitAnnotated})
    analysis = Ref{Cint}(0)
    preference = Ref{Cint}(0)
    anns_ref = Ref(map(a->a.ptr, anns))
    anns_size = Csize_t(length(anns))
    @check ccall((:git_merge_analysis, :libgit2), Cint,
                  (Ptr{Cint}, Ptr{Cint}, Ptr{Void}, Ptr{Ptr{Void}}, Csize_t),
                   analysis, preference, repo.ptr, anns_ref, anns_size)
    return analysis[], preference[]
end

"""Fastforward merge changes into current head """
function ffmerge!(repo::GitRepo, ann::GitAnnotated)
    cmt = GitCommit(repo, GitHash(ann))

    checkout_tree(repo, cmt)
    with(head(repo)) do head_ref
        cmt_oid = GitHash(cmt)
        msg = "libgit2.merge: fastforward $(string(cmt_oid)) into $(name(head_ref))"
        new_head_ref = if reftype(head_ref) == Consts.REF_OID
            target!(head_ref, cmt_oid, msg=msg)
        else
            GitReference(repo, cmt_oid, fullname(head_ref), msg=msg)
        end
        close(new_head_ref)
    end
    return true
end

""" Merge changes into current head """
function merge!(repo::GitRepo, anns::Vector{GitAnnotated};
                merge_opts::MergeOptions = MergeOptions(),
                checkout_opts::CheckoutOptions = CheckoutOptions())
    anns_size = Csize_t(length(anns))
    @check ccall((:git_merge, :libgit2), Cint,
                  (Ptr{Void}, Ptr{Ptr{Void}}, Csize_t,
                   Ptr{MergeOptions}, Ptr{CheckoutOptions}),
                   repo.ptr, map(x->x.ptr, anns), anns_size,
                   Ref(merge_opts), Ref(checkout_opts))
    info("Review and commit merged changes.")
    return true
end

"""Internal implementation of merge.
Returns `true` if merge was successful, otherwise `false`
"""
function merge!(repo::GitRepo, anns::Vector{GitAnnotated}, fastforward::Bool;
                merge_opts::MergeOptions = MergeOptions(),
                checkout_opts::CheckoutOptions = CheckoutOptions())
    ma, mp = merge_analysis(repo, anns)
    if isset(ma, Cint(Consts.MERGE_ANALYSIS_UP_TO_DATE))
        return true # no merge - everything is up to date
    end

    ffpref = if fastforward
        Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY
    elseif isset(mp, Cint(Consts.MERGE_PREFERENCE_NONE))
        Consts.MERGE_PREFERENCE_NONE
    elseif isset(mp, Cint(Consts.MERGE_PREFERENCE_NO_FASTFORWARD))
        Consts.MERGE_PREFERENCE_NO_FASTFORWARD
    elseif isset(mp, Cint(Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY))
        Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY
    else
        throw(ArgumentError("unknown merge preference: $(mp)."))
    end

    merge_result = if ffpref == Consts.MERGE_PREFERENCE_NONE
        if isset(ma, Cint(Consts.MERGE_ANALYSIS_FASTFORWARD))
            if length(anns) > 1
                warn("Unable to perform Fast-Forward merge with mith multiple merge heads.")
                false
            else
                ffmerge!(repo, anns[1])
            end
        elseif isset(ma, Cint(Consts.MERGE_ANALYSIS_NORMAL))
            merge!(repo, anns,
                   merge_opts=merge_opts,
                   checkout_opts=checkout_opts)
        end
    elseif ffpref == Consts.MERGE_PREFERENCE_FASTFORWARD_ONLY
        if isset(ma, Cint(Consts.MERGE_ANALYSIS_FASTFORWARD))
            if length(anns) > 1
                warn("Unable to perform Fast-Forward merge with mith multiple merge heads.")
                false
            else
                ffmerge!(repo, anns[1])
            end
        else
            warn("Cannot perform fast-forward merge.")
            false
        end
    elseif ffpref == Consts.MERGE_PREFERENCE_NO_FASTFORWARD
        if isset(ma, Cint(Consts.MERGE_ANALYSIS_NORMAL))
            merge!(repo, anns,
                   merge_opts=merge_opts,
                   checkout_opts=checkout_opts)
        end
    else
        throw(ArgumentError("unknown merge analysis result: $(ma)"))
    end
    return merge_result
end

function merge_base(repo::GitRepo, one::AbstractString, two::AbstractString)
    oid1_ptr = Ref(GitHash(one))
    oid2_ptr = Ref(GitHash(two))
    moid_ptr = Ref(GitHash())
    moid = try
        @check ccall((:git_merge_base, :libgit2), Cint,
                (Ptr{GitHash}, Ptr{Void}, Ptr{GitHash}, Ptr{GitHash}),
                moid_ptr, repo.ptr, oid1_ptr, oid2_ptr)
        moid_ptr[]
    catch e
        #warn("Pkg:",path(repo),"=>",e.msg)
        GitHash()
    end
    return moid
end

max_stars_repo_name: liveinjs/julia

max_issues_repo_name: liveinjs/julia

max_forks_repo_name: liveinjs/julia


text:
filename = replace(basename(@__FILE__), ".jl" => "")
filepath = joinpath(@__DIR__, "data/", filename, "input")
input = readlines(filepath)

# part 1

let h = 0, d = 0
    for line = input
        cmd, val = split(line, " ")
        x = parse(Int, val)
        if cmd == "forward"
            h += x
        elseif cmd == "down"
            d += x
        elseif cmd == "up"
            d -= x
        end
    end
    part_one_answer = h * d
    println(part_one_answer)
end

# part 2

let a = 0, h = 0, d = 0
    for line = input
        cmd, val = split(line, " ")
        x = parse(Int, val)
        if cmd == "forward"
            h += x
            d += a * x
        elseif cmd == "down"
            a += x
        elseif cmd == "up"
            a -= x
        end
    end
    part_two_answer = h * d
    println(part_two_answer)
end


text:
using Compose, DataFrames, Gadfly

set_default_plot_size(6.6inch,3inch)

df = DataFrame(
    A=["A","B","C","D","A","C","D","B"],
    B=[1,1,1,1,2,2,2,2],
    C=[ 12059, 57263, 79003, 13125, 25268, 13291, 17365, 30154]
)
df.xmin = df.B .- 0.4
df.xmax = df.B .+ 0.4

theme1 = Theme(bar_spacing=0.2cx)
p1 = plot(df, xmin=:xmin, xmax=:xmax, y=:C, color=:A, Geom.bar(position=:stack))
p2 = plot(df, x=:B, y=:C, color=:A, Geom.bar(position=:stack), theme1)
hstack(p1, p2)


text:
module GridapDistributed

using Gridap
using Gridap.Helpers
using Gridap.Algebra
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.CellData
using Gridap.Visualization
using Gridap.FESpaces
using Gridap.MultiField

using PartitionedArrays
const PArrays = PartitionedArrays

using SparseArrays
using WriteVTK
using FillArrays

import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
import LinearAlgebra: det, tr, cross, dot, ⋅
import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj
import Gridap.Fields: grad2curl

export FullyAssembledRows
export SubAssembledRows

include("Algebra.jl")

include("Geometry.jl")

include("CellData.jl")

include("Visualization.jl")

include("FESpaces.jl")

include("DivConformingFESpaces.jl")

include("MultiField.jl")

end # module


text:
using FASTX

DATA_DIR = "../data"
SEQ_DIR = "$DATA_DIR/sequences"
INT_DIR = "$DATA_DIR/interactions"
DSCRIPT_DIR = "$DATA_DIR/dscript"

TRAINING_SET_PATH = "$INT_DIR/training.tsv"
TEST_SET_PATH = "$INT_DIR/test.tsv"
SEQUENCES_PATH = "$SEQ_DIR/human_sequences_clustered.fasta"
ALL_SEQUENCES_PATH = "$SEQ_DIR/human_sequences.fasta"
PEPTIDES_PATH = "$SEQ_DIR/peptide_sequences.tsv"
FORMATTED_SEQUENCE_PATH = "$DSCRIPT_DIR/formatted_sequences_with_peptides.fasta"

println("Load the protein sequences...")
proteins = []
for record in FASTA.Reader(open(SEQUENCES_PATH, "r"))
    protein_id = FASTA.identifier(record)
    push!(proteins, ">$(protein_id)\n$(FASTA.sequence(record))")
end

println("Loading the peptide sequences...")
peptides = [split(p, "\t") for p in readlines(PEPTIDES_PATH)]
all_proteins = [proteins; [">$(peptide[1])\n$(peptide[2])" for peptide in peptides]]

println("Saving the sequences...")
formatted_sequences = join(all_proteins, "\n")
write(open(FORMATTED_SEQUENCE_PATH, "w"), formatted_sequences)

train_target_path = "$DSCRIPT_DIR/train.tsv"
cp(TRAINING_SET_PATH, train_target_path, force=true)

test_target_path = "$DSCRIPT_DIR/test.tsv"
cp(TEST_SET_PATH, test_target_path, force=true)

peptide_interactions = "$DSCRIPT_DIR/peptide_interactions.tsv"
open(peptide_interactions, "w") do io
    write(io, "p1\tp2\tlabel\n")
    for peptide in peptides
        for protein in proteins
            peptide_id = peptide[1]
            protein_id = split(protein, "\t")[1]
            write(io, "$(peptide_id)\t$(protein_id)\n")
        end
    end
end

max_stars_repo_name: GreenCUBIC/PeptideScreening

max_issues_repo_name: GreenCUBIC/PeptideScreening

max_forks_repo_name: GreenCUBIC/PeptideScreening


text:
import Primes

@testset "Constants" begin
    @test Groebner.FIRST_COMPUTE_PRIME == 2^31-1
    @test Groebner.FIRST_CHECK_PRIME == 2^30 + 3
end

@testset "Lucky primes" begin
    v = [[BigInt(1), BigInt(2), BigInt(7)*123, BigInt(9)*(2^31-1)], [BigInt(11)*4]]
    vc = deepcopy(v)

    @test Groebner.isluckyprime(v, 5)
    @test !Groebner.isluckyprime(v, 11)
    @test !Groebner.isluckyprime(v, 2)
    @test !Groebner.isluckyprime(v, 7)
    @test !Groebner.isluckyprime(v, 2^31-1)
    @test Groebner.isluckyprime(v, 13)

    @test Groebner.nextluckyprime([BigInt[]]) == Groebner.FIRST_COMPUTE_PRIME
    @test Groebner.nextluckyprime(v) == Groebner.nextluckyprime(v, 2^31-1) == Primes.nextprime(2^31)
    @test Groebner.nextluckyprime([[Primes.nextprime(BigInt(2^31))]]) == 2^31-1
    @test Groebner.nextluckyprime([[Primes.nextprime(BigInt(2^31))]], 2^31-1) == Primes.nextprime(1+Primes.nextprime(BigInt(2^31)))

    @test v == vc
end

@testset "Good primes" begin
    v = [[BigInt(1), BigInt(2), BigInt(7)*123, BigInt(9)*(2^31-1)], [BigInt(11)*4]]
    vc = deepcopy(v)

    moduli = [17]
    @test Groebner.nextgoodprime(v, moduli, 16) == 19
    @test Groebner.nextgoodprime(v, [17, 19], 16) == 23
    @test Groebner.nextgoodprime(v, Int[], 17) == 19
    @test Groebner.nextgoodprime(v, Int[]) == Groebner.FIRST_CHECK_PRIME

    @test v == vc
end

@testset "Scaling coefficients" begin
    v = [BigInt(1) // 3, BigInt(2) // 6, BigInt(-5) // 12, BigInt(-3)//12]
    vv = Groebner.scale_denominators(v)
    @test vv == BigInt[4, 4, -5, -3]

    v = [
            [BigInt(1)//5, BigInt(2)//3, BigInt(7)//11],
            [BigInt(11), 1],
            [BigInt(-2)//1, BigInt(-1)//2],
            [BigInt(-1)//7, BigInt(3)//14, BigInt(-12)//7],
            [BigInt(1)//2],
            [BigInt(4115)//4115226, 1]
    ]
    vc = deepcopy(v)
    vv = Groebner.scale_denominators(v)

    @test vv == [
            [BigInt(33), BigInt(55)*2, BigInt(15)*7],
            [BigInt(11), 1],
            [BigInt(-4), BigInt(-1)],
            [BigInt(-2), BigInt(3), BigInt(-24)],
            [BigInt(1)],
            [BigInt(4115), 4115226]
    ]
    @test vc == v
end

@testset "Reconstructing coefficients modulo" begin
    v = [[BigInt(3), BigInt(4)], [BigInt(5), BigInt(6)]]
    qq = Vector{Vector{Rational{BigInt}}}(undef, 0)
    modulo = BigInt(11)
    vc = deepcopy(v)

    Groebner.reconstruct_modulo!(qq, v, modulo)
    @test v == vc
    @test qq == [[-2//3, 1//3], [-1//2, 1//2]]

    Groebner.reconstruct_modulo!(qq, v, modulo)
    @test v == vc
    @test qq == [[-2//3, 1//3], [-1//2, 1//2]]

    v = [[BigInt(9), BigInt(10)], [BigInt(11), BigInt(12)]]
    modulo = BigInt(37)
    vc = deepcopy(v)
    Groebner.reconstruct_modulo!(qq, v, modulo)
    @test v == vc
    @test qq == [[-1//4, 3//4], [-4//3, -1//3]]

end

@testset "Reconstructing coefficients crt" begin
    modulo = BigInt(1)
    ring_ff = Groebner.PolyRing(2, 3, :deglex, 5, :input)
    coeffs_ff = [[UInt64(1), UInt64(2)], [UInt64(4), UInt64(3)]]
    gbcoeffs_accum = Vector{Vector{BigInt}}(undef, 0)
    coeffs_ff_c = deepcopy(coeffs_ff)

    Groebner.reconstruct_crt!(gbcoeffs_accum, modulo, coeffs_ff, ring_ff)

    @test coeffs_ff_c == coeffs_ff
    @test gbcoeffs_accum == [[BigInt(1), 2], [BigInt(4), 3]]
    @test modulo == 5

    ring_ff = Groebner.PolyRing(2, 3, :deglex, 7, :input)
    coeffs_ff = [[UInt64(6), UInt64(1)], [UInt64(5), UInt64(3)]]
    coeffs_ff_c = deepcopy(coeffs_ff)

    Groebner.reconstruct_crt!(gbcoeffs_accum, modulo, coeffs_ff, ring_ff)
    @test coeffs_ff_c == coeffs_ff
    @test gbcoeffs_accum == [[BigInt(6), 22], [BigInt(19), 3]]
    @test modulo == 35
end

@testset "Reducing coefficients" begin
    coeffs = [
        BigInt(0), BigInt(100), BigInt(99), BigInt(6),
        BigInt(3), BigInt(4), BigInt(5), BigInt(-1),
        BigInt(-2), BigInt(-3), BigInt(-4), BigInt(-5),
        BigInt(-6), BigInt(-7), BigInt(-8), BigInt(-9),
        BigInt(-10), BigInt(-11), BigInt(-12), BigInt(-13),
        BigInt(-14), BigInt(-15), BigInt(-16), BigInt(-17)
    ]
    coeffs_v = deepcopy(coeffs)
    coeffs_ff = Vector{UInt64}(undef, length(coeffs))
    prime = 5

    Groebner.reduce_modulo!(coeffs, prime, coeffs_ff)
    @test coeffs == coeffs_v
    @test coeffs_ff == UInt64[
        0, 0, 4, 1,
        3, 4, 0, 4,
        3, 2, 1, 0,
        4, 3, 2, 1,
        0, 4, 3, 2,
        1, 0, 4, 3
    ]

    ring, coeffs_ff = Groebner.reduce_modulo([coeffs], Groebner.PolyRing(2, 3, :a, 4, :b), prime)
    @test coeffs == coeffs_v
    @test ring.ch == 5
    @test coeffs_ff == [UInt64[
        0, 0, 4, 1,
        3, 4, 0, 4,
        3, 2, 1, 0,
        4, 3, 2, 1,
        0, 4, 3, 2,
        1, 0, 4, 3
    ]]

end


text:
using Why
using Test

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


text:
"""
#/* compute Laplace coefficients and Leverrier derivative/
#          j
#     j   d     i
#    a   ---   b (a)
#          j    s
#        da
#
#   by series summation */
#/* Code due to Jack Wisdom */
"""
function laplace_wisdom(s::Rational,i::Integer,j::Integer,a::Number)
# function laplace_wisdom,s,i,j,a  IDL
# double laplace(double s, int i, int j, double a);  c

##define LAPLACE_EPS 1.0e-12
const LAPLACE_EPS = convert(eltype(a),1.0e-12)

#if (i lt 0) then i = -i
i=abs(i)

if j <= i     #/* compute first term in sum */
  factor4 = one(a)
  for k=0:j-1
    factor4 *= i - k
  end
  lap_coef_sum = factor4
  q0 = 0
else
  q0 = fld(j + 1 - i,2)
  lap_coef_sum = zero(a)
  factor4 = one(a)
end

#  /* compute factors for terms in lap_coef_sum */

factor1 = s
factor2 = s + i
factor3 = i + 1
for q=1:q0-1       #/* no contribution for q = 0 */
  factor1 *= s + q
  factor2 *= s + i + q
  factor3 *= i+q+1
end
if q0 > 1
  q=q0
else
  q=1
end
#println(j+1-i,q0)
term = a*a * factor1 * factor2 / (factor3 * q)

#  /* sum series */

while  (term*factor4) > LAPLACE_EPS
  factor4 = one(a)
  for k=0:j-1
    factor4 *= (2*q + i - k)
  end
  lap_coef_sum += term * factor4
  factor1 += 1
  factor2 += 1
  factor3 += 1
  q = q+1
  term *= a*a * factor1 * factor2 / (factor3 * q)
end

#  /* fix coefficient */

for k=0:i-1
  lap_coef_sum *= (s+k)/(k+1)
end

apower = (q0 <= 0) ?  i : 2*q0 + i - 2
lap_coef_sum *= 2 * a^apower
# Return the Laplace Coefficient:
return lap_coef_sum
end

max_stars_repo_name: ericagol/TTVFaster

max_issues_repo_name: ericagol/TTVFaster

max_forks_repo_name: ericagol/TTVFaster


text:
################################################################################
# --------------- Base type for gradient-based generator:
################################################################################
"""
    AbstractGradientBasedGenerator

An abstract type that serves as the base type for gradient-based counterfactual generators.
"""
abstract type AbstractGradientBasedGenerator <: AbstractGenerator end

# ----- Julia models -----
"""
    ∂ℓ(generator::AbstractGradientBasedGenerator, M::Union{Models.LogisticModel, Models.BayesianLogisticModel}, counterfactual_state::CounterfactualState)

The default method to compute the gradient of the loss function at the current counterfactual state for gradient-based generators. It assumes that `Zygote.jl` has gradient access.
"""
function ∂ℓ(generator::AbstractGradientBasedGenerator, M::Models.Models.AbstractDifferentiableJuliaModel, counterfactual_state::CounterfactualState)
    gradient(() -> ℓ(generator, counterfactual_state), Flux.params(counterfactual_state.x′))[counterfactual_state.x′]
end

# ----- RTorch model -----
using RCall
"""
    ∂ℓ(generator::AbstractGradientBasedGenerator, M::Models.RTorchModel, counterfactual_state::CounterfactualState)

The default method to compute the gradient of the loss function at the current counterfactual state for gradient-based generators. It assumes that `Zygote.jl` has gradient access.
"""
function ∂ℓ(generator::AbstractGradientBasedGenerator, M::Models.RTorchModel, counterfactual_state::CounterfactualState)
    nn = M.nn
    x_cf = counterfactual_state.x′
    t = counterfactual_state.target_encoded
    R"""
    x <- torch_tensor($x_cf, requires_grad=TRUE)
    output <- $nn(x)
    obj_loss <- nnf_binary_cross_entropy_with_logits(output,$t)
    obj_loss$backward()
    """
    grad = rcopy(R"as_array(x$grad)")
    return grad
end

# ----- PyTorch model -----
using PyCall
"""
    ∂ℓ(generator::AbstractGradientBasedGenerator, M::Models.RTorchModel, counterfactual_state::CounterfactualState)

The default method to compute the gradient of the loss function at the current counterfactual state for gradient-based generators. It assumes that `Zygote.jl` has gradient access.
"""
function ∂ℓ(generator::AbstractGradientBasedGenerator, M::Models.PyTorchModel, counterfactual_state::CounterfactualState)
    py"""
    import torch
    from torch import nn
    """
    nn = M.nn
    x′ = counterfactual_state.x′
    t = counterfactual_state.target_encoded
    x = reshape(x′, 1, length(x′))
    py"""
    x = torch.Tensor($x)
    x.requires_grad = True
    t = torch.Tensor($[t]).squeeze()
    output = $nn(x).squeeze()
    obj_loss = nn.BCEWithLogitsLoss()(output,t)
    obj_loss.backward()
    """
    grad = vec(py"x.grad.detach().numpy()")
    return grad
end

"""
    ∂h(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState)

The default method to compute the gradient of the complexity penalty at the current counterfactual state for gradient-based generators. It assumes that `Zygote.jl` has gradient access.
"""
∂h(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState) = gradient(() -> h(generator, counterfactual_state), Flux.params(counterfactual_state.x′))[counterfactual_state.x′]

# Gradient:
"""
    ∇(generator::AbstractGradientBasedGenerator, M::Models.AbstractDifferentiableModel, counterfactual_state::CounterfactualState)

The default method to compute the gradient of the counterfactual search objective for gradient-based generators. It simply computes the weighted sum over partial derivates. It assumes that `Zygote.jl` has gradient access.
"""
∇(generator::AbstractGradientBasedGenerator, M::Models.AbstractDifferentiableModel, counterfactual_state::CounterfactualState) = ∂ℓ(generator, M, counterfactual_state) + generator.λ * ∂h(generator, counterfactual_state)

"""
    generate_perturbations(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState)

The default method to generate feature perturbations for gradient-based generators through simple gradient descent.
"""
function generate_perturbations(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState)
    𝐠ₜ = ∇(generator, counterfactual_state.M, counterfactual_state) # gradient
    Δx′ = - (generator.ϵ .* 𝐠ₜ) # gradient step
    return Δx′
end

"""
    mutability_constraints(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState)

The default method to return mutability constraints that are dependent on the current counterfactual search state. For generic gradient-based generators, no state-dependent constraints are added.
"""
function mutability_constraints(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState)
    mutability = counterfactual_state.params[:mutability]
    return mutability # no additional constraints for GenericGenerator
end

"""
    conditions_satisified(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState)

The default method to check if the all conditions for convergence of the counterfactual search have been satisified for gradient-based generators.
"""
function conditions_satisified(generator::AbstractGradientBasedGenerator, counterfactual_state::CounterfactualState)
    𝐠ₜ = ∇(generator, counterfactual_state.M, counterfactual_state)
    status = all(abs.(𝐠ₜ) .< generator.τ)
    return status
end

##################################################
# Specific Generators
##################################################

include("GenericGenerator.jl") # Wachter et al. (2017)
include("GreedyGenerator.jl") # Schut et al. (2021)
include("REVISEGenerator.jl") # Joshi et al. (2019)


text:
# Autogenerated wrapper script for rtmidi_jll for i686-w64-mingw32-cxx03
export rtmidi

JLLWrappers.@generate_wrapper_header("rtmidi")
JLLWrappers.@declare_library_product(rtmidi, "librtmidi-6.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        rtmidi,
        "bin\\librtmidi-6.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using FractionalDiffEq, Plots
α=[0.94, 0.94, 0.94]; ϕ=[0.2, 0, 0.5]; τ=0.009; T=1.4; h=0.001
function delaychen!(dy, y, ϕ, t)
	a=35; b=3; c=27
	dy[1] = a*(y[2]-ϕ[1])
	dy[2] = (c-a)*ϕ[1]-y[1]*y[3]+c*y[2]
	dy[3] = y[1]*y[2]-b*ϕ[3]
end
prob = FDDESystem(delaychen!, ϕ, α, τ, T)
y, x=solve(prob, h, DelayABM())
plot(x[:, 1], x[:, 2], x[:, 3], title="Fractional Order Chen Delayed System")


text:
# Meshgrid function (ndgrid)
function meshgrid{T}(vx::Array{T}, vy::Array{T})
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

function smart_innerY{T<:Signed, F<:AbstractFloat}(Y::Array{F}, N::T, U::Array{F}, L::Array{F}, σ1::F, μ1::F)
    ## Here, we assume that the vector Y is already normalized.

    ## The computation of lambda max involves taking all inner products of Y.
    ## The smart method allows for the computation of the i+1_th parcel based on
    ## the result of the i_th parcel.

    arraytype = typeof(Y[1])

    ## Preallocating the vector IP_vector for speed
    IP_vector = Array{arraytype}(zeros(N,1))

    index_array = collect(1:N)'
    IP_vector[1]= sum( Y.* (index_array'-μ1) / σ1 )
    sumY = 0

    for i = 2:N
        sumY = sumY + Y[i-1]
        IP_vector[i] = sumY * ( U[i-1] - L[i-1] )
    end

    return IP_vector
end

function calculate_matrices{T<:Int}(N::T)
    ## index
    i = collect(1:N-1)'

    ## calculate arrays mean and variance
    array_μ = (N-i')/N
    array_σ = sqrt( (1/(N-1)) *( (i').*array_μ .^2 + (N-i').*(1-array_μ).^2 ) )

    ## With these parameters, compose the Lower and Upper diagonal matrices
    U = - array_μ ./ array_σ
    L = (1 - array_μ) ./ array_σ

    ## Calculate μ1 and σ1
    μ1 = (N+1)/2
    σ1 = sqrt( (1/6)*(N+1)*(2*N+1)-μ1^2 )

    ## Concatenate
    array_σ2 = [σ1; array_σ]

    ## Build X
    repU, repLtranspose = meshgrid(U, L)
    #x1 = ( ((1:N)' - μ1) / σ1 )'
    #x2 = triu(repU) + tril(repLtranspose') - diagm(L[:])
    #x3 = [x1 [x2; L']]
    X = [( ((1:N)' - μ1) / σ1 )' [triu(repU) + tril(repLtranspose') - diagm(L[:]); L']]

    ## Calculate Gram Matrix (inner products matrix)
    GMat = triu( GramMatrixTM(N, U, L, σ1, μ1) );
    GMat = GMat + GMat' - N*eye(N);

    return U, L, μ1, σ1, X, GMat, array_σ2
end

function GramMatrixTM{T<:Int, F<:AbstractFloat}(N::T, U::Array{F}, L::Array{F}, σ1::F, μ1::F)
    ## builds the inner products matrix in the smart way
    ## (much faster than X'*X)
    I, J = meshgrid(collect(1:(N-1)), collect(1:(N-1)))
    Uaux, Laux = meshgrid(U, L)

    minimo = min(I, J);
    maximo = max(I, J);

    GMat1 = minimo .* (U*U') + (maximo-minimo) .* (Laux.*Uaux) + (N-maximo) .* (L*L')

    ## reduce the memory usage, since we don't need theses variables anymore...
    I, J, minimo, maximo = 0, 0, 0, 0
    gc()

    j = collect(1:N-1)';
    IP = (1/σ1) * ( j .* U' .* ((j+1)/2 - μ1) + (N-j) .* L' .* ((j+1+N)/2 - μ1) )

    GMat = [[N; IP'] [IP; GMat1]]
end

function calculate_μ{T<:Int}(N::T)
    i = 1:N-1
    (N-i')/(N)
end

function calculate_μ1{T<:Int}(N::T)
    (N+1)/2
end

function calculate_σ{T<:Int, F<:AbstractFloat}(N::T,array_μ::Array{F})
    i = 1:N-1
    sqrt( (1/(N- 1)) *( (i').*array_μ .^2 + (N-i').*(1-array_μ).^2 ) )
end

function calculate_σ1{T<:Int, F<:AbstractFloat}(N::T,μ1::F)
    sqrt( (1/6)*(N+1)*(2*N+1)-μ1^2 )
end

function CoordinateDescentCore{T<:Signed, F<:AbstractFloat}(A_actual::Array{T}, GMat::Array{F, 2}, j::T, β_tilde::Array{F}, β_ols::Array{F}, N::T, arrayInnerY::Array{F}, λ::F, array_σ2::Array{F})
    # aux = 0
    #
    # if size(A_actual)[1] >= 1
    #     aux = ( GMat[j,A_actual]' * β_tilde[A_actual] )[1]
    # end
    aux = sum(GMat[j,A_actual]' * β_tilde[A_actual])

    ##  Computation of the OLS estimator as simple regression
    β_ols[j] = 1 / N * (arrayInnerY[j] - aux) + β_tilde[j]

    ##  Soft-thresholding operator
    if j == 1
        λEff = λ
    else
        λEff = λ * array_σ2[j]
    end

    return aux, λEff, β_ols
end

function updateCandidates{T<:Signed, F<:AbstractFloat}(λEff::F, β_ols::Array{F}, j::T, A_actual::Array{T}, lastchange::T, β_tilde::Array{F}; w = [])
    breakThis = 0

    if size(w,1) == 0
        w = ones( size(β_ols) )
    end

    if λEff*w[j] >= abs(β_ols[j])
        if in(j, A_actual) # any(j == A_actual)
            splice!(A_actual,findin(A_actual,j)[1])
            lastchange = j
        elseif j == lastchange
            breakThis = 1
        end
        β_tilde[j] = 0
    else
        β_tilde[j] = sign(β_ols[j]) * (abs(β_ols[j]) - λEff*w[j])
        ## Includes an index
        if !in(j, A_actual) #not(any(j==A_actual))
            push!(A_actual, Int(j))
            sort!(A_actual)
            lastchange = j
        elseif j == lastchange
            breakThis = 1
        end
    end

    return A_actual, β_tilde, lastchange, breakThis
end

function convCheck{T<:Signed}(A_reference::Array{T}, A_actual::Array{T})
    if length(A_reference) != length(A_actual)
        ##  Did not converge
        convergenceFlag = 0
    elseif A_reference == A_actual && length(A_reference) == length(A_actual)
        ##  Converged!
        convergenceFlag = 2
    end

    return convergenceFlag
end

function calcBestB{T<:Signed, F<:AbstractFloat}(N::T, Y::Array{F}, X::Array{F}, β_lasso::Array{F}, array_σ2::Array{F})
    bestB = -1
    BICold = 1000000

    for i = 1:size(β_lasso, 2)
        BICnew = N * ( log(var( Y - X*β_lasso[:,i] )) ) + sum( abs(β_lasso[:,i]./array_σ2) .> 1e-3 ) * log(N)
        if BICnew < BICold
            BICold = BICnew
            bestB = i
        end
    end

    return bestB
end


text:
"""
Abstract supertype for parameters. Currently, the subtypes are `Parameters` and `BoundedParameter`.
"""
abstract type AbstractParameter end
"""
Abstract supertype for different definitions of parameter sets. Currently the subtypes are `ParameterGroup` and `ParameterSpace`.
"""
abstract type AbstractParameterSet end
"""
Abstract supertype for structs that contain multiple `Parameters` and/or `BoundedParameters`. Useful for printing a single struct that holds parameters.
"""
abstract type ParameterGroup <: AbstractParameterSet end
"""
Abstract supertype that contains all the parameters in the model. Useful for printing the parameter space.
"""
abstract type ParameterSpace <: AbstractParameterSet end
"""
A simple wrapper that contains the value of the parameter and a description.
"""
struct Parameter{T,S<:String} <: AbstractParameter
    value::T
    description::S
end
"""
A simple wrapper that contains the value of the parameter, the lower and upper bounds, and a description.
"""
struct BoundedParameter{T,S<:String} <: AbstractParameter
    value::T
    lb::T
    ub::T
    description::S
end
"""
    Base.getproperty(name::ParameterGroup, val::Symbol)

Extension of `getproperty` function in order to call the value of an `AbstractParameter` within a `ParameterGroup` as T.x instead of T.x.value
"""
function Base.getproperty(name::ParameterGroup, val::Symbol)
    if typeof(getfield(name, val)) <: AbstractParameter
        return getfield(getfield(name, val), :value)
    else
        return getfield(name, val)
    end
end
"""
A function that returns the description of an `AbstractParameter` within a `ParameterGroup`. This is necessary because of the definition of getproperty for `ParameterGroup` that returns the value of the parameter instead of `AbstractParameter`. For example, T.x returns x.value and not x, which contains the description in the wrapper.
"""
function description(name::ParameterGroup, val::Symbol)
    if typeof(getfield(name, val)) <: AbstractParameter
        return getfield(getfield(name, val), :description)
    else
        @warn "No description provided"
        return nothing
    end
end


text:
using Documenter, Bigsimr
using Distributions

DocMeta.setdocmeta!(Bigsimr, :DocTestSetup, :(using Bigsimr); recursive=true)

# Gotta set this environment variable when using the GR run-time on CI machines.
# This happens as examples will use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"

makedocs(
    sitename = "Bigsimr.jl",
    modules  = [Bigsimr],
    doctest  = false,
    pages    = [
        "Bigsimr.jl" => "index.md",
        "Guides" => [
            "Getting Started" => "getting_started.md",
            "Pearson Matching" => "pearson_matching.md",
            "Nearest Correlation Matrix" => "nearest_correlation_matrix.md"
        ],
        "API Reference" => "main_functions.md",
        "Details" => "details.md",
        "Index" => "function_index.md"
    ]
)

deploydocs(
    repo = "github.com/SchisslerGroup/Bigsimr.jl.git",
    versions = ["stable" => "v^", "v#.#", "dev" => "master"]
)


text:
function newton_solve(f, g, H, x0; ftol=-Inf, gtol=1e-4, iter=100,
        lower=-Inf, upper=Inf, reg = 1e-3)
    x = copy(x0)
    X = [copy(x0)]
    for k = 1:iter
		fe = f(x)
		ge = g(x)
		He = H(x)
        He += reg * norm(He, Inf) * I
        ((norm(ge, Inf) < gtol) || (fe < ftol)) && break
        p = - He \ ge
        α = clamped_linesearch(f, x, p, lower=lower, upper=upper)
        x = clamp.(x + α*p, lower, upper)
        push!(X, copy(x))
        println("k:", k,
            "   f:", scn(fe, digits=3),
            "   ∇:", scn.(norm(ge, Inf)),
            "   α:", scn.(α),
            # "   g:", scn.(-g(x, rot=rot)),
            # "   p:", scn.(p),
            # "   x:", scn.(x),
            # "   H:", scn.(H),
            )
    end
    return x, X
end

function clamped_linesearch(f, x, p; iter=10,
        lower=-Inf, upper=Inf)
    α = 1.0
    fprev = f(x)
    for k = 1:iter
        xc = clamp.(x + α*p, lower, upper)
        (f(xc) <= fprev) && break
        α /= 2
    end
    return α
end

################################################################################
# Generate & Save Dataset
################################################################################
generate_hardware_dataset(N=50, sleep_ratio=0.01)

################################################################################
# Load Dataset
################################################################################
params0, trajs0, pairs0 = open_dataset(:hardwarebox; N=50)

pairs0
params0
data0 = params0[:data]

function quad_loss(x, timestep, x0, v0)
	loss = 0.0
	for k = 1:7
		t = (k-1) * timestep
		Δx = x[k][3] - (x0 + v0*t - t^2/2 * 9.81)
		loss += Δx^2
	end
	return loss
end

function pairs_loss(pairs, θ)
	loss = 0.0
	for i = 1:50
		ps = pairs[(i-1)*6 .+ (1:6)]
		x = [[ps[1][1][1:3]]; [p[2][1:3] for p in ps]]
		loss += quad_loss(x, θ[1], θ[1 + (i-1)*2 .+ (1:2)]...)
	end
	return loss / 50
end

function pairs_grad(pairs, θ)
	grad = FiniteDiff.finite_difference_gradient(θ -> pairs_loss(pairs, θ), θ)
	return grad
end

function pairs_hess(pairs, θ)
	hess = FiniteDiff.finite_difference_hessian(θ -> pairs_loss(pairs, θ), θ)
	return hess
end

θ0 = ones(101)
pairs_loss(pairs0, θ0)
pairs_grad(pairs0, θ0)
plot(Gray.(abs.(pairs_hess(pairs0, θ0))))

f00(θ) = pairs_loss(pairs0, θ)
g00(θ) = pairs_grad(pairs0, θ)
H00(θ) = pairs_hess(pairs0, θ)

θsol, θSOL = newton_solve(f00, g00, H00, θ0, gtol=1e-6, reg = 0.)
θsol


text:
x, y = rand(200), rand(200)
jtest = JointDistanceDistributionTest(B = 20)

@test causality(x, y, jtest) isa Vector{T} where T


text:
function SeisSmooth2(A;L1=7,L2=7,Nrepeat=1)
    a1 = isodd(L1) ? 0 : 1
    a2 = isodd(L2) ? 0 : 1
    kernel = ones(eltype(A),L1,L2)/(L1*L2);
    for irepeat = 1 : Nrepeat
	A = conv2(A,kernel)
	A = A[floor(Int,L1/2) + 1:end-floor(Int,L1/2) + a1,floor(Int,L2/2) + 1:end-floor(Int,L2/2) + a2]
	end
    return A
end


text:
export exportImage, exportMovie, exportMovies

exportImage(filename, im::ImageMeta; kargs...) = exportImage(filename, im.data; kargs...)

function exportImage(filename, im::AbstractMatrix{T}; vmin=0.0, vmax=1.0,
                              colormap="gray", normalize=true, kargs...) where {T<:Real}
  imC = colorize( im, vmin, vmax, cmap(colormap), normalize=normalize )
  exportImage(filename, imC; kargs...)
end

function exportImage(filename, im::AbstractMatrix{T}; pixelResizeFactor=1) where {T<:Colorant}
  file, ext = splitext(filename)

  imR = repeat(im, inner=[pixelResizeFactor,pixelResizeFactor])

  minPxSpacing = minimum(pixelspacing(im))
  newSize = ceil.(Int64, collect(pixelspacing(im)) / minPxSpacing .* collect(size(imR)) )

  dataResized = Images.imresize(imR,(newSize[1],newSize[2]))
  rgbdata = convert(Array{RGB},dataResized)
  ImageMagick.save(filename, rgbdata)
end

function exportImage(filename, data::Vector{T}; kargs...) where {T<:AbstractMatrix}
    file, ext = splitext(filename)
    exportImage(file*"_xy.png", data[1]; kargs...)
    exportImage(file*"_xz.png", data[2]; kargs...)
    exportImage(file*"_yz.png", data[3]; kargs...)
end

### export movies ###

exportMovie(filename, data::ImageMeta; kargs...) = exportMovie(filename, data.data; kargs...)

function exportMovie(filename, data::AbstractArray{T,3}; vmin=0.0, vmax=1.0,
                              colormap="gray", normalize=true, kargs...) where {T<:Real}
  imC = colorize( data, vmin, vmax, cmap(colormap), normalize=normalize )
  exportMovie(filename, imC; kargs...)
end

function exportMovie(filename, data::AbstractArray{T,3}; pixelResizeFactor=1) where {T<:Colorant}
  file, ext = splitext(filename)

  if pixelResizeFactor > 1
    data = repeat(data, inner=[pixelResizeFactor,pixelResizeFactor,1])
  end

  minPxSpacing = minimum(pixelspacing(data))
  newSize = ceil.(Int64, collect(pixelspacing(data)[1:2]) / minPxSpacing .* collect(size(data)[1:2]) )

  datai = similar(data, (newSize[1],newSize[2],size(data,3)))

  for l=1:size(data,3)
    datai[:,:,l] = Images.imresize(data[:,:,l],(newSize[1],newSize[2]))
  end

  rgbdata = convert(Array{RGB},datai)
  @debug "saving $filename"
  ImageMagick.save(file*".gif", rgbdata)
end

function exportMovies(filename, data::Vector; kargs...)
  file, ext = splitext(filename)

  exportMovie(file*"_xy.gif", data[1]; kargs...)
  exportMovie(file*"_xz.gif", data[2]; kargs...)
  exportMovie(file*"_yz.gif", data[3]; kargs...)
end


text:
using Test, ForwardDiff, Distributions, FDM, Flux.Tracker
using StatsFuns: binomlogpdf

# Real

x_real = randn(5)

dists = [Normal(0, 1)]

for dist in dists

    f(x::Vector) = sum(logpdf.(Ref(dist), x))

    ForwardDiff.gradient(f, x_real)

end

# Postive

x_positive = randn(5).^2

dists = [Gamma(2, 3)]

for dist in dists

    f(x::Vector) = sum(logpdf.(Ref(dist), x))

    g = x -> ForwardDiff.gradient(f, x)

end

let
    foo = p->binomlogpdf(10, p, 3)
    @test isapprox(
        Tracker.gradient(foo, 0.5)[1],
        central_fdm(5, 1)(foo, 0.5);
        rtol=1e-8,
        atol=1e-8,
    )
    @test isapprox(
        Tracker.gradient(foo, 0.5)[1],
        ForwardDiff.derivative(foo, 0.5);
        rtol=1e-8,
        atol=1e-8,
    )

    bar = p->logpdf(Binomial(10, p), 3)
    @test isapprox(
        Tracker.gradient(bar, 0.5)[1],
        central_fdm(5, 1)(bar, 0.5);
        rtol=1e-8,
        atol=1e-8,
    )
    @test isapprox(
        Tracker.gradient(bar, 0.5)[1],
        ForwardDiff.derivative(bar, 0.5),
        rtol=1e-8,
        atol=1e-8,
    )
end

let
    foo = p->Turing.poislogpdf(p, 1)
    @test isapprox(
        Tracker.gradient(foo, 0.5)[1],
        central_fdm(5, 1)(foo, 0.5);
        rtol=1e-8,
        atol=1e-8,
    )
    @test isapprox(
        Tracker.gradient(foo, 0.5)[1],
        ForwardDiff.derivative(foo, 0.5);
        rtol=1e-8,
        atol=1e-8,
    )

    bar = p->logpdf(Poisson(p), 3)
    @test isapprox(
        Tracker.gradient(bar, 0.5)[1],
        central_fdm(5, 1)(bar, 0.5);
        rtol=1e-8,
        atol=1e-8,
    )
    @test isapprox(
        Tracker.gradient(bar, 0.5)[1],
        ForwardDiff.derivative(bar, 0.5);
        rtol=1e-8,
        atol=1e-8,
    )
end


text:
"""
$(SIGNATURES)

Compute the area of the given `ConvexHull` using the
[shoelace formula](https://en.wikipedia.org/wiki/Shoelace_formula).
"""
function area(hull::ConvexHull)
    T = eltype(hull)
    vertices = hull.vertices
    n = length(vertices)
    n <= 2 && return zero(arithmetic_closure(T))
    @inbounds begin
        ret = cross2(vertices[n], vertices[1])
        @simd for i in Base.OneTo(n - 1)
            ret += cross2(vertices[i], vertices[i + 1])
        end
        return abs(ret) / 2
    end
end


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

_sec2rad{T<:AbstractFloat}(sec::T) = deg2rad(sec/3600.0)

"""
    sec2rad(sec) -> radians

### Purpose ###

Convert from seconds to radians.

### Argument ###

* `sec`: number of seconds.  It can be either a scalar or a vector.

### Output ###

The number of radians corresponding to `sec`.  If `sec` is an array, an array of
the same length is returned.

### Example ###

``` julia
sec2rad(3600*30)
# => 0.5235987755982988
```

### Notes ###

Use `rad2sec` to convert radians to seconds.
"""
sec2rad(sec::Real) = _sec2rad(float(sec))

@vectorize_1arg Real sec2rad


text:
"""CellList type. The `indices` dictionary maps the `CartesianIndex` of each cell to the indices of points in that cell.
"""
struct CellList{d}
    data::Dict{CartesianIndex{d}, Vector{Int}}
end

"""Construct CellList from points `p` and radius `r>0`.

# Examples
```julia
n, d, r = 100, 2, 0.01
p = rand(n, d)
c = CellLists(p, r)
```
"""
function CellList(p::Array{T, 2}, r::T; offset::Int=0) where T <: AbstractFloat
    @assert r > 0
    n, d = size(p)
    cells = @. Int(fld(p, r))
    data = Dict{CartesianIndex{d}, Vector{Int}}()
    for j in 1:n
        @inbounds cell = CartesianIndex(cells[j, :]...)
        if haskey(data, cell)
            @inbounds push!(data[cell], j + offset)
        else
            data[cell] = [j + offset]
        end
    end
    CellList{d}(data)
end

"""Compute offsets to neighboring cells in `d` dimensions."""
@inline function neighbors(d::Int)
    n = CartesianIndices((fill(-1:1, d)...,))
    return n[1:fld(length(n), 2)]
end

@inline function distance_condition(p1::Vector{T}, p2::Vector{T}, r::T) where T <: AbstractFloat
    sum(@. (p1 - p2)^2) ≤ r^2
end

@inline function brute_force!(ps::Vector{Tuple{Int, Int}}, is::Vector{Int}, p::Array{T, 2}, r::T) where T <: AbstractFloat
    for (k, i) in enumerate(is[1:(end-1)])
        for j in is[(k+1):end]
            if @inbounds distance_condition(p[i, :], p[j, :], r)
                push!(ps, (i, j))
            end
        end
    end
end

@inline function brute_force!(ps::Vector{Tuple{Int, Int}}, is::Vector{Int}, js::Vector{Int}, p::Array{T, 2}, r::T) where T <: AbstractFloat
    for i in is
        for j in js
            if @inbounds distance_condition(p[i, :], p[j, :], r)
                push!(ps, (i, j))
            end
        end
    end
end

"""Return a vector of all pairs that are in neighboring cells in the cell list.

# Examples
```julia-repl
julia> near_neighbors(c, p, r)
[(1, 4), (3, 11), ...]
```
"""
function near_neighbors(c::CellList{d}, p::Array{T, 2}, r::T) where d where T <: AbstractFloat
    ps = Vector{Tuple{Int, Int}}()
    offsets = neighbors(d)
    # Iterate over non-empty cells
    for (cell, is) in c.data
        # Pairs of points within the cell
        brute_force!(ps, is, p, r)
        # Pairs of points with non-empty neighboring cells
        for offset in offsets
            neigh_cell = cell + offset
            if haskey(c.data, neigh_cell)
                @inbounds js = c.data[neigh_cell]
                brute_force!(ps, is, js, p, r)
            end
        end
    end
    return ps
end