text:
"""
    getOfficialGeneSymbol(queryId::String; sourceDb::String, chromosomal::Bool=true)
    getOfficialGeneSymbol(queryId::Vector{String}; sourceDb::String, chromosomal::Bool=true)

Compute the official gene symbol for a single or multiple `queryId`(s).

If `sourceDb` is missing, the database will automatically try to detect the queryId's
source database. Per default `chromosomal` is set to `true`.

### Arguments
- `queryId` : can either be of type `String` or `Vector{String}`
- `sourceDb::String` : source database of the `queryId`
- `chromosomal::Bool` : if `true`, only features located on one of the 23 chromosomes ({1,2,...,X,Y})
              are considered. If `false`, features located on human alternative chromosomes will also
              be considered. Alternative sequences are sequences regions which differ from the
              genomic DNA on the primary assembly, such as HSCHR17_2_CTG4 [info on ensembl blog](http://www.ensembl.info/2011/05/20/accessing-non-reference-sequences-in-human/).

### Examples
```julia-repl
julia> getOfficialGeneSymbol("596")
WARNING: Some input IDs have multiple mappings
2×3 DataFrames.DataFrame
│ Row │ InputId │ InputSourceDb │ OfficialGeneSymbol │
├─────┼─────────┼───────────────┼────────────────────┤
│ 1   │ 596     │ NCBI gene     │ BCL2               │
│ 2   │ 596     │ WikiGene      │ BCL2               │

julia> getOfficialGeneSymbol("596", sourceDb = "WikiGene")
1×3 DataFrames.DataFrame
│ Row │ InputId │ InputSourceDb │ OfficialGeneSymbol │
├─────┼─────────┼───────────────┼────────────────────┤
│ 1   │ 596     │ WikiGene      │ BCL2               │

julia> getOfficialGeneSymbol(["BCL2", "Bcl-2", "AMPK", "NotIn", "GJE1"])
WARNING: No values found for:
NotIn
WARNING: Some input IDs have multiple mappings
WARNING: Some input IDs result in multiple official gene symbols
5×3 DataFrames.DataFrame
│ Row │ InputId │ InputSourceDb        │ OfficialGeneSymbol │
├─────┼─────────┼──────────────────────┼────────────────────┤
│ 1   │ GJE1    │ Official Gene Symbol │ GJE1               │
│ 2   │ GJE1    │ Gene Symbol Alias    │ GJC3               │
│ 3   │ Bcl-2   │ Gene Symbol Alias    │ BCL2               │
│ 4   │ BCL2    │ Official Gene Symbol │ BCL2               │
│ 5   │ AMPK    │ Gene Symbol Alias    │ PRKAA2             │
```
"""
function getOfficialGeneSymbol(
        queryId::S;
        sourceDb::S = "",
        chromosomal::B = true) where {B<:Bool, S<:String}

   return getOfficialGeneSymbol([queryId]; sourceDb = sourceDb, chromosomal = chromosomal);
end

function getOfficialGeneSymbol(
        queryId::T;
        sourceDb::S = "",
        chromosomal::B = true) where {B<:Bool, S<:String, T<:Vector{String}}

   q = "CALL rcsi.convert.table.getOfficialGeneSymbol(" *
      "{ids}, " *
      (sourceDb == "" ? "null, " : ("'" * sourceDb * "', ")) *
      "{chr}) YIELD value " *
      "RETURN " *
      "   value.InputId AS InputId, " *
      "   value.InputSourceDb AS InputSourceDb, " *
      "   value.OfficialGeneSymbol AS OfficialGeneSymbol";

   x = Neo4j.cypherQuery(
      NEO4J_CONNECTION,
      q,
      "ids" => queryId,
      "chr" => chromosomal
   );
   if DataFrames.nrow(x) > 0
      xDiff = setdiff(queryId, x[:InputId]);

      checkReturnedQueryIds(x, queryId);

      in2multiOut = @> begin
         x
         @select(:InputId, :OfficialGeneSymbol)
         unique()
         @by(:InputId, n = length(:OfficialGeneSymbol))
         @where(:n .> 1)
      end
      if !isempty(in2multiOut)
         warn("Some input IDs result in multiple official gene symbols", prefix = "WARNING: ");
      end
   end

   return x;
end

"""
    convertFromToExtended(queryId::String; targetDb::Union{String, Vector{String}} attributesUnion{String, Vector{String}}, sourceDb::String, longFormat::Bool, chromosomal::Bool)
    convertFromToExtended(queryId::Vector{String}; targetDb::Union{String, Vector{String}}, attributesUnion{String, Vector{String}}, sourceDb::String, longFormat::Bool, chromosomal::Bool)

Retrieve molecular identifier from other databases, `targetDb`, for single or mulitple query IDs, `queryId`,
and moreover information on Ensembl gene, transcript and peptide IDs, such as ID and genomic loation.

### Arguments
- `queryId` : either be of type `String` or `Vector{String}`
- `targetDb` : either of type `String` or `Vector{String}` (default: `[GeneSymbolDB]`) (keyword argument).
- `attributes` : either be of type `String` or `Vector{String}`. Please see `getValidAttributes()`
              for the kind of attributes available.
- `sourceDb::String` : source database of the `queryId` (keyword argument)
- `longFormat::Bool` : `true` -> target values are stored in two columns `TargetDb` and `TargetId`.
                    If `false`, each `TargetDb` value will be in a separate column. In the case
                    when there are multiple `TargetId` values for one `TargetDb` values will be
                    concatenated and separated by a comma (only for databases other than Ensembl) (keyword argument).
- `chromosomal::Bool` : if `true`, only features located on one of the 23 chromosomes ({1,2,...,X,Y})
              are considered. If `false`, features located on human alternative chromosomes will also
              be considered. Alternative sequences are sequences regions which differ from the
              genomic DNA on the primary assembly, such as HSCHR17_2_CTG4
              [info on ensembl blog](http://www.ensembl.info/2011/05/20/accessing-non-reference-sequences-in-human/). (keyword argument)

**See also**
`getValidAttributes()`
`getValidDatabases()`

### Examples
```julia-repl
julia>

```

"""
function convertFromToExtended(
queryId::S;
targetDb::U = ["GeneSymbolDB"],
attributes::U = Vector{String}(),
sourceDb::S = "",
longFormat::B = true,
chromosomal::B = true) where {S<:String, U<:Union{String, Vector{String}}, B<:Bool}

return convertFromToExtended(
  [queryId],
  targetDb = isa(targetDb, String) ? String[targetDb] : targetDb,
  attributes = isa(attributes, String) ? [attributes] : attributes,
  sourceDb = sourceDb,
  longFormat = longFormat,
  chromosomal = chromosomal);
end

function convertFromToExtended(
   queryId::T;
   targetDb::U = ["GeneSymbolDB"],
   attributes::U = Vector{String}(),
   sourceDb::S = "",
   longFormat::B = true,
   chromosomal::B = true) where {S<:String, T<:Vector{String}, U<:Union{String, Vector{String}}, B<:Bool}

   targetDb = isa(targetDb, String) ? String[targetDb] : targetDb;
   attributes = isa(attributes, String) ? [attributes] : attributes;

   tDb = copy(targetDb);
   append!(tDb, collect(keys(addDbColDict)));

   q = "CALL rcsi.convert.json.convertIdsExtended(" *
         "{queryId}, {tDb}, {attributes}, " *
         (sourceDb == "" ? "null, " : ("'" * sourceDb * "', ")) *
         "{chromosomal});";

   tx = Neo4j.transaction(NEO4J_CONNECTION)
   tx(q,
      "queryId" => queryId,
      "tDb" => tDb,
      "attributes" => isempty(attributes) ? String[] : attributes,
      "chromosomal" => chromosomal);
   c = Neo4j.commit(tx);

   if isempty(c.results)
      warn("message: " * c.errors[1]["message"]);
      error("code: " * c.errors[1]["code"]);
   end

   res = c.results[end];

   # Parse returned JSON object
   if isempty(res["data"])
      return DataFrames.DataFrame();
   end
   x = parseGngJson(res["data"]);

   if DataFrames.nrow(x) > 0
   # Check if some symbols were not found or returned multiple values
   keysUq = unique(x[[:InputId, :InputSourceDb]]);

   gd = DataFramesMeta.groupby(keysUq, :InputId)
   for group in gd
      DataFrames.nrow(group) == 1 && continue;
      g = unique(group[:InputSourceDb]);
      if length(g) == 1
         if g[1] == "Official Gene Symbol" && !in("ArrayExpress", targetDb)
            warn("Input ID '$(group[:InputId][1])' is an official gene symbol that codes for multiple Ensembl gene IDs. " *
                  "Possible reasons:" *
                  "\n* One of the ENSG-IDs is located on an 'Alternative Human chromosome' or " *
                  "\n* the ENSG-IDs encode different biotypes" *
                  "\n\nConsider " *
                  ((!chromosomal && !in("ArrayExpress", targetDb)) ?
                     "setting `chromosomal = false` and adding 'ArrayExpress' to targetDb parameter" :
                        (!chromosomal ? "setting `chromosomal = false`" : "adding 'ArrayExpress' to targetDb parameter")));
         elseif g[1] == "Gene Symbol Alias"
            warn("Input ID '$(group[:InputId][1])' represents multiple official gene symbols.")
         end
      elseif all(map(i->in(i,g),String["Gene Symbol Alias", "Official Gene Symbol"]))
         warn("Input ID '$(group[:InputId][1])' represents a gene symbol alias as well as an official gene symbol.")
      end
   end

   if !longFormat
      x = unstackDf(x);
   end

   # Remove id columns (id columns that are not in the targetDb are removed)
   ## it always assumes it is unstacked!!!
   tDbDiff = setdiff(tDb, targetDb);
   !isempty(tDbDiff) && map(i->delete!(x, addDbColDict[i]), tDbDiff);

   x = unique(x);

   checkReturnedQueryIds(x, queryId);
   else
      return DataFrames.DataFrame();
   end

   return x;

end

"""
    convertFromTo(queryId::String; targetDb::Union{String, Vector{String}}, sourceDb::String, longFormat::Bool=true, chromosomal::Bool=true)
    convertFromTo(queryId::Vector{String}; targetDb::Union{String, Vector{String}}, sourceDb::String, longFormat::Bool=true, chromosomal::Bool=true)

Given a single or multiple input IDs, `queryId`, convert to a single or multiple target databases, `targetDb`.
If `sourceDb` is missing, the database will automatically try to detect the queryId's
source database. Per default `chromosomal` is set to `true`.

### Arguments
- `queryId` : either be of type `String` or `Vector{String}`
- `targetDb` : either of type `String` or `Vector{String}` (default: `[GeneSymbolDB]`) (keyword argument).
- `sourceDb::String` : source database of the `queryId` (keyword argument)
- `longFormat::Bool` : `true` -> target values are stored in two columns `TargetDb` and `TargetId`.
                    If `false`, each `TargetDb` value will be in a separate column. In the case
                    when there are multiple `TargetId` values for one `TargetDb` values will be
                    concatenated and separated by a comma (only for databases other than Ensembl) (keyword argument).
- `chromosomal::Bool` : if `true`, only features located on one of the 23 chromosomes ({1,2,...,X,Y})
              are considered. If `false`, features located on human alternative chromosomes will also
              be considered. Alternative sequences are sequences regions which differ from the
              genomic DNA on the primary assembly, such as HSCHR17_2_CTG4
              [info on ensembl blog](http://www.ensembl.info/2011/05/20/accessing-non-reference-sequences-in-human/). (keyword argument)

**See also**
`getValidDatabases()`

### Examples
```julia-repl
julia> convertFromTo(["IGHV1-2", "Bcl-2"])
2×4 DataFrames.DataFrame
│ Row │ InputId │ InputSourceDb        │ TargetDb             │ TargetId │
├─────┼─────────┼──────────────────────┼──────────────────────┼──────────┤
│ 1   │ IGHV1-2 │ Official Gene Symbol │ Official Gene Symbol │ IGHV1-2  │
│ 2   │ Bcl-2   │ Gene Symbol Alias    │ Official Gene Symbol │ BCL2     │

julia> convertFromTo(["IGHV1-2", "Bcl-2"], chromosomal = false)
WARNING:root:Input ID 'IGHV1-2' is an official gene symbol that codes for multiple Ensembl gene IDs. Possible reasons:
* One of the ENSG-IDs is located on an 'Alternative Human chromosome' or
* the ENSG-IDs encode different biotypes

Consider setting `chromosomal = false` and adding 'ArrayExpress' to targetDb parameter
2×4 DataFrames.DataFrame
│ Row │ InputId │ InputSourceDb        │ TargetDb             │ TargetId │
├─────┼─────────┼──────────────────────┼──────────────────────┼──────────┤
│ 1   │ IGHV1-2 │ Official Gene Symbol │ Official Gene Symbol │ IGHV1-2  │
│ 2   │ Bcl-2   │ Gene Symbol Alias    │ Official Gene Symbol │ BCL2     │

# Once we run this again we can see that IGHV1-2 also encodes for a ENSG-ID on an alternative sequence
julia> convertFromTo(["IGHV1-2", "Bcl-2"], ["EntrezGene", "GeneSymbolDB", "ArrayExpress"], chromosomal = false)
6×4 DataFrames.DataFrame
│ Row │ InputId │ InputSourceDb        │ TargetDb             │ TargetId        │
├─────┼─────────┼──────────────────────┼──────────────────────┼─────────────────┤
│ 1   │ IGHV1-2 │ Official Gene Symbol │ Ensembl Human Gene   │ ENSG00000211934 │
│ 2   │ IGHV1-2 │ Official Gene Symbol │ Official Gene Symbol │ IGHV1-2         │
│ 3   │ IGHV1-2 │ Official Gene Symbol │ Ensembl Human Gene   │ ENSG00000282550 │
│ 4   │ Bcl-2   │ Gene Symbol Alias    │ Ensembl Human Gene   │ ENSG00000171791 │
│ 5   │ Bcl-2   │ Gene Symbol Alias    │ Official Gene Symbol │ BCL2            │
│ 6   │ Bcl-2   │ Gene Symbol Alias    │ NCBI gene            │ 596             │

# Using the above query and setting `longFormat = false`:
julia> convertFromTo(["IGHV1-2", "Bcl-2"], ["EntrezGene", "GeneSymbolDB", "ArrayExpress"], chromosomal = false, longFormat=false)
3×5 DataFrames.DataFrame
│ Row │ InputId │ InputSourceDb        │ Ensembl_Human_Gene │ Official_Gene_Symbol │ NCBI_gene │
├─────┼─────────┼──────────────────────┼────────────────────┼──────────────────────┼───────────┤
│ 1   │ Bcl-2   │ Gene Symbol Alias    │ ENSG00000171791    │ BCL2                 │ 596       │
│ 2   │ IGHV1-2 │ Official Gene Symbol │ ENSG00000211934    │ IGHV1-2              │ missing   │
│ 3   │ IGHV1-2 │ Official Gene Symbol │ ENSG00000282550    │ IGHV1-2              │ missing   │

```

"""
function convertFromTo(
        queryId::S;
        targetDb::U = ["GeneSymbolDB"],
        sourceDb::S = "",
        longFormat::B = true,
        chromosomal::B = true) where {S<:String, U<:Union{String, Vector{String}}, B<:Bool}

return convertFromTo(
           [queryId],
           targetDb = isa(targetDb, String) ? [targetDb] : targetDb,
           sourceDb = sourceDb,
           longFormat = longFormat,
           chromosomal = chromosomal);

end

function convertFromTo(
        queryId::T;
        targetDb::U = ["GeneSymbolDB"],
        sourceDb::S = "",
        longFormat::B = true,
        chromosomal::B = true) where {S<:String, T<:Vector{String}, U<:Union{String, Vector{String}}, B<:Bool}

   targetDb = isa(targetDb, String) ? String[targetDb] : targetDb;

   tDb = copy(targetDb);
   append!(tDb, collect(keys(addDbColDict)));

   q = "CALL rcsi.convert.table.convertIds({ids}, {dbs}, " *
   (sourceDb == "" ? "null, " : ("'" * sourceDb * "', ")) *
   "{chromosomal}) " *
   "YIELD value " *
   "RETURN value.InputId AS InputId, " *
   "value.InputSourceDb AS InputSourceDb, " *
   "value.EnsemblGeneId AS EnsemblGeneId, " *
   "value.TargetDb AS TargetDb, " *
   "value.TargetId AS TargetId;"

   x = Neo4j.cypherQuery(
      NEO4J_CONNECTION,
      q,
      "ids" => queryId,
      "dbs" => targetDb,
      "chromosomal" => chromosomal);

   if DataFrames.nrow(x) > 0
      # Check if some symbols were not found or returned multiple values
      keysUq = unique(x[[:InputId, :InputSourceDb, :EnsemblGeneId]]);

      gd = DataFramesMeta.groupby(keysUq, :InputId)
      for group in gd
         DataFrames.nrow(group) == 1 && continue;
         g = unique(group[:InputSourceDb]);
         if length(g) == 1
            if g[1] == "Official Gene Symbol" && !in("ArrayExpress", targetDb)
               warn("Input ID '$(group[:InputId][1])' is an official gene symbol that codes for multiple Ensembl gene IDs. " *
                     "Possible reasons:" *
                     "\n* One of the ENSG-IDs is located on an 'Alternative Human chromosome' or " *
                     "\n* the ENSG-IDs encode different biotypes" *
                     "\n\nConsider " *
                     ((!chromosomal && !in("ArrayExpress", targetDb)) ?
                        "setting `chromosomal = false` and adding 'ArrayExpress' to targetDb parameter" :
                           (!chromosomal ? "setting `chromosomal = false`" : "adding 'ArrayExpress' to targetDb parameter")));
            elseif g[1] == "Gene Symbol Alias"
               warn("Input ID '$(group[:InputId][1])' represents multiple official gene symbols.")
            end
         elseif all(map(i->in(i,g),String["Gene Symbol Alias", "Official Gene Symbol"]))
            warn("Input ID '$(group[:InputId][1])' represents a gene symbol alias as well as an official gene symbol.")
         end
      end

      if !longFormat
         x = unstackDf(x);
      end

      delete!(x, :EnsemblGeneId);

      x = unique(x);

      checkReturnedQueryIds(x, queryId);
   end

   return x;

end

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


text:
"""
Module: SymbolicDiff (Symbolic Operation for Arithmetic)
"""

"""
seval(f, env, cache)
Return the value for expr f
"""

function seval(f)
    seval(f, globalenv, SymbolicCache())
end

function seval(f, env::SymbolicEnv)
    seval(f, env, SymbolicCache())
end

function seval(f, cache::SymbolicCache)
    seval(f, globalenv, cache)
end

function seval(f::SymbolicValue{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    f.val
end

function seval(f::SymbolicVariable{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    env[f.var]
end

function seval(f::AbstractNumberSymbolic{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    get(cache, f) do
        retval = _eval(Val(f.op), f, env, cache)
        cache[f] = retval
    end
end

"""
_eval(::Val{xx}, f, env, cache)

Dispached function to evaluate the expr f
"""

function _eval(::Val{:+}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, env, cache) for x = f.args]
    +(args...)
end

function _eval(::Val{:-}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, env, cache) for x = f.args]
    -(args...)
end

function _eval(::Val{:*}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, env, cache) for x = f.args]
    *(args...)
end

function _eval(::Val{:/}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    x/y
end

function _eval(::Val{:^}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    x^y
end

function _eval(::Val{:exp}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    exp(x)
end

function _eval(::Val{:log}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    @assert length(f.args) == 1
    x, = [seval(x, env, cache) for x = f.args]
    log(x)
end

function _eval(::Val{:sqrt}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    sqrt(x)
end

function _eval(::Val{:sum}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    sum(x)
end

function _eval(::Val{:dot}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    dot(x,y)
end

ext: jl


text:
"Simulate state trajectory with feedback control"
function rollout!(prob::Problem{T},solver::iLQRSolver{T},alpha::T=1.0) where T
    X = prob.X; U = prob.U
    K = solver.K; d = solver.d; X̄ = solver.X̄; Ū = solver.Ū
    X̄[1] = prob.x0

    for k = 2:prob.N
        # Calculate state trajectory difference
        δx = state_diff(X̄[k-1],X[k-1],prob,solver)

        # Calculate updated control
        Ū[k-1] = U[k-1] + K[k-1]*δx + alpha*d[k-1]

        # Propagate dynamics
        evaluate!(X̄[k], prob.model, X̄[k-1], Ū[k-1], get_dt(prob,Ū[k-1]))

        # Check that rollout has not diverged
        if ~(norm(X̄[k],Inf) < solver.opts.max_state_value && norm(Ū[k-1],Inf) < solver.opts.max_control_value)
            return false
        end
    end
    return true
end

function rollout!(prob::Problem{T}) where T
    N = prob.N
    if !all(isfinite.(to_array(prob.X)))
        prob.X[1] = prob.x0
        rollout!(prob.X, prob.model, prob.U, get_dt_traj(prob))
    end
end

function rollout!(X::AbstractVectorTrajectory, model::Model{M,Discrete}, U::AbstractVectorTrajectory, dt) where {M,T}
    N = length(X)
    for k = 1:N-1
        evaluate!(X[k+1], model, X[k], U[k], dt[k])
    end
end

function rollout(model::Model{M,Discrete}, x0::Vector, U::AbstractVectorTrajectory, dt) where M
    n = model.n
    N = length(U)+1
    X = [zero(x0) for k = 1:N]
    X[1] = x0
    rollout!(X, model, U, dt)
    return X
end
rollout(prob::Problem{T,Discrete}) where T = rollout(prob.model, prob.x0, prob.U, prob.dt)

function state_diff(x̄::Vector{T},x::Vector{T},prob::Problem{T},solver::iLQRSolver{T}) where T
    if true
        x̄ - x
    else
        nothing #TODO quaternion
    end
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using Test, Dates, DataFrames
import BLPData

if Sys.islinux() && haskey(ENV, "SKIP_LINUX_TESTS")
    @warn("Skipping tests for Linux.")
    exit()
end

@testset "Error Codes" begin
    error_string = BLPData.blpapi_getLastErrorDescription(Int32(1))
    @test error_string == "BLPAPI_ERROR_UNKNOWN"
end

@testset "Version Info" begin
    vinfo = BLPData.get_version_info()
    @test vinfo.major == 3
    println("BLPAPI Version: $vinfo")
end

@testset "BLPDateTime" begin
    today_dt = Dates.today()
    @test today_dt == Date(BLPData.BLPDateTime(today_dt))
end

@testset "Session Options" begin
    @testset "Default" begin
        opt = BLPData.SessionOptions()
        @test opt.handle != C_NULL
        server_host = BLPData.get_server_host(opt)
        server_port = BLPData.get_server_port(opt)
        client_mode = BLPData.get_client_mode(opt)
        println("Default Server Host: $server_host")
        println("Default Server Port: $server_port")
        println("Default client mode: $client_mode")
        BLPData.destroy!(opt)
        @test opt.handle == C_NULL
    end

    @testset "Constructor" begin
        opt = BLPData.SessionOptions(host="my_host", port=9999, client_mode=BLPData.BLPAPI_CLIENTMODE_DAPI)
        @test BLPData.get_server_host(opt) == "my_host"
        @test BLPData.get_server_port(opt) == 9999
        @test BLPData.get_client_mode(opt) == BLPData.BLPAPI_CLIENTMODE_DAPI
    end
end

@testset "BLPName" begin
    name1 = BLPData.BLPName(:HEY_YOU) # uses blpapi_Name_create
    name2 = BLPData.BLPName(:HEY_YOU) # uses blpapi_Name_create
    @test name1.handle == name2.handle

    # blpapi_Name_findName works only if blpapi_Name_create was used previously
    name3 = BLPData.BLPName(BLPData.blpapi_Name_findName(String(:HEY_YOU)))
    @test name3.handle == name1.handle

    BLPData.destroy!(name1)
    @test name1.handle == C_NULL

    name4 = BLPData.BLPName(BLPData.blpapi_Name_findName(String(:HEY_YOU)))
    @test name2.handle == name4.handle
end

@testset "CorrelationId" begin
    @testset "header" begin
        sz = UInt32(sizeof(BLPData.CorrelationId))
        tp = BLPData.BLPAPI_CORRELATION_TYPE_INT
        class = UInt32(10)
        reserved = UInt32(0)
        header = BLPData.CorrelationIdHeader(sz, tp, class, reserved)
        @test header == BLPData.CorrelationIdHeader(BLPData.encode_header(header))
    end

    @testset "corr" begin
        corr = BLPData.CorrelationId(10, class_id=UInt32(255))
        header = BLPData.CorrelationIdHeader(corr)
        @test header.struct_size == sizeof(BLPData.CorrelationId)
        @test header.correlation_type == BLPData.BLPAPI_CORRELATION_TYPE_INT
        @test header.class_id == 0xff
        @test header.reserved == 0
    end
end

@testset "SubscriptionTopic" begin
    @testset "eq default" begin
        topic1 = BLPData.SubscriptionTopic(BLPData.CorrelationId(), "hey you")
        topic2 = BLPData.SubscriptionTopic(BLPData.CorrelationId(), "hey you")
        @test topic1 == topic2
        @test hash(topic1) == hash(topic2)
    end

    @testset "eq custom value" begin
        topic1 = BLPData.SubscriptionTopic(BLPData.CorrelationId(5), "hey you")
        topic2 = BLPData.SubscriptionTopic(BLPData.CorrelationId(5), "hey you")
        @test topic1 == topic2
        @test hash(topic1) == hash(topic2)
    end

    @testset "diff custom value" begin
        topic1 = BLPData.SubscriptionTopic(BLPData.CorrelationId(4), "hey you")
        topic2 = BLPData.SubscriptionTopic(BLPData.CorrelationId(5), "hey you")
        @test topic1 != topic2
    end
end

@testset "SubscriptionList" begin
    sublist = BLPData.SubscriptionList()
    @test isempty(sublist)
    result = push!(sublist, "hey you out there")
    @test !isempty(sublist)
    @test length(sublist) == 1
    @test result == sublist[1]
    result2 = push!(sublist, "can you feel me")
    @test length(sublist) == 2

    for topic in sublist
        @test isa(topic, BLPData.SubscriptionTopic)
        @test topic == result || topic == result2
        @test topic.correlation_id != BLPData.CorrelationId()
    end
end

if haskey(ENV, "NO_BLP_SERVICE")
    @info("NO_BLP_SERVICE env variable is set. Skipping Terminal tests.")
else
    include("terminal_testsets.jl")
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 5
n = 5
k = 10
blockdim = 5

for (func,typ) in ((:test_geam_csr,:CudaSparseMatrixCSR),
                   (:test_geam_csc,:CudaSparseMatrixCSC))
    @eval begin
        function $func(elty)
            A = sparse(rand(elty,m,n))
            B = sparse(rand(elty,m,n))
            alpha = rand(elty)
            beta = rand(elty)
            C = alpha * A + beta * B
            d_A = $typ(A)
            d_B = $typ(B)
            d_C = CUSPARSE.geam(alpha,d_A,beta,d_B,'O','O','O')
            h_C = to_host(d_C)
            @test_approx_eq(C,h_C)
            d_C = CUSPARSE.geam(d_A,beta,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = A + beta * B
            @test_approx_eq(C,h_C)
            d_C = CUSPARSE.geam(alpha,d_A,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = alpha * A + B
            @test_approx_eq(C,h_C)
            d_C = CUSPARSE.geam(d_A,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = A + B
            @test_approx_eq(C,h_C)
            d_C = d_A + d_B
            h_C = to_host(d_C)
            C = A + B
            @test_approx_eq(C,h_C)
            d_C = d_A - d_B
            h_C = to_host(d_C)
            C = A - B
            @test_approx_eq(C,h_C)
            B = sparse(rand(elty,k,n))
            d_B = $typ(B)
            @test_throws(DimensionMismatch,CUSPARSE.geam(d_B,d_A,'O','O','O'))
        end
    end
end

function test_geam_mix(elty)
    A = spdiagm(rand(elty,m))
    B = spdiagm(rand(elty,m)) + sparse(diagm(rand(elty,m-1),1))
    d_A = CudaSparseMatrixCSR(A)
    d_B = CudaSparseMatrixCSC(B)
    d_C = d_B + d_A
    h_C = to_host(d_C)
    @test h_C ≈ A + B
    d_C = d_A + d_B
    h_C = to_host(d_C)
    @test h_C ≈ A + B
    d_C = d_A - d_B
    h_C = to_host(d_C)
    @test h_C ≈ A - B
    d_C = d_B - d_A
    h_C = to_host(d_C)
    @test h_C ≈ B - A
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_geam_csr(elty)
    test_geam_csc(elty)
    test_geam_mix(elty)
    println("geam took ", toq(), " for ", elty)
end


text:
"""
    struct CLagrangianFESpace{S} <: SingleFieldFESpace
      grid::Grid
      dof_to_node::Vector{Int}
      dof_to_comp::Vector{Int8}
      node_and_comp_to_dof::Vector{S}
      # + private fields
    end
"""
struct CLagrangianFESpace{S} <: SingleFieldFESpace
  space::UnconstrainedFESpace
  grid::Grid
  dof_to_node::Vector{Int}
  dof_to_comp::Vector{Int8}
  node_and_comp_to_dof::Vector{S}

  @doc """
      CLagrangianFESpace(::Type{T},grid::Grid) where T
  """
  function CLagrangianFESpace(::Type{T},grid::Grid) where T
    space, dof_to_node, dof_to_comp, node_and_comp_to_dof = _generate_clagrangian_fespace(T,grid)
    S = eltype(node_and_comp_to_dof)
    new{S}(space,grid,dof_to_node,dof_to_comp,node_and_comp_to_dof)
  end
end

# FESpace interface

constraint_style(::Type{<:CLagrangianFESpace}) = Val{false}()

function num_free_dofs(f::CLagrangianFESpace)
  num_free_dofs(f.space)
end

function get_cell_basis(f::CLagrangianFESpace)
  get_cell_basis(f.space)
end

function zero_free_values(f::CLagrangianFESpace)
  zero_free_values(f.space)
end

# SingleFieldFESpace interface

function get_cell_dofs(f::CLagrangianFESpace)
  get_cell_dofs(f.space)
end

function get_cell_dof_basis(f::CLagrangianFESpace)
  get_cell_dof_basis(f.space)
end

function num_dirichlet_dofs(f::CLagrangianFESpace)
  num_dirichlet_dofs(f.space)
end

function num_dirichlet_tags(f::CLagrangianFESpace)
  num_dirichlet_tags(f.space)
end

function zero_dirichlet_values(f::CLagrangianFESpace)
  zero_dirichlet_values(f.space)
end

function get_dirichlet_dof_tag(f::CLagrangianFESpace)
  get_dirichlet_dof_tag(f.space)
end

function scatter_free_and_dirichlet_values(f::CLagrangianFESpace,free_values,dirichlet_values)
  scatter_free_and_dirichlet_values(f.space,free_values,dirichlet_values)
end

function gather_free_and_dirichlet_values!(free_values,dirichlet_values,f::CLagrangianFESpace,cell_vals)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f.space,cell_vals)
end

function gather_dirichlet_values!(dirichlet_vals,f::CLagrangianFESpace,cell_vals)
    gather_dirichlet_values!(dirichlet_vals,f.space,cell_vals)
end

# Helpers

function _generate_clagrangian_fespace(T,grid)

    grid_reffes = get_reffes(grid)
    function fun(reffe)
      p = get_polytope(reffe)
      orders = get_orders(reffe)
      LagrangianRefFE(T,p,orders)
    end
    reffes = map(fun,grid_reffes)

    nnodes = num_nodes(grid)
    dof_to_node, dof_to_comp, node_and_comp_to_dof = _generate_dof_layout_component_major(T,nnodes)

    cell_nodes = get_cell_nodes(grid)
    cell_to_ctype = get_cell_type(grid)

    cell_dofs = _generate_cell_dofs_clagrangian_fespace(
      T,
      cell_nodes,
      reffes,
      cell_to_ctype,
      node_and_comp_to_dof)

    cell_map = get_cell_map(grid)
    cell_shapefuns, cell_dof_basis = compute_cell_space(reffes, cell_to_ctype, cell_map)

    nfree = length(dof_to_node)
    ndirichlet = 0

    dirichlet_dof_tag = Int8[]
    dirichlet_cells = Int[]
    ntags = 0

    space = UnconstrainedFESpace(
      nfree,
      ndirichlet,
      cell_dofs,
      cell_shapefuns,
      cell_dof_basis,
      dirichlet_dof_tag,
      dirichlet_cells,
      ntags)

    space, dof_to_node, dof_to_comp, node_and_comp_to_dof

end

function _generate_dof_layout_component_major(::Type{<:Real},nnodes::Integer)
  ndofs = nnodes
  dof_to_comp = ones(Int8,ndofs)
  dof_to_node = collect(1:nnodes)
  node_and_comp_to_dof = dof_to_node
  (dof_to_node, dof_to_comp, node_and_comp_to_dof)
end

function _generate_dof_layout_component_major(::Type{T},nnodes::Integer) where T
  ncomps = num_components(T)
  V = change_eltype(T,Int)
  ndofs = ncomps*nnodes
  dof_to_comp = zeros(Int8,ndofs)
  dof_to_node = zeros(Int,ndofs)
  node_and_comp_to_dof = zeros(V,nnodes)
  m = zero(mutable(V))
  for node in 1:nnodes
    o = ncomps*(node-1)
    for comp in 1:ncomps
      dof = comp+o
      dof_to_comp[dof] = comp
      dof_to_node[dof] = node
      m[comp] = dof
    end
    node_and_comp_to_dof[node] = m
  end
  (dof_to_node, dof_to_comp, node_and_comp_to_dof)
end

function _generate_cell_dofs_clagrangian_fespace(
  ::Type{<:Real},
  cell_nodes,
  reffes,
  cell_to_ctype,
  node_and_comp_to_dof)

  cell_nodes
end

function _generate_cell_dofs_clagrangian_fespace(
  ::Type{T},
  cell_nodes,
  reffes,
  cell_to_ctype,
  node_and_comp_to_dof) where T

  ncomps = num_components(T)

  ctype_to_lnode_to_comp_to_ldof = map(get_node_and_comp_to_dof,reffes)
  ctype_to_num_ldofs = map(num_dofs,reffes)

  cell_dofs = _allocate_cell_dofs_clagrangian_fespace(ctype_to_num_ldofs,cell_to_ctype)

  cache = array_cache(cell_nodes)
  _fill_cell_dofs_clagrangian_fespace!(
    cell_dofs,
    cache,
    cell_nodes,
    ctype_to_lnode_to_comp_to_ldof,
    cell_to_ctype,
    node_and_comp_to_dof,
    ncomps)

  cell_dofs
end

function _allocate_cell_dofs_clagrangian_fespace(ctype_to_num_ldofs,cell_to_ctype)
  ncells = length(cell_to_ctype)
  ptrs = zeros(Int32,ncells+1)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    nldofs = ctype_to_num_ldofs[ctype]
    ptrs[cell+1] = nldofs
  end
  length_to_ptrs!(ptrs)
  ndata = ptrs[end]-1
  data = zeros(Int,ndata)
  Table(data,ptrs)
end

function _fill_cell_dofs_clagrangian_fespace!(
  cell_dofs,
  cache,
  cell_nodes,
  ctype_to_lnode_and_comp_to_ldof,
  cell_to_ctype,
  node_and_comp_to_dof,
  ncomps)

  ncells = length(cell_to_ctype)

  for cell in 1:ncells
    nodes = getindex!(cache,cell_nodes,cell)
    ctype = cell_to_ctype[cell]
    lnode_and_comp_to_ldof = ctype_to_lnode_and_comp_to_ldof[ctype]
    p = cell_dofs.ptrs[cell]-1
    for (lnode, node) in enumerate(nodes)
      for comp in 1:ncomps
        ldof = lnode_and_comp_to_ldof[lnode][comp]
        dof = node_and_comp_to_dof[node][comp]
        cell_dofs.data[p+ldof] = dof
      end
    end
  end

end


text:
"""
    t_matrix(particle, medium, ω, order)

Returns a finite T-matrix, with size depending on `order`, for a specific `particle` within a `medium` with specific physical properties.
"""
function t_matrix(p::AbstractParticle{T,Dim}, medium::PhysicalMedium{T,Dim}, ω::T, order::Integer)::AbstractMatrix{T} where {T<:AbstractFloat,Dim}

    error("T-matrix function is not yet written for $(name(p.medium)) $(name(p.shape)) in a $(name(medium)) medium")
end

"""
    get_t_matrices(PhysicalMedium, Vector{Particles}, ω, basis_order::Integer)

Returns vector of T-matrices from a vector of particles in a specific domain.
Can save computation if multiple of the same kind of particle are present in the
vector.
"""
function get_t_matrices(medium::PhysicalMedium, particles::AbstractParticles, ω::AbstractFloat, basis_order::Integer)::Vector

    t_matrices = Vector{AbstractMatrix}(undef,length(particles))

    # Vector of particles unique up to congruence, and the respective T-matrices
    unique_particles = Vector{AbstractParticle}(undef,0)
    unique_t_matrices = Vector{AbstractMatrix}(undef,0)

    for p_i in eachindex(particles)
        p = particles[p_i]

        # If we have calculated this T-matrix before, just point to that one
        found = false
        for cp_i in eachindex(unique_particles)
            cp = unique_particles[cp_i]
            if iscongruent(p, cp)
                t_matrices[p_i] = unique_t_matrices[cp_i]
                found = true
                break
            end
        end

        # Congruent particle was not found, we must calculate this t-matrix
        if !found
            t_matrices[p_i] = t_matrix(p, medium, ω, basis_order)
            push!(unique_particles, particles[p_i])
            push!(unique_t_matrices, t_matrices[p_i])
        end

    end

    return t_matrices
end


text:
using WriteVTK

read_triangles(; input::AbstractString, output::AbstractString) = open(input) do f
    num_nodes = read(f, Int)
    num_elements = read(f, Int)
    coords = Matrix{Float64}(undef, 2, num_nodes)
    read!(f, coords)
    vertices = Matrix{Int64}(undef, 3, num_elements)
    read!(f, vertices)
    vertices .+= 1

    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, vertices[:, j])
             for j in 1:num_elements]

    vtk_grid(output, coords, cells)
end

max_stars_repo_name: slmcbane/Elasticity

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: slmcbane/Elasticity

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: slmcbane/Elasticity

max_forks_repo_licenses: MIT Unlicense


text:
test = """
    ..#.#..#####.#.#.#.###.##.....###.##.#..###.####..#####..#....#..#..##..###..######.###...####..#..#####..##..#.#####...##.#.#..#.##..#.#......#.###.######.###.####...#.##.##..#..#..#####.....#.#....###..#.##......#.....#..#..#..##..#...##.######.####.####.#.#...#.......#..#.#.#...####.##.#......#..#...##.#.##..#...##.#.##..###.#......#.#.......#.#.#.####.###.##...#.....####.#..#..#.##.#....##..#.####....##...##..#...#......#.#.......#.......##..####..#...#.#.#...##..#.#..###..#####........#..####......#..#

    #..#.
    #....
    ##..#
    ..#..
    ..###
    """

function parse_input(io::IO)
    cnv(c) = c == "#" ? 1 : 0
    alg = cnv.(split(readline(io), ""))
    @assert length(alg) == 512
    @assert alg[1] & alg[end] == 0
    @assert isempty(readline(io))
    lines = [cnv.(split(line, "")) for line in readlines(io)]
    image = permutedims(hcat(lines...))
    return (alg, pad(image))
end

function pad(M::AbstractMatrix; padding::Integer = 2)::AbstractMatrix
    (h, w) = size(M)
    return [fill(0, padding, w + 2padding)
            fill(0, h, padding) M fill(0, h, padding)
            fill(0, padding, w + 2padding)]
end

function str(M::AbstractMatrix)::String
    return join([join(map(v -> v == 1 ? '#' : '.', row))
                 for row in eachrow(M)], "\n")
end

function decode(M::AbstractMatrix, ci::CartesianIndex)::Int
    bits = vec(permutedims(M[ci - one(ci):ci + one(ci)]))
    return parse(Int, join(bits, ""); base = 2) + 1
end

function enhance(alg::AbstractVector, M::AbstractMatrix)::AbstractMatrix
    R = fill(M[1] == 0 ? alg[1] : alg[end], size(M) .+ (2, 2))
    for ci in CartesianIndex(2, 2):CartesianIndex(size(M) .- (1, 1))
        R[ci + one(ci)] = alg[decode(M, ci)]
    end
    return clip(R)
end

function clip(M::AbstractMatrix; padding::Integer = 2)::AbstractMatrix
    (h, w) = size(M)
    foreground = 1 - M[1]
    (top, bottom) = extrema([r for r in 1:h if foreground ∈ M[r, :]])
    (left, right) = extrema([c for c in 1:w if foreground ∈ M[:, c]])
    return M[top - padding:bottom + padding, left - padding:right + padding]
end

function q1(alg, img)
    return sum(enhance(alg, enhance(alg, img)))
end

function q2(alg, img)
    return sum(reduce((img, _) -> enhance(alg, img), 1:50; init = img))
end

let input = parse_input(IOBuffer(test))
    @assert q1(input...) == 35
    @assert q2(input...) == 3351
end

open("day20.in", "r") do io
    input = parse_input(io)
    println("Q1: ", q1(input...))
    println("Q2: ", q2(input...))
end

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
Random.seed!(0)

# Linear Dynamic model
# = x = [p; v]
dt = 0.01
A = [1.0 dt; 0.0 1.0]
B = zeros(2, 2)
W = 0.01*Matrix{Float64}(I, 2, 2)

dyn = LinearDynamicModel(A, B, W)

# Linear Observation Model
# observe postition
H = [1 0]
obs = LinearObservationModel(H)

# define linear kalman filter
kf = KalmanFilter(dyn, obs)

# simulation
R_gt = 0.5*Matrix{Float64}(I, 1, 1) # ground truth noise covariance
time_step = 0.1:dt:1000
s0 = State([0.0, 1.0], Matrix{Float64}(I, 2, 2)) # initial state belief


text:
"""
    create_HO(dim::Int)

This function creates the operators for a Harmonic oscillator of the specified dimension.

## args
* dim: Hilbert space dimesnion of the oscillator

## returns
(lower, raise, N): tuple of the lowering, raising, and photon number operator for the HO
"""
function create_HO(dim::Int)

    lower = zeros(dim,dim)
    for ii = 1:(dim-1)
        @inbounds lower[ii,ii+1] = sqrt(ii)
    end

    raise = lower'

    N = raise*lower

    return lower, raise, N
end

"""
    create_HO_network(dim::Array{Int},freqs::Vector{Float64},couple::Array{T,2}) where {T <: Number}

This function creates the Hamiltonian for a Harmonic oscillator network.

## args
* dim:      array of the Hilbert space dimensions of the oscillators
* freqs:    frequencies of the oscillators
* couple:   matrix of tunnel couplings (beamsplitter interaction) between the HOs. Only the upper triangle is used, and the diagonal should be zero. For row index j and column index k, couple[j,k] is the rate for the operator raise_j ⊗ lower_k

## returns
* H:        full Hamiltonian of the network
"""
function create_HO_network(dim::Array{Int},freqs::Vector{Float64},couple::Array{T,2}) where {T <: Number}

    num = length(freqs)

    H = zeros(prod(dim),prod(dim))

    # Tunnel coupling (note the normal ordering)
    for jj = 1:1:num
        ~, raise1, ~ = create_HO(dim[jj])
        for kk = (jj+1):1:num
            if abs(couple[jj,kk]) > 1e-10
                lower2, ~, ~ = create_HO(dim[kk])
                # raise_jj = kron(Matrix{ComplexF64}(I,dim^(jj-1),dim^(jj-1)),kron(raise,Matrix{ComplexF64}(I,dim^(num-jj),dim^(num-jj))))
                raise_jj = kron(Matrix{ComplexF64}(I,prod(dim[1:jj-1]),prod(dim[1:jj-1])),kron(raise1,Matrix{ComplexF64}(I,prod(dim[(jj+1):end]),prod(dim[(jj+1):end]))))
                # lower_kk = kron(Matrix{ComplexF64}(I,dim^(kk-1),dim^(kk-1)),kron(lower,Matrix{ComplexF64}(I,dim^(num-kk),dim^(num-kk))))
                lower_kk = kron(Matrix{ComplexF64}(I,prod(dim[1:kk-1]),prod(dim[1:kk-1])),kron(lower2,Matrix{ComplexF64}(I,prod(dim[(kk+1):end]),prod(dim[(kk+1):end]))))

                @inbounds H = H + couple[jj,kk]*raise_jj*lower_kk
            end
        end
    end

    H = H + H'

    # Self Hamiltonians
    for ii = 1:1:num
        ~, ~, N1 = create_HO(dim[ii])
        @inbounds H = H + freqs[ii]*kron(Matrix{ComplexF64}(I,prod(dim[1:ii-1]),prod(dim[1:ii-1])),kron(N1,Matrix{ComplexF64}(I,prod(dim[(ii+1):end]),prod(dim[(ii+1):end]))))
    end

    return H

end

# To be updated
# In place version of the above
"""
    create_HO_network!(H::Array{T1,2},dim::Int,freqs::Vector{Float64},couple::Array{T2,2}) where {T1 <: Number, T2 <: Number}

In place version of `create_HO_network` where the Hamiltonian H is modified

## args
* H:      full Hamiltonian of the network
* dim:    array of the Hilbert space dimensions of the oscillators
* freqs:  frequencies of the oscillators
* couple: matrix of tunnel couplings (beamsplitter interaction) between the
            HOs. Only the upper triangle is used, and the diagonal should be zero.
            For row index j and column index k, couple[j,k] is the rate for the
            operator raise_j ⊗ lower_k
## returns
nothing
"""
function create_HO_network!(H::Array{T1,2},dim::Int,freqs::Vector{Float64},couple::Array{T2,2}) where {T1 <: Number, T2 <: Number}

    num = length(freqs)
    (lower, raise, N) = create_HO(dim)

    # Tunnel coupling (note the normal ordering)
    for jj = 1:1:num
        for kk = (jj+1):1:num
            if abs(couple[jj,kk]) > 1e-10
                raise_jj = kron(Matrix{ComplexF64}(I,dim^(jj-1),dim^(jj-1)),kron(raise,Matrix{ComplexF64}(I,dim^(num-jj),dim^(num-jj))))
                lower_kk = kron(Matrix{ComplexF64}(I,dim^(kk-1),dim^(kk-1)),kron(lower,Matrix{ComplexF64}(I,dim^(num-kk),dim^(num-kk))))

                H[:,:] = H[:,:] + couple[jj,kk]*raise_jj*lower_kk
            end
        end
    end

    H[:,:] = H[:,:] + H[:,:]'

    # Self Hamiltonians
    for ii = 1:1:num
        H[:,:] = H[:,:] + freqs[ii]*kron(Matrix{ComplexF64}(I,dim^(ii-1),dim^(ii-1)),kron(N,Matrix{ComplexF64}(I,dim^(num-ii),dim^(num-ii))))
    end

    nothing

end

"""
    create_KPO_network(dim::Array{Int},freqs::Vector{Float64},couple::Array{T1,2},sqz::Array{T2,2},kerr::Array{Float64,2}) where {T1 <: Number, T2 <: Number}

This function creates the Hamiltonian for a Kerr parametric oscillator network.

## args
* dim:    array of the Hilbert space dimensions of the oscillators
* freqs:  frequencies of the oscillators
* couple: matrix of tunnel couplings (beamsplitter interaction) between the
            HOs. Only the upper triangle is used, and the diagonal should be zero.
            For row index j and column index k, couple[j,k] is the rate for the
            operator raise_j ⊗ lower_k
* sqz:    matrix of squeezing interaction rates. Only the upper triangle is
            used, with the diagonal describing single-mode, and the off-diagonal
            two-mode squeezing. For row index j and column index k, sqz[j,k] is
            the rate for the operator raise_j ⊗ raise_k
* kerr:   matrix of the self and cross Kerr interaction rates. Only the upper
            triangle is used with the diagonal describing self, and the
            off-diagonal cross Kerr. For row index j and column index k, kerr[j,k]
            is the rate for the operator N_j ⊗ N_k. Note that this is NOT
            NORMAL ORDERED for j = k

## returns
* H:      full Hamiltonian of the network
"""
function create_KPO_network(dim::Array{Int},freqs::Vector{Float64},couple::Array{T1,2},sqz::Array{T2,2},kerr::Array{Float64,2}) where {T1 <: Number, T2 <: Number}

    num = length(freqs)
    # (lower, raise, N) = create_HO(dim)

    H = zeros(prod(dim),prod(dim))

    # Squeezing and Kerr (factor of 1/2 in Kerr term is to account for doing H + H' later)
    for jj = 1:1:num
        ~, raise1, N1 = create_HO(dim[jj])
        for kk = jj:1:num
            if abs(sqz[jj,kk]) > 1e-10
                ~, raise2, ~ = create_HO(dim[kk])
                # raise_jj = kron(Matrix{ComplexF64}(I,dim^(jj-1),dim^(jj-1)),kron(raise,Matrix{ComplexF64}(I,dim^(num-jj),dim^(num-jj))))
                raise_jj = kron(Matrix{ComplexF64}(I,prod(dim[1:jj-1]),prod(dim[1:jj-1])),kron(raise1,Matrix{ComplexF64}(I,prod(dim[(jj+1):end]),prod(dim[(jj+1):end]))))
                # raise_kk = kron(Matrix{ComplexF64}(I,dim^(kk-1),dim^(kk-1)),kron(raise,Matrix{ComplexF64}(I,dim^(num-kk),dim^(num-kk))))
                raise_kk = kron(Matrix{ComplexF64}(I,prod(dim[1:kk-1]),prod(dim[1:kk-1])),kron(raise2,Matrix{ComplexF64}(I,prod(dim[(kk+1):end]),prod(dim[(kk+1):end]))))

                H = H + sqz[jj,kk]*raise_jj*raise_kk
            end

            if abs(kerr[jj,kk]) > 1e-10

                ~, ~, N2 = create_HO(dim[kk])
                N_jj = kron(Matrix{ComplexF64}(I,prod(dim[1:jj-1]),prod(dim[1:jj-1])),kron(N1,Matrix{ComplexF64}(I,prod(dim[(jj+1):end]),prod(dim[(jj+1):end]))))
                N_kk = kron(Matrix{ComplexF64}(I,prod(dim[1:kk-1]),prod(dim[1:kk-1])),kron(N2,Matrix{ComplexF64}(I,prod(dim[(kk+1):end]),prod(dim[(kk+1):end]))))

                H = H + kerr[jj,kk]*N_jj*N_kk/2  # This expression is NOT NORMAL ORDERED for jj = kk
            end
        end
    end

    H = H + H' + create_HO_network(dim,freqs,couple)

    return H

end

# To be updated
# In place version of the above
"""
    create_KPO_network!(H::Array{T1,2},dim::Int,freqs::Vector{Float64},couple::Array{T2,2},sqz::Array{T3,2},kerr::Array{Float64,2}) where {T1 <: Number, T2 <: Number, T3 <: Number}

In place version of `create_KPO_network` where the Hamiltonian H is update.

This function creates the Hamiltonian for a Kerr parametric oscillator network.

## args
* dim:    array of the Hilbert space dimensions of the oscillators
* freqs:  frequencies of the oscillators
* couple: matrix of tunnel couplings (beamsplitter interaction) between the
            HOs. Only the upper triangle is used, and the diagonal should be zero.
            For row index j and column index k, couple[j,k] is the rate for the
            operator raise_j ⊗ lower_k
* sqz:    matrix of squeezing interaction rates. Only the upper triangle is
            used, with the diagonal describing single-mode, and the off-diagonal
            two-mode squeezing. For row index j and column index k, sqz[j,k] is
            the rate for the operator raise_j ⊗ raise_k
* kerr:   matrix of the self and cross Kerr interaction rates. Only the upper
            triangle is used with the diagonal describing self, and the
            off-diagonal cross Kerr. For row index j and column index k, kerr[j,k]
            is the rate for the operator N_j ⊗ N_k. Note that this is NOT
            NORMAL ORDERED for j = k

## returns
nothing
"""
function create_KPO_network!(H::Array{T1,2},dim::Int,freqs::Vector{Float64},couple::Array{T2,2},sqz::Array{T3,2},kerr::Array{Float64,2}) where {T1 <: Number, T2 <: Number, T3 <: Number}

    num = length(freqs)
    (lower, raise, N) = create_HO(dim)

    # Squeezing and Kerr (factor of 1/2 in Kerr term is to account for doing H + H' later)
    for jj = 1:1:num
        for kk = jj:1:num
            if abs(sqz[jj,kk]) > 1e-10
                raise_jj = kron(Matrix{ComplexF64}(I,dim^(jj-1),dim^(jj-1)),kron(raise,Matrix{ComplexF64}(I,dim^(num-jj),dim^(num-jj))))
                raise_kk = kron(Matrix{ComplexF64}(I,dim^(kk-1),dim^(kk-1)),kron(raise,Matrix{ComplexF64}(I,dim^(num-kk),dim^(num-kk))))

                N_jj = kron(Matrix{ComplexF64}(I,dim^(jj-1),dim^(jj-1)),kron(N,Matrix{ComplexF64}(I,dim^(num-jj),dim^(num-jj))))
                N_kk = kron(Matrix{ComplexF64}(I,dim^(kk-1),dim^(kk-1)),kron(N,Matrix{ComplexF64}(I,dim^(num-kk),dim^(num-kk))))

                H[:,:] = H[:,:] + sqz[jj,kk]*raise_jj*raise_kk
            end

            if abs(kerr[jj,kk]) > 1e-10

                N_jj = kron(Matrix{ComplexF64}(I,dim^(jj-1),dim^(jj-1)),kron(N,Matrix{ComplexF64}(I,dim^(num-jj),dim^(num-jj))))
                N_kk = kron(Matrix{ComplexF64}(I,dim^(kk-1),dim^(kk-1)),kron(N,Matrix{ComplexF64}(I,dim^(num-kk),dim^(num-kk))))

                H[:,:] = H[:,:] + kerr[jj,kk]*N_jj*N_kk/2  # This expression is NOT NORMAL ORDERED for jj = kk
            end
        end
    end

    H[:,:] = H[:,:] + H[:,:]' + create_HO_network(dim,freqs,couple)

    nothing

end

"""
    create_Qubit()

This function creates the operators for a single qubit. Note the convention
used that (1,0)^T is the excited state.

## returns
* X, Y, Z:  Pauli matrices
* lower:    lowering operator
* raise:    raising operator

as a tuple : (X, Y, Z, lower, raise)
"""
function create_Qubit()

    Z = [1. 0.;0. -1.]
    X = [0. 1.;1. 0.]
    Y = -1im*Z*X
    lower = (X - 1im*Y)/2
    raise = lower'

    return X, Y, Z, lower, raise

end

"""
    create_Qubit(conv::Int)

This function creates the operators for a single qubit. Note the convention
used that (1,0)^T is the excited state.

## args
* conv:     +1 or -1 indicating the convention that either (1,0)^T or (0,1)^T is
                the excited state. +1 and -1 imply a qubit self Hamiltonian of
                +Z or -Z respectively.

## returns:
* X, Y, Z:  Pauli matrices
* lower:    lowering operator
* raise:    raising operator

as a tuple : (X, Y, Z, lower, raise)
"""
function create_Qubit(conv::Int)

    Z = [1. 0.;0. -1.]
    X = [0. 1.;1. 0.]
    Y = -1im*Z*X
    lower = (X - conv*1im*Y)/2
    raise = lower'

    return X, Y, Z, lower, raise

end

"""
    create_qubit_network(freqs::Vector{Float64},couple::Array{T,2};Ctype::String = "flipflop",conv::Int = 1) where {T <: Number}

This function creates a network of qubits with the coupling type specified and returns its Hamiltonian.
Input arguments:

## args
* freqs:    vector of frequencies of the qubits
* couple:   matrix of coupling rates between the qubits. Only the upper triangle
                is used, and the diagonal should be zero. For row index j and
                column index k, couple[j,k] is the rate for the operator
                raise_j ⊗ lower_k for the (default) RWA flipflop interaction
* Ctype:    one of {"XX","YY","ZZ","flipflop"} describing the type of the
                interaction. The first three possibilites describe Pauli
                matrix coupling between the qubits. If type = "flipflop" or is
                unset, then the interaction is the RWA, excitation number
                conserving flip-flop interaction raise_j ⊗ lower_k + h.c..
                Note that the type is consisnent across the qubit network.
* conv:     +1 or -1 indicating the convention that either (1,0)^T or (0,1)^T
                is the excited state. +1 and -1 imply a qubit self Hamiltonian
                of +Z or -Z respectively.

## returns
* H:        full Hamiltonian of the network
"""
function create_qubit_network(freqs::Vector{Float64},couple::Array{T,2};Ctype::String = "flipflop",conv::Int = 1) where {T <: Number}

    num = length(freqs)
    X, Y, Z, lower, raise = create_Qubit(conv)

    H = zeros(2^num,2^num)

    # Coupling
    if Ctype == "XX"
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10 && imag(couple[jj,kk]) < 1e-10
                    # X_jj = kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(X,Matrix{ComplexF64}(I,2^(num-jj),2^(num-jj))))
                    # X_kk = kron(Matrix{ComplexF64}(I,2^(kk-1),2^(kk-1)),kron(X,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))
                    # H = H + couple[jj,kk]*X_jj*X_kk
                    H = H + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(X,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(X,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end

    elseif Ctype == "YY"
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10 && imag(couple[jj,kk]) < 1e-10
                    # Y_jj = kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(Y,Matrix{ComplexF64}(I,2^(num-jj),2^(num-jj))))
                    # Y_kk = kron(Matrix{ComplexF64}(I,2^(kk-1),2^(kk-1)),kron(Y,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))
                    # H = H + couple[jj,kk]*Y_jj*Y_kk
                    H = H + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(Y,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(Y,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end

    elseif Ctype == "ZZ"
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10 && imag(couple[jj,kk]) < 1e-10
                    # Y_jj = kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(Y,Matrix{ComplexF64}(I,2^(num-jj),2^(num-jj))))
                    # Y_kk = kron(Matrix{ComplexF64}(I,2^(kk-1),2^(kk-1)),kron(Y,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))
                    # H = H + couple[jj,kk]*Y_jj*Y_kk
                    H = H + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(Z,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(Z,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end

    else
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10
                    # raise_jj = kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(raise,Matrix{ComplexF64}(I,2^(num-jj),2^(num-jj))))
                    # lower_kk = kron(Matrix{ComplexF64}(I,2^(kk-1),2^(kk-1)),kron(lower,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))
                    # H = H + couple[jj,kk]*raise_jj*lower_kk
                    H = H + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(raise,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(lower,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end
        H = H + H'
    end

    # Self Hamiltonians
    for ii = 1:1:num
        H = H + conv*freqs[ii]*kron(Matrix{ComplexF64}(I,2^(ii-1),2^(ii-1)),kron(Z,Matrix{ComplexF64}(I,2^(num-ii),2^(num-ii))))
    end

    return H
end

# In place version of above
"""
    create_qubit_network!(H::Array{T1,2},freqs::Vector{Float64},couple::Array{T2,2};Ctype::String = "flipflop",conv::Int = 1) where {T1 <: Number, T2 <: Number}

In place version of `create_Qubit_network` where the Hamiltonian H is updated.

## args
* freqs:    vector of frequencies of the qubits
* couple:   matrix of coupling rates between the qubits. Only the upper
                triangle is used, and the diagonal should be zero. For
                row index j and column index k, couple[j,k] is the rate
                for the operator raise_j ⊗ lower_k for the (default) RWA
                flipflop interaction
* Ctype:    one of {"XX","YY","ZZ","flipflop"} describing the type of the
                interaction. The first three possibilites describe Pauli
                matrix coupling between the qubits. If type = "flipflop"
                or is unset, then the interaction is the RWA, excitation
                number conserving flip-flop interaction
                raise_j ⊗ lower_k + h.c.. Note that the type is consisnent
                across the qubit network.
* conv:     +1 or -1 indicating the convention that either (1,0)^T or
                (0,1)^T is the excited state. +1 and -1 imply a qubit self
                Hamiltonian of +Z or -Z respectively.

## returns
* nothing
"""
function create_qubit_network!(H::Array{T1,2},freqs::Vector{Float64},couple::Array{T2,2};Ctype::String = "flipflop",conv::Int = 1) where {T1 <: Number, T2 <: Number}

    num = length(freqs)
    X, Y, Z, lower, raise = create_Qubit(conv)

    # Coupling
    if Ctype == "XX"
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10 && imag(couple[jj,kk]) < 1e-10
                    H[:,:] = H[:,:] + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(X,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(X,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end

    elseif Ctype == "YY"
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10 && imag(couple[jj,kk]) < 1e-10
                    H[:,:] = H[:,:] + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(Y,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(Y,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end

    elseif Ctype == "ZZ"
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10 && imag(couple[jj,kk]) < 1e-10
                    H[:,:] = H[:,:] + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(Z,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(Z,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end

    else
        for jj = 1:1:num
            for kk = (jj+1):1:num
                if abs(couple[jj,kk]) > 1e-10
                    H[:,:] = H[:,:] + couple[jj,kk]*kron(Matrix{ComplexF64}(I,2^(jj-1),2^(jj-1)),kron(raise,kron(Matrix{ComplexF64}(I,2^(kk-jj-1),2^(kk-jj-1)),kron(lower,Matrix{ComplexF64}(I,2^(num-kk),2^(num-kk))))))
                end
            end
        end
        H[:,:] = H[:,:] + H[:,:]'
    end

    # Self Hamiltonians
    for ii = 1:1:num
        H[:,:] = H[:,:] + conv*freqs[ii]*kron(Matrix{ComplexF64}(I,2^(ii-1),2^(ii-1)),kron(Z,Matrix{ComplexF64}(I,2^(num-ii),2^(num-ii))))
    end

    nothing
end


text:
function DAMP(ALPHA)
      A=[-0.2670    0.8820   -0.1080   -8.8000   -0.1260   -0.3600   -7.2100   -0.3800    0.0610
         -0.1100    0.8520   -0.1080  -25.8000   -0.0260   -0.3590   -0.5400   -0.3630    0.0520
          0.3080    0.8760   -0.1880  -28.9000    0.0630   -0.4430   -5.2300   -0.3780    0.0520
          1.3400    0.9580    0.1100  -31.4000    0.1130   -0.4200   -5.2600   -0.3860   -0.0120
          2.0800    0.9620    0.2580  -31.2000    0.2080   -0.3830   -6.1100   -0.3700   -0.0130
          2.9100    0.9740    0.2260  -30.7000    0.2300   -0.3750   -6.6400   -0.4530   -0.0240
          2.7600    0.8190    0.3440  -27.7000    0.3190   -0.3290   -5.6900   -0.5500    0.0500
          2.0500    0.4830    0.3620  -28.2000    0.4370   -0.2940   -6.0000   -0.5820    0.1500
          1.5000    0.5900    0.6110  -29.0000    0.6800   -0.2300   -6.2000   -0.5950    0.1300
          1.4900    1.2100    0.5290  -29.8000    0.1000   -0.2100   -6.4000   -0.6370    0.1580
          1.8300   -0.4930    0.2980  -38.3000    0.4470   -0.1200   -6.6000   -1.0200    0.2400
          1.2100   -1.0400   -2.2700  -35.3000   -0.3300   -0.1000   -6.0000   -0.8400    0.1500]
      #
      S= 0.2 * ALPHA
      K= floor(S)+3
      if(K <= 1)
            K= 2
      end
      if(K >=  12)
            K=  11
      end
      DA= S - K +3
      L = K + fix( 1.1*sign(DA) )
      #
      D=zeros(9,1)
      K=convert(Int64, K)
      L=convert(Int64, L)
      for I= 1:9
            D[I]= A[K,I] + abs(DA) * (A[L,I] - A[K,I])
      end
      return D
end


text:
#=
This page contains various utility functions.
=#

########################################
# Size of layer - overload Base.size
########################################

size(ℓ::NetLayer) = size(ℓ.A)

########################################
# Network Display
########################################

# Display of the NeuralNet type
function show(io::IO, net::NeuralNet)
    println(io, "Neural Network")
    println(io, "Training Data Dimensions - $(size(net.data.X_train))")
    println(io, "Layers:")
    num_layers = length(net.layers)
    for (n, ℓ) in enumerate(net.layers)
        println(io, "Layer $n - ", displaystr(ℓ))
    end
end

# Display of hidden layers
function show(io::IO, ℓ::NetLayer)
    println(io, displaystr(ℓ))
end

# Display string for layers with activation functions
function displaystr(ℓ::NetLayer)
    str = string(typeof(ℓ))[7:end]
    activation = "$(ℓ.f!)"[7:end-1]
    str *= ", Activation - $activation"
    str *= ", Dimensions - $(size(ℓ))"
end

# Display string for layers without activation functions
function displaystr{L <: Union{InputLayer, WordEmbeddingLayer, PoolLayer, OutputLayer}}(ℓ::L)
    str = string(typeof(ℓ))[7:end]
    str *= ", Dimensions - $(size(ℓ))"
end

########################################
# Batch Views
########################################

# View of a batch - outer dimension contains the observations (e.g. images)
viewbatch{R<:Real, T<:Integer}(X::Array{R,1}, ids::Vector{T}) = view(X, ids)
viewbatch{R<:Real, T<:Integer}(X::Array{R,2}, ids::Vector{T}) = view(X, :, ids)
viewbatch{R<:Real, T<:Integer}(X::Array{R,3}, ids::Vector{T}) = view(X, :, :, ids)
viewbatch{R<:Real, T<:Integer}(X::Array{R,4}, ids::Vector{T}) = view(X, :, :, :, ids)
viewbatch{R<:Real, T<:Integer}(X::Array{R,1}, start::T, stop::T) = view(X, start:stop)
viewbatch{R<:Real, T<:Integer}(X::Array{R,2}, start::T, stop::T) = view(X, :, start:stop)
viewbatch{R<:Real, T<:Integer}(X::Array{R,3}, start::T, stop::T) = view(X, :, :, start:stop)
viewbatch{R<:Real, T<:Integer}(X::Array{R,4}, start::T, stop::T) = view(X, :, :, :, start:stop)

########################################
# Count all params in net
########################################

function num_params(net::NeuralNet)
    count = 0
    for ℓ in net.layers[2:end]
        typeof(ℓ) <: PoolLayer && continue
        count += length(ℓ.b)
        count += length(ℓ.W)
    end
    count
end

########################################
# Word Embedding Layer functions
########################################

function string_id(word::String, vocab::Vector{String})
    for i = 1:length(vocab)
        vocab[i] == word && return i
    end
    "No match"
end

function display_nearest_words(embed::WordEmbeddingLayer, vocab::Vector{String},
                               word::String, k::Int)

    # Number of words in vocab
    vocab_size = length(vocab);

    # Find id of word
    id = string_id(word, vocab)

    # Compute distance to other words
    word_rep = embed.W[:, id]
    diff = embed.W .- word_rep
    dist = vec(sqrt(sum(diff.^2, 1)))

    # Sort and remove top item
    table = [dist 1:vocab_size]
    table = sortrows(table)[2:k+1, :]

    # Print list
    @printf("%-10s %s\n", "word", "distance")
    @printf("%s\n", "--------   --------")
    for w = 1:k
        @printf("%-10s %.2f\n", vocab[Int(table[w, 2])], table[w, 1])
    end
end

function predict_next_word(net::NeuralNet, vocab::Vector{String}, input::Tuple,
                           k::Int)

    # Find indices of the inputs words and place them in input layer
    num_inputs = length(input)
    ids = zeros(Int, num_inputs)
    for i in 1:num_inputs
        ids[i] = string_id(input[i], vocab)
    end
    net.layers[1].A[:, 1] = ids

    # Forward prop
    fwd_prop!(net)
    net.layers[end].A[:, 1]

    # Predicted next word
    vocab_size = length(vocab)
    output = [net.layers[end].A[:, 1] 1:vocab_size]
    output = sortrows(output, rev=true)[1:k, :]

    # Print completed sentences
    space = length(prod(input)) + length(input) + 9
    println("string" * " "^space * "probability")
    println("-"^space * " "^6 * "-"^11)
    for i in 1:k
        for w in input
            print("$w ")
            #@printf("%s ", w)
        end
        @printf("%-14s %.3f\n", vocab[Int(output[i, 2])], output[i, 1])
    end
end

########################################
# Demo Data Paths
########################################

# Root and dataset paths
rootpath = dirname(@__FILE__)[1:end-4]

# ngram demo data
ngrampath = joinpath(rootpath, "demo", "ngrams", "ngramdata.jld")
load_ngrams() = load(ngrampath, "train_data", "valid_data", "test_data", "vocab")

# stl10 demo data
stl10path = joinpath(rootpath, "demo", "stl10")

function load_patches()
    patches1 = load(joinpath(stl10path, "stl_sampled_patches1.jld"), "patches1")
    patches2 = load(joinpath(stl10path, "stl_sampled_patches2.jld"), "patches2")
    patches3 = load(joinpath(stl10path, "stl_sampled_patches3.jld"), "patches3")
    patches4 = load(joinpath(stl10path, "stl_sampled_patches4.jld"), "patches4")
    patches5 = load(joinpath(stl10path, "stl_sampled_patches5.jld"), "patches5")
    patches6 = load(joinpath(stl10path, "stl_sampled_patches6.jld"), "patches6")
    patches7 = load(joinpath(stl10path, "stl_sampled_patches7.jld"), "patches7")
    patches8 = load(joinpath(stl10path, "stl_sampled_patches8.jld"), "patches8")
    return hcat(patches1, patches2, patches3, patches4, patches5, patches6,
                patches7, patches8)
end

function load_train_subset()
    train_images1, train_labels1 = load(joinpath(stl10path, "stl_train_subset1.jld"), "train_images1", "train_labels1")
    train_images2, train_labels2 = load(joinpath(stl10path, "stl_train_subset2.jld"), "train_images2", "train_labels2")
    train_images3, train_labels3 = load(joinpath(stl10path, "stl_train_subset3.jld"), "train_images3", "train_labels3")
    train_images4, train_labels4 = load(joinpath(stl10path, "stl_train_subset4.jld"), "train_images4", "train_labels4")
    train_images5, train_labels5 = load(joinpath(stl10path, "stl_train_subset5.jld"), "train_images5", "train_labels5")
    train_images6, train_labels6 = load(joinpath(stl10path, "stl_train_subset6.jld"), "train_images6", "train_labels6")
    train_images7, train_labels7 = load(joinpath(stl10path, "stl_train_subset7.jld"), "train_images7", "train_labels7")
    train_images8, train_labels8 = load(joinpath(stl10path, "stl_train_subset8.jld"), "train_images8", "train_labels8")
    train_images = cat(4, train_images1, train_images2, train_images3, train_images4,
                          train_images5, train_images6, train_images7, train_images8)
    train_labels = vcat(train_labels1, train_labels2, train_labels3, train_labels4,
                        train_labels5, train_labels6, train_labels7, train_labels8)
    return train_images, train_labels
end

function load_test_subset()
    test_images1, test_labels1 = load(joinpath(stl10path, "stl_test_subset1.jld"), "test_images1", "test_labels1")
    test_images2, test_labels2 = load(joinpath(stl10path, "stl_test_subset2.jld"), "test_images2", "test_labels2")
    test_images3, test_labels3 = load(joinpath(stl10path, "stl_test_subset3.jld"), "test_images3", "test_labels3")
    test_images4, test_labels4 = load(joinpath(stl10path, "stl_test_subset4.jld"), "test_images4", "test_labels4")
    test_images5, test_labels5 = load(joinpath(stl10path, "stl_test_subset5.jld"), "test_images5", "test_labels5")
    test_images6, test_labels6 = load(joinpath(stl10path, "stl_test_subset6.jld"), "test_images6", "test_labels6")
    test_images7, test_labels7 = load(joinpath(stl10path, "stl_test_subset7.jld"), "test_images7", "test_labels7")
    test_images8, test_labels8 = load(joinpath(stl10path, "stl_test_subset8.jld"), "test_images8", "test_labels8")
    test_images9, test_labels9 = load(joinpath(stl10path, "stl_test_subset9.jld"), "test_images9", "test_labels9")
    test_images10, test_labels10 = load(joinpath(stl10path, "stl_test_subset10.jld"), "test_images10", "test_labels10")
    test_images11, test_labels11 = load(joinpath(stl10path, "stl_test_subset11.jld"), "test_images11", "test_labels11")
    test_images12, test_labels12 = load(joinpath(stl10path, "stl_test_subset12.jld"), "test_images12", "test_labels12")
    test_images13, test_labels13 = load(joinpath(stl10path, "stl_test_subset13.jld"), "test_images13", "test_labels13")
    test_images = cat(4, test_images1, test_images2, test_images3, test_images4,
                         test_images5, test_images6, test_images7, test_images8,
                         test_images9, test_images10, test_images11, test_images12,
                         test_images13)
    test_labels = vcat(test_labels1, test_labels2, test_labels3, test_labels4,
                       test_labels5, test_labels6, test_labels7, test_labels8,
                       test_labels9, test_labels10, test_labels11, test_labels12,
                       test_labels13)
    return test_images, test_labels
end

load_features() = load(joinpath(stl10path, "stl_features.jld"), "W", "b")


text:
# Intercepts for `mapreduce`, `mapfoldl` and `mapfoldr` under `op` `+`.
const plustype = Union{typeof(+), typeof(Base.add_sum)}
const type_tuple = :(Tuple{Any, $plustype, ∇ArrayOrScalar})
for f in (:mapfoldl, :mapfoldr)
    @eval begin
        import Base: $f
        @explicit_intercepts $f $type_tuple [false, false, true] #(init=0,)
        ∇(::typeof($f), ::Type{Arg{3}}, p, y, ȳ, f, ::$plustype, A::∇ArrayOrScalar) =
            hasmethod(∇, Tuple{typeof(f), Type{Arg{1}}, Real}) ?
                broadcast(An->ȳ * ∇(f, Arg{1}, An), A) :
                broadcast(An->ȳ * fmad(f, (An,), Val{1}), A)
    end
end


text:
abstract type AxisObject end

mutable struct Title <: AxisObject
    text::String
    #
    textstyle::TextStyle
    dist     ::Option{F64}  # distance labels - title
end
Title(t) = default(Title, t)

mutable struct TicksLabels <: AxisObject
    names::Vector{String}
    #
    textstyle::TextStyle
    angle    ::Option{F64}     # rotation of labels
    format   ::Option{String}  # format of the ticks labels
    shift    ::Option{F64}     # move labels to left/right
    dist     ::Option{F64}     # ⟂ distance to spine
    off      ::Bool            # whether to suppress the labels
end
TicksLabels(n) = default(TicksLabels, n)

mutable struct Ticks <: AxisObject
    places::Vector{F64}     # where the ticks are
    labels::TicksLabels     # their label
    #
    linestyle::LineStyle    # how the ticks marks look
    length   ::Option{F64}  # how long the ticks spine (negative for outside)
    symticks ::Bool         # draws ticks on 2 sides of
    off      ::Bool         # whether to suppress them
    grid     ::Bool         # ticks increased to mirrorred axis
end
Ticks(p=Float64[], l=fl2str.(p)) = default(Ticks, p, TicksLabels(l))

abstract type AxesObject end

mutable struct Legend <: AxesObject
    handles::Vector{DrawingHandle}
    labels::Vector{Union{String,Vector{String}}}
    #
    position::Option{String}
    textstyle::TextStyle
    offset  ::T2F
    bgcolor ::Option{Colorant}
    margins ::Option{T2F}
    nobox   ::Bool
    off     ::Bool
end
Legend(h, l) = default(Legend, h, l)


text:
#=
word -> length of word
* easy to generate data
* (should be) easy to train
* (should be) easy to test
=#

using Flux

data = [[1, 2, 3], [2, 3, 4], [1, 2, 3]]
labels = [3.0, 3.0, 2.0]

simple_net = Chain(
  Dense(3, 10, σ),
  Dense(10, 1))

loss(x, y) = Flux.mse(simple_net(x), y)

ps = Flux.params(simple_net)

# just use gradient descent
opt = Flux.Descent(0.01)

Flux.train!(loss, ps, zip(data, labels), opt)

# save the weights at the end
weights = Tracker.data.(params(simple_net));
using BSON: @save
@save "mymodel.bson" weights

# loading the model weights
# @load "mymodel.bson" weights
# Flux.loadparams!(simple_net, weights)


text:
function solve_one()
    data = Nothing()
    tot = 0

    for line in readlines("data/day_six.txt")
        if length(line) == 0
            tot += length(data)
            data = Nothing()
        else
            if isnothing(data)
                data = Set(line)
            else
                data = intersect(data, line)
            end
        end
    end
    return tot + length(data)
end

println(solve_one())

max_stars_repo_name: ColCarroll/adventofcode

max_issues_repo_name: ColCarroll/adventofcode

max_forks_repo_name: ColCarroll/adventofcode


text:
module Normalization
export
norm_libsize

# using StatsBase

function norm_libsize(X)
    lbsz=sum(X,dims=1)
    X=(X./lbsz)*1e4;
    return X
end

end


text:
####################################################################################
#               Constains activation functions and its derivatives
####################################################################################
export CosineSquasher
export HardTanh
export @IndicatorFunction
export @LReLU
export RampFunction
export ReLU
export Sigmoid
export @ThresholdFunction
# derivatives
export derivativeRampFunction
export derivativeReLU
"""
Threshold(polynomial, t)
Return a Threshold Function defined by `polynomial`and `t`.
Output in {-1,1}

# Arguments
- `polynomial`: A functions from real to real numbers
- `t`: A real number
"""
macro ThresholdFunction(polynomial, t::Real)
    return :( f(x)=sign($(esc(polynomial))(x) - $(esc(t))) )
end

"""
    CosineSquasher(x::Real) :: Real
Return activation function Cosine Squasher
"""
function CosineSquasher(x::Real) :: Real
    salida = 1
    if x < - pi/2
        salida =  0
    elseif x <=  pi/2
        salida =  1/2 * (1 + cos( x - pi/2))
    end
    return salida
end

"""
    IndicatorFunction(t::Real)
Return an indicator function giving its threshold `t`.
"""
macro IndicatorFunction(t::Real)
    return :( f(x)=(($(esc(t)) < x) ? 1 : 0 ) )
end

"""
    RampFunction(x::Real)
Return Ramp  function a bounded ReLU function
"""
function RampFunction(x::Real)
    return min(1,max(0,x))
end
"""
    derivateveRampFunction(x::Real)
Return the derivate of the Ramp function
"""
function derivativeRampFunction(x::Real)
    return (0<=x<=1) ? 1 : 0
end

"""
    ReLU(x::Real)
ReLU function
"""
function ReLU(x::Real)
    return max(0,x)
end
"""
    derivativeReLU(x::Real)
ReLU function
"""
function derivativeReLU(x::Real)
    return (x<0) ? 0 : 1
end

"""
    Sigmoid(x::Real)
Sigmoid activation function
"""
function Sigmoid(x::Real)
    return 1/(1+ exp(-x))
end

"""
    HardTanh(x::Real)
Hard Hyperbolic Function
"""
function HardTanh(x::Real)
    if x <= -1
        return -1
    elseif x <= 1
        return x
    else
        return 1
    end
end

"""
    LReLU(a::Real)
Leaky ReLU
"""
macro LReLU(a::Real)
    return :( f(x::Real)=(x<0) ? $(esc(a))*x : x )
end


text:
"""
    plan_route!(agent, dest, pathfinder::AStar{D})
Calculate and store the shortest path to move the agent from its current position to
`dest` (a position e.g. `(1, 5)` or `(1.3, 5.2)`) using the provided `pathfinder`.

Use this method in conjuction with [`move_along_route!`](@ref).
"""
function Agents.plan_route!(
    agent::A,
    dest::Dims{D},
    pathfinder::AStar{D},
) where {D,A<:AbstractAgent}
    path = find_path(pathfinder, agent.pos, dest)
    isnothing(path) && return
    pathfinder.agent_paths[agent.id] = path
end

"""
    plan_best_route!(agent, dests, pathfinder::AStar{D}; kwargs...)
Calculate, store, and return the best path to move the agent from its current position to
a chosen destination taken from `dests` using `pathfinder`.

The `condition = :shortest` keyword retuns the shortest path which is shortest out of the
possible destinations. Alternatively, the `:longest` path may also be requested.

Return the position of the chosen destination. Return `nothing` if none of the supplied
destinations are reachable.
"""
function Agents.plan_best_route!(
    agent::A,
    dests,
    pathfinder::AStar{D,P,M,Int64};
    condition::Symbol = :shortest,
) where {A<:AbstractAgent,D,P,M}
    @assert condition ∈ (:shortest, :longest)
    compare = condition == :shortest ? (a, b) -> a < b : (a, b) -> a > b
    best_path = Path{D,Int64}()
    best_target = nothing
    for target in dests
        path = find_path(pathfinder, agent.pos, target)
        isnothing(path) && continue
        if isempty(best_path) || compare(length(path), length(best_path))
            best_path = path
            best_target = target
        end
    end

    isnothing(best_target) && return
    pathfinder.agent_paths[agent.id] = best_path
    return best_target
end

"""
    move_along_route!(agent, model::ABM{<:GridSpace{D}}, pathfinder::AStar{D})
Move `agent` for one step along the route toward its target set by [`plan_route!`](@ref)

For pathfinding in models with [`GridSpace`](@ref).

If the agent does not have a precalculated path or the path is empty, it remains stationary.
"""
function Agents.move_along_route!(
    agent::A,
    model::ABM{<:GridSpace{D},A},
    pathfinder::AStar{D}
) where {D,A<:AbstractAgent}
    isempty(agent.id, pathfinder) && return

    move_agent!(agent, first(pathfinder.agent_paths[agent.id]), model)
    popfirst!(pathfinder.agent_paths[agent.id])
end

"""
    Pathfinding.nearby_walkable(position, model::ABM{<:GridSpace{D}}, pathfinder::AStar{D}, r = 1)
Return an iterator over all [`nearby_positions`](@ref) within "radius" `r` of the given
`position` (excluding `position`), which are walkable as specified by the given `pathfinder`.
"""
nearby_walkable(position, model::ABM{<:GridSpace{D}}, pathfinder::AStar{D}, r = 1) where {D} =
    Iterators.filter(x -> pathfinder.walkmap[x...] == 1, nearby_positions(position, model, r))

"""
    Pathfinding.random_walkable(model, pathfinder::AStar{D})
Return a random position in the given `model` that is walkable as specified by the given
`pathfinder`.
"""
random_walkable(model::ABM{<:GridSpace{D}}, pathfinder::AStar{D}) where {D} =
    Tuple(rand(model.rng, filter(x -> pathfinder.walkmap[x], CartesianIndices(model.space.s))))


text:
# compute self energy derivative
function compute_dΣ!(
    Λ     :: Float64,
    r     :: Reduced_lattice,
    m     :: Mesh,
    a     :: Action,
    da    :: Action,
    Σ_tol :: NTuple{2, Float64}
    )     :: Nothing

    # compute self energy derivative for all frequencies
    @sync for i in 2 : length(m.σ)
        Threads.@spawn begin
            integrand = v -> compute_dΣ_kernel(Λ, m.σ[i], v, r, m, a)
            da.Σ[i]   = quadgk(integrand, -Inf, -2.0 * Λ, 0.0, 2.0 * Λ, Inf, atol = Σ_tol[1], rtol = Σ_tol[2])[1]
        end
    end

    return nothing
end

# compute corrections to self energy derivative
function compute_dΣ_corr!(
    Λ     :: Float64,
    r     :: Reduced_lattice,
    m     :: Mesh,
    a     :: Action,
    da    :: Action,
    da_Σ  :: Action,
    Σ_tol :: NTuple{2, Float64}
    )     :: Nothing

    # compute first correction
    @sync for i in 2 : length(m.σ)
        Threads.@spawn begin
            integrand = v -> compute_dΣ_kernel_corr1(Λ, m.σ[i], v, r, m, a, da_Σ)
            da_Σ.Σ[i] = quadgk(integrand, -Inf, -2.0 * Λ, 0.0, 2.0 * Λ, Inf, atol = Σ_tol[1], rtol = Σ_tol[2])[1]
        end
    end

    # compute second correction and parse to da
    @sync for i in 2 : length(m.σ)
        Threads.@spawn begin
            integrand  = v -> compute_dΣ_kernel_corr2(Λ, m.σ[i], v, r, m, a, da_Σ)
            da.Σ[i]   += da_Σ.Σ[i]
            da.Σ[i]   += quadgk(integrand, -Inf, -2.0 * Λ, 0.0, 2.0 * Λ, Inf, atol = Σ_tol[1], rtol = Σ_tol[2])[1]
        end
    end

    return nothing
end


text:
# ---
# title: 1492. The kth Factor of n
# id: problem1492
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/the-kth-factor-of-n/description/>
# hidden: true
# ---
#
# Given two positive integers `n` and `k`.
#
# A factor of an integer `n` is defined as an integer `i` where `n % i == 0`.
#
# Consider a list of all factors of `n` sorted in **ascending order** , return
# _the_`kth` _factor_ in this list or return **-1** if `n` has less than `k`
# factors.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 12, k = 3
#     Output: 3
#     Explanation: Factors list is [1, 2, 3, 4, 6, 12], the 3rd factor is 3.
#
#
# **Example 2:**
#
#
#
#     Input: n = 7, k = 2
#     Output: 7
#     Explanation: Factors list is [1, 7], the 2nd factor is 7.
#
#
# **Example 3:**
#
#
#
#     Input: n = 4, k = 4
#     Output: -1
#     Explanation: Factors list is [1, 2, 4], there is only 3 factors. We should return -1.
#
#
# **Example 4:**
#
#
#
#     Input: n = 1, k = 1
#     Output: 1
#     Explanation: Factors list is [1], the 1st factor is 1.
#
#
# **Example 5:**
#
#
#
#     Input: n = 1000, k = 3
#     Output: 4
#     Explanation: Factors list is [1, 2, 4, 5, 8, 10, 20, 25, 40, 50, 100, 125, 200, 250, 500, 1000].
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= n <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=
Created on 03/02/2021 09:20:46
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Example on uncertainty propagation.
=#

struct Measurement{T<:AbstractFloat} <: Number
    x::T
    σ::T
    function Measurement(x::T, σ::T) where {T<:AbstractFloat}
        if σ < zero(σ)
            error("Measurement error should be non-zero")
        end
        new{T}(x, σ)
    end
end

val(m::Measurement) = m.x;
err(m::Measurement) = m.σ;

Base.show(io::IO, measurement::Measurement) = print(io, "$(measurement.x) ± $(measurement.σ)");

±(x, σ) = Measurement(x, σ);

# scalar multiplication
Base.:*(a::Real, m::Measurement) = a * m.x ± abs(a) * m.σ;
Base.:/(m::Measurement, a::Real) = inv(a) * m;
# adding and substracting measurements
Base.:+(m1::Measurement, m2::Measurement) = (m1.x + m2.x) ± √(m1.σ^2 + m2.σ^2);
Base.:-(m1::Measurement, m2::Measurement) = (m1.x - m2.x) ± √(m1.σ^2 + m2.σ^2);
Base.:-(m::Measurement) = -m.x ± m.σ;
# adding a constant
Base.:+(m::Measurement, a::Real) = m + (a ± zero(a));
Base.:+(a::Real, m::Measurement) = m + a;
# multiplying two measurments
Base.:*(m1::Measurement, m2::Measurement) = m1.x * m2.x ± (m1.x * m2.x) * √((m1.σ / m1.x)^2 + (m2.σ / m2.x)^2);

m1 = 9.1 ± 0.3
m2 = 8.0 ± 0.2
m3 = 10.2 ± 0.3


text:
## Kernel machines figures
using Plots
pyplot()
theme(:dark)

##
x1 = [1.0,5]
x2 = [2,-1.0]
y1 = 1.0
y2 = -1.0

v = [x1'; x2'] \ [y1,y2]

plot(framestyle=:origin,aspect_ratio=:equal,xlim=(-3,3),ylim=(-1.5,5.5))
scatter!([x1[1]],[x1[2]],label="")
annotate!([x1[1]],[x1[2]],text("+1 ",:white,:right,12))
scatter!([x2[1]],[x2[2]],label="")
annotate!([x2[1]],[x2[2]],text("-1 ",:white,:right,12))

##
plot!(size=(200,200))
savefig("kernel-interpolation-2d-setup.pdf")
##
scatter!([v[1]],[v[2]],label="")
plot!(v[1]*[-5,5], v[2]*[-5,5],label="")
annotate!([v[1]],[v[2]], text("v ",:white,:right,12))
savefig("kernel-interpolation-2d-soluiton.pdf")
##
x1 = [1.0,5,1]
x2 = [2,-1.0,-1]
y1 = 1.0
y2 = -1.0

v = [x1'; x2'] \ [y1,y2]

plot(aspect_ratio=:equal,xlim=(-3,3),ylim=(-1.5,5.5))
scatter!([x1[1]],[x1[2]],[x1[3]],label="")
annotate!([x1[1]],[x1[2]],[x1[3]],"+1 ")
scatter!([x2[1]],[x2[2]],[x2[3]],label="")
##
plot!(size=(200,200),dpi=300)
anim = @animate for i=1:180
  plot!(camera=(i,30))
end
gif(anim, "kernel-interpolation-3d-setup.gif")

##
#annotate!([x2[1]],[x2[2]],text("-1 ",:white,:right,12))
scatter!([v[1]],[v[2]],[v[3]], label="")
plot!(v[1]*[-5,5], v[2]*[-5,5], v[3]*[-5,5], label="",color=5)

## Other solutions
using LinearAlgebra
v = nullspace([x1'; x2'])
x1'*a
x1'*(a+v[:,1])

b = a+v[:,1]
scatter!([b[1]],[b[2]],[b[3]], label="")
plot!(b[1]*[-5,5], b[2]*[-5,5], b[3]*[-5,5], label="", color=8)

##
plot!(size=(200,200),dpi=300)
anim = @animate for i=1:180
  plot!(camera=(i,30))
end
gif(anim, "kernel-interpolation-3d-solution.gif")