text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ 000c1c07-fd7a-4e08-9e36-3f0aa6069309
begin
	using PlutoUI

	using CSV, DataFrames, Pipe, TabularDisplay
	using Plots, StatsPlots
end

# ╔═╡ 578530b0-9b50-11eb-07e1-5db7ceb1b1c7
md"""
## Data wrangling in Julia

Based on presentations by Tom Kwong:
  - [Data Wrangling Techniques in Julia - Part 1](https://www.youtube.com/watch?v=txme9o0EdLk)
  - [Data Wrangling Techniques in Julia - Part 2](https://www.youtube.com/watch?v=NbqQZq42gLc)

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ 6bcfba84-98ab-4d26-b98c-d86cb65dc277
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ f88a1737-7a84-4103-82c8-1920bae6b2e8
md"""
### What is data wrangling?

Take raw data (and/or unorganized data) and turn them (data as plural of datum) into something useful.

A possible data wrangling is as follows:

 read | $\Rightarrow$ | clean |  $\Rightarrow$ | Tidy
  --- | --- | --- | --- | ---
      |     |         $\Uparrow$ |           | $\Downarrow$
      |              | **Visualize** | $\Leftarrow$ | **Analyze**

**Julia ecosystem**:

Read | Clean, Tidy, Analyze | Visualize
--- | --- | --- |
CSV | DataFrames| Plots, StatsPlots
Parquet | DataFramesMeta | Makie, StatsMakie
SASLib | Query| Gadfly
ReadStat | JuliaDB, JuliaDBMeta | Gaston
JSON, JSON2/3 | | VegaLite
LightXML, EzXML | | UnicodePlots
XLSX, ExcelReader, Taro | |
PDFIO | |
"""

# ╔═╡ 02ac758d-3706-4bf4-8a36-83aff4422d0e
md"""
### Read

Utilize automated tools to read data into memory as cleanly as possible.

`CSV.jl`:
  - Skip leading/trailing rows
  - Parse missing values e.g "NA"
  - Parse date/time values
  - Normalize columns
  - Select/drop columns on read
  - User specified column types
  - Auto-delimiter detection
  - Automatic pooled string columns
  - Parse booleans
  - Transpose data
"""

# ╔═╡ fb1f4efc-6b89-48a6-95d5-8e3661edf53d
md"""
### Clean

  - Delete unwanted columns
  - Delete rows with missing data
  - Rename columns
  - Fix column data type
  - ...

Examples
```julia
describe(df)          # take a quick peek

select!(df, 1:5)      # keep "good" columns / delete unwanted ones
select!(df, Not(["garbaggecol1", "garbaggecol2"]))

rename!(df, "Region Name" => :region_name)

dropmissing!(df)      # Remove rows with missing data
dropmissing!(df, Between(:x3, :x5))
```

"""

# ╔═╡ f22c3d88-e0e6-4d1d-92d4-c99f80aa4733
md"""
### Tidy Data

  - Each variable forms a *column*
  - Each observation forms a *row*
  - Each type of observational unit forms a *table*

**Tidy Recipes, examples:**

```julia
# Stack all date columns into rows
# The column names are stored ion a new column "Date"
sdf = stack(df, Not(:Region_name); variable_name=:Date)
sdf = stack(df, Not(1); variable_name=:Date)
sdf = stack(df, 2:5; variable_name=:Date)

# Date column is a CategoricalArray type. Make it Date type
sdf.Date = [Date(get(x)) for x ∈ sdf.Date]
# or maybe using broadcast:
sdf.Date .= get.(sdf.Data) .|> Date

# Hoe to turn it back to a wide format?
df2 = unstack(sdf, :date, :value)

# It may introduce Union{Missing, T} column type. Let's fix it.
disallowmissing!(df2, 2:5)
```
"""

# ╔═╡ 2dc07426-ad94-46a0-8da2-68d687ad87c2
md"""
### Analyze

Get some insights from data.

**Analuyze Recipe**:

```julia
selct(df, :region_name, 4:5)                   # select by column

filter(:region_name => ==("Abilene, TX"), df)  # note df is the 2nd arg here
filter("2020-01-31" => >(400_000), df)

sort(df, :region_name)
sort(df, :region_name, rev=true)

# Transform means adding new columns - many more ways of using Transform
transform(df, :region_name => ByRow(length) => :region_name_len)

# Group
groupby(sdf, :Date)

# Summarize the grouped data:
combine(groupby(sdf, :Date), :value => mean => :avg)

# Using Pipe.jl to build a transformation pipeline
@pipe sdf |>
    groupby(_, :Date) |>
	combine(_, :value => mean => :avg)

country = DataFrame(Name=["United States"])

# Joining data is easy
innerjoin(df, country, on=[:region_name => :name])
leftjoin(df, country, on=[:region_name => :name])
rightjoin(df, country, on=[:region_name => :name])
outerjoin(df, country, on=[:region_name => :name])

# Also
semijoin(df, country, on=[:region_name => :name])
antijoin(df, country, on=[:region_name => :name])
```
"""

# ╔═╡ cb11e3e6-f385-44db-b739-cf056a7f8841
md"""
### Visualize

Gain more insight by looking at the data in a graphical form.

**Visualize Recipes:**

`Plots.jl` | `StatsPlots.jl`
--- | ---
`plot` | `groupedbar`
`scatter` | `corrplot`
`histogram` | `marginhist`
`heatmap` | `boxplot`
`bar` | `violin`
... | ...

"""

# ╔═╡ 944f1843-1802-45b7-a505-8bb32ad09902
md"""
## Let's practice

### Tidying
"""

# ╔═╡ f9d0761d-a859-404a-812f-f2a424c7b191
begin
	df = CSV.File("data/youth_suicide.csv", header=true) |> DataFrame;
	size(df)
end

# ╔═╡ 6790c878-d494-4247-ab1a-d5122d3d845b
first(df, 5)

# ╔═╡ ff4cc1c3-5967-4614-82eb-cc8771e962dc
names(df)

# ╔═╡ 5315e6b5-344e-42e7-a9a7-b53cfaf2b5b5
with_terminal() do
	displaytable(names(df), index=true)
end

# ╔═╡ 8083c3d3-51f6-437e-ba12-38f8b5814e1e
describe(df, :eltype, :nmissing, :first => first)

# ╔═╡ e70b9da8-34f9-4193-a86e-5e7dc9ed23aa
## turn colums 2 to 19 to row (observations)
sdf = stack(df, 2:19; variable_name=:type_year);

# ╔═╡ e783b82c-052e-49d3-87c8-70142b990d5f
size(sdf)

# ╔═╡ e9d507d7-df67-4bd4-b80b-02385b2f9be5
begin
	## Now spread content of column `:type_year` into 2 separate columns
	sdf[!, :ctype] .= split.(sdf.type_year, " ") .|> a -> getindex(a, 1)
	sdf[!, :year] .= split.(sdf.type_year, " ") .|> a -> getindex(a, 2)
	select!(sdf, Not(:type_year))
	first(sdf, 3)
end

# ╔═╡ 723dc0d9-96c4-43c0-81fd-611209e10106
## Let's check if the Total type really contains the total
## 1 - Pick a county
sdf[sdf.County .== "Yakima", :]

# ╔═╡ 0031eb09-8e5f-47fe-b04d-1600fad8f994
## Remove the total row which is redundant

sdf₁ = @pipe sdf |>
	filter(:ctype => ≠("Total"), _) |>
	filter(:year => ≠("(2008-2012)"), _);

# ╔═╡ f8ee2398-2b69-4133-b2d6-11a2c12b9d27
sdf₁[sdf₁.County .== "Yakima", :]

# ╔═╡ c52c8ce7-c68b-4b7b-b507-5a813d32cf93
describe(sdf₁)

# ╔═╡ 4b14a6ad-6dc5-4319-95bf-007e84ae6877
unique(sdf₁.year)
## Now we want integer instead of string and we want to get rid of the parentheses.

# ╔═╡ 54a8cfd0-5daa-40e3-8a0e-12413cc36e83
## mutate column year from String to Int
sdf₁[!, :year] .= replace.(sdf₁.year, r"[\(\)]" => "") |> a -> parse.(Int, a)

# ╔═╡ 3378787b-c9cf-4d47-a119-2a35abc3b0b4
describe(sdf₁, :eltype, :min, :max, :nunique, :nmissing)

# ╔═╡ 53aa180c-7260-46cf-a2a3-4a49e0c30844
begin
	# Male and Female variables can be put back as columns now
	df₁ = unstack(sdf₁, :ctype, :value)
	df₁[df₁.County .== "Yakima", :]
end

# ╔═╡ 56d8dcc1-8b76-4332-bbe4-f7081e178300
md"""
### Visualize
"""

# ╔═╡ 73a99bae-e8a6-44d1-af74-9778e17df7b9
# Plot King county's yearly rate
filter(:County => ==("King"), df₁)

# ╔═╡ 5eaa88b0-95e4-4410-bd99-cb9ad40fdb79
# bar chart
@pipe df₁ |>
	filter(:County => ==("King"), _) |>
	bar(_.year, _.Male,
        title = "King County Youth Suicides (Male)",
        legend = :none,
        size = (450, 300))

# ╔═╡ 5548eaac-d390-4cb1-bf49-0b8985c1f0f0
begin
	# Prepare to plot both Male and Female together
	df_stacked = stack(df₁, 3:4; variable_name = "gender");
	names(df_stacked)
end

# ╔═╡ bb947d94-62f3-45f6-a0d1-9ed596e289c8
@pipe df_stacked |>
	filter(:County => ==("King"), _) |>
	 groupedbar(    # StatsPlots.jl
		_.year,     # x-axis
	    _.value;    # y-axis
        group = _.gender,
        bar_position = :stack,
        bar_width = 0.7,
        title = "King County Suicides",
        size = (450, 300),
        legend = :topleft)

# ╔═╡ 5fadd975-cbcf-4a0e-b493-0382875354ba
@pipe df₁ |>
    groupby(_, :County)

# ╔═╡ a92c1c18-39bd-4dcf-81a7-f2f62577e156
# 5 year totals
@pipe df₁ |>
    groupby(_, :County) |>
    combine(_, :Female => sum, :Male => sum)

# ╔═╡ 2c21d894-7e57-4f92-b08d-a4d527f5c25a
# 5 year totals with male/female combined
@pipe df₁ |>
    groupby(_, :County) |>
    combine(_, :Female => sum, :Male => sum) |>
    select(_, :County, [:Female_sum, :Male_sum] => ByRow(+) => :Total)

# ╔═╡ 72aeaa39-6529-4fc3-ba86-f2b113a5a075
let
    data₁ = @pipe df₁ |>
        groupby(_, :year) |>
        combine(_, :Female => sum => :Female, :Male => sum => :Male) |>
        stack(_, 2:3; variable_name = :gender, value_name = :suicides)

    plot(data₁.year, data₁.suicides; groups = data₁.gender,
        title = "Youth Suicides Trend",
        legend = :topleft,
        linewidth = 3,
        size = (450, 300))
end

# ╔═╡ a7802df0-6819-4650-a539-ccbb99333e92
let
    data₂ = @pipe df₁ |>
        select(_, :County, :year, [:Female, :Male] => ByRow(+) => :total) |>
        unstack(_, :year, :total)

    values = Matrix(data₂[:, 2:end])

    StatsPlots.heatmap(2008:2012, data₂.County, values;
        title = "Youth Suicide Heatmap",
        xticks = :all,
        yticks = :all,
        size = (400, 600))
end

# ╔═╡ 8018394c-a123-43c4-aaaf-3d942284e182
html"""

"""

# ╔═╡ Cell order:
# ╟─578530b0-9b50-11eb-07e1-5db7ceb1b1c7
# ╠═000c1c07-fd7a-4e08-9e36-3f0aa6069309
# ╟─6bcfba84-98ab-4d26-b98c-d86cb65dc277
# ╟─f88a1737-7a84-4103-82c8-1920bae6b2e8
# ╟─02ac758d-3706-4bf4-8a36-83aff4422d0e
# ╟─fb1f4efc-6b89-48a6-95d5-8e3661edf53d
# ╟─f22c3d88-e0e6-4d1d-92d4-c99f80aa4733
# ╟─2dc07426-ad94-46a0-8da2-68d687ad87c2
# ╟─cb11e3e6-f385-44db-b739-cf056a7f8841
# ╟─944f1843-1802-45b7-a505-8bb32ad09902
# ╠═f9d0761d-a859-404a-812f-f2a424c7b191
# ╠═6790c878-d494-4247-ab1a-d5122d3d845b
# ╠═ff4cc1c3-5967-4614-82eb-cc8771e962dc
# ╠═5315e6b5-344e-42e7-a9a7-b53cfaf2b5b5
# ╠═8083c3d3-51f6-437e-ba12-38f8b5814e1e
# ╠═e70b9da8-34f9-4193-a86e-5e7dc9ed23aa
# ╠═e783b82c-052e-49d3-87c8-70142b990d5f
# ╠═e9d507d7-df67-4bd4-b80b-02385b2f9be5
# ╠═723dc0d9-96c4-43c0-81fd-611209e10106
# ╠═0031eb09-8e5f-47fe-b04d-1600fad8f994
# ╠═f8ee2398-2b69-4133-b2d6-11a2c12b9d27
# ╠═c52c8ce7-c68b-4b7b-b507-5a813d32cf93
# ╠═4b14a6ad-6dc5-4319-95bf-007e84ae6877
# ╠═54a8cfd0-5daa-40e3-8a0e-12413cc36e83
# ╠═3378787b-c9cf-4d47-a119-2a35abc3b0b4
# ╠═53aa180c-7260-46cf-a2a3-4a49e0c30844
# ╟─56d8dcc1-8b76-4332-bbe4-f7081e178300
# ╠═73a99bae-e8a6-44d1-af74-9778e17df7b9
# ╠═5eaa88b0-95e4-4410-bd99-cb9ad40fdb79
# ╠═5548eaac-d390-4cb1-bf49-0b8985c1f0f0
# ╠═bb947d94-62f3-45f6-a0d1-9ed596e289c8
# ╠═5fadd975-cbcf-4a0e-b493-0382875354ba
# ╠═a92c1c18-39bd-4dcf-81a7-f2f62577e156
# ╠═2c21d894-7e57-4f92-b08d-a4d527f5c25a
# ╠═72aeaa39-6529-4fc3-ba86-f2b113a5a075
# ╠═a7802df0-6819-4650-a539-ccbb99333e92
# ╟─8018394c-a123-43c4-aaaf-3d942284e182

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
export MutablePose2Pose2Gaussian, PackedMutablePose2Pose2Gaussian
# import IncrementalInference: getSample

"""
    $TYPEDEF

Specialized Pose2Pose2 factor type (Gaussian), which allows for rapid accumulation of odometry information as a branch on the factor graph.
"""
mutable struct MutablePose2Pose2Gaussian  <: IIF.FunctorPairwise
  Zij::MvNormal
end
function getSample(fct::MutablePose2Pose2Gaussian, N::Int=100)
  return (rand(fct.Zij, N), )
end

"""
    $SIGNATURES

Residual function for MutablePose2Pose2Gaussian.

Related

Pose2Pose2, Pose3Pose3, InertialPose3, DynPose2Pose2, Point2Point2, VelPoint2VelPoint2
"""
function (::MutablePose2Pose2Gaussian)(
                 res::Vector{Float64},
                 userdata,
                 idx::Int,
                 meas::Tuple,
                 wxi::Array{Float64,2},
                 wxj::Array{Float64,2}  )::Nothing
  #
  wXjhat = SE2(wxi[1:3,idx])*SE2(meas[1][1:3,idx])
  jXjhat = SE2(wxj[1:3,idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end

## Serialization support

"""
$(TYPEDEF)
"""
mutable struct PackedMutablePose2Pose2Gaussian  <: IIF.PackedInferenceType
  datastr::String
  PackedMutablePose2Pose2Gaussian() = new()
  PackedMutablePose2Pose2Gaussian(x::String) = new(x)
end
function convert(::Type{MutablePose2Pose2Gaussian}, d::PackedMutablePose2Pose2Gaussian)
  return MutablePose2Pose2Gaussian(extractdistribution(d.datastr))
end
function convert(::Type{PackedMutablePose2Pose2Gaussian}, d::MutablePose2Pose2Gaussian)
  return PackedMutablePose2Pose2Gaussian(string(d.Zij))
end

ext: jl


text:
function render(mdp::SimpleGridWorld, step::Union{NamedTuple,Dict};
                color = s->reward(mdp, s),
                policy::Union{Policy,Nothing} = nothing
               )

    nx, ny = mdp.size
    cells = []
    for x in 1:nx, y in 1:ny
        cell = cell_ctx((x,y), mdp.size)
        if policy !== nothing
            a = action(policy, GWPos(x,y))
            txt = compose(context(), text(0.5, 0.5, aarrow[a], hcenter, vcenter), stroke("black"))
            compose!(cell, txt)
        end
        clr = tocolor(color(GWPos(x,y)))
        compose!(cell, rectangle(), fill(clr), stroke("gray"))
        push!(cells, cell)
    end
    grid = compose(context(), linewidth(0.5mm), cells...)
    outline = compose(context(), linewidth(1mm), rectangle(), stroke("gray"))

    if haskey(step, :s)
        agent_ctx = cell_ctx(step[:s], mdp.size)
        agent = compose(agent_ctx, circle(0.5, 0.5, 0.4), fill("orange"))
    else
        agent = nothing
    end

    sz = min(w,h)
    return compose(context((w-sz)/2, (h-sz)/2, sz, sz), agent, grid, outline)
end

function cell_ctx(xy, size)
    nx, ny = size
    x, y = xy
    return context((x-1)/nx, (ny-y)/ny, 1/nx, 1/ny)
end

tocolor(x) = x
function tocolor(r::Float64)
    minr = -10.0
    maxr = 10.0
    frac = (r-minr)/(maxr-minr)
    return get(ColorSchemes.redgreensplit, frac)
end

const aarrow = Dict(:up=>'↑', :left=>'←', :down=>'↓', :right=>'→')


text:
module Planar_truss_examples
using FinEtools
using FinEtools.FENodeSetModule
using FinEtools.MeshExportModule

function Planar_truss()
    # Planar truss  structure loaded with concentrated forces  at some nodes.
    # The nodal displacements should be
    # 0         0
    # 0         0
    # 0.0213    0.0408
    # -0.0160    0.0462
    # 0.0427    0.1501
    # -0.0053    0.1661
    Area = 1.5
    E = 1.0e7 # Young's modulus
    nu = 0.0
    alpha = 0.0
    fens = FENodeSetModule.FENodeSet([0.0 0;
    0 40;
    40 0;
    40 40;
    80 0;
    80 40] )
    fes = FESetL2([1     3
    1     4
    2     4
    3     4
    3     5
    5     4
    6     4
    5     6])

    MR = DeforModelRed1D
    material = MatDeforElastIso(MR,  0.0, E, nu, alpha)
    # display(material )

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field
    setebc!(u, 1)
    setebc!(u, 2)
    applyebc!(u)
    numberdofs!(u)
    display(u)

    integdata = IntegData(fes, GaussRule(1, 1), (loc, conn, N) -> Area, false)
    femm = FEMMDeforLinear(MR, integdata, CSys(2, 1), material)
    display(femm.mcsys)
    K = stiffness(femm,  geom,  u)

    fi = ForceIntensity(vec([0 -2000.0]))
    lfemm = FEMMBase(IntegData(FESetP1(reshape([3], 1,1)), PointRule()))
    F = distribloads(lfemm,  geom,  u,  fi,  3);
    fi = ForceIntensity(vec([+2000.0 0]))
    lfemm = FEMMBase(IntegData(FESetP1(reshape([5], 1,1)), PointRule()))
    F = F + distribloads(lfemm,  geom,  u,  fi,  3);
    fi = ForceIntensity(vec([+4000.0 +6000.0]))
    lfemm = FEMMBase(IntegData(FESetP1(reshape([6], 1,1)), PointRule()))
    F = F + distribloads(lfemm,  geom,  u,  fi,  3);

    K = cholesky(K)
    U=  K\F
    scattersysvec!(u, U[:])

    sfld =  elemfieldfromintegpoints(femm, geom, u, :Cauchy, 1)
    display(sfld)
    println("Cauchy = $(sfld.values)")
    vfld =  elemfieldfromintegpoints(femm, geom, u, :vm, 1)
    display(vfld)

    File = "Planar_truss.vtk"
    MeshExportModule.vtkexportmesh(File, fens, fes;
    scalars=[("sx", sfld.values), ("vm", vfld.values)])
    @async run(`"paraview.exe" $File`)
    # try rm(File) catch end

end # Planar_truss

function allrun()
    println("#####################################################")
    println("# Planar_truss ")
    Planar_truss()
    return true
end # function allrun

end # module Planar_truss_examples


text:
simple = @stochastic_model begin
    @stage 1 begin
        @decision(model, x₁ >= 40)
        @decision(model, x₂ >= 20)
        objective = @expression(model, 100*x₁ + 150*x₂)
        @objective(model, Min, objective)
        @constraint(model, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @uncertain q₁ q₂ d₁ d₂
        @recourse(model, 0 <= y₁ <= d₁)
        @recourse(model, 0 <= y₂ <= d₂)
        @objective(model, Max, q₁*y₁ + q₂*y₂)
        @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end

ξ₁ = @scenario q₁ = 24.0 q₂ = 28.0 d₁ = 500.0 d₂ = 100.0 probability = 0.4
ξ₂ = @scenario q₁ = 28.0 q₂ = 32.0 d₁ = 300.0 d₂ = 300.0 probability = 0.6

simple_res = SPResult([46.67, 36.25], Dict(1 => [300., 100.], 2 => [300., 100.]), -855.83, -1518.75, 662.92, 286.92, -1445.92, -568.92)
push!(problems, (simple, [ξ₁,ξ₂], simple_res, "Simple"))


text:
using ArgParse

function getclargs()
    settings = ArgParseSettings()

    @add_arg_table! settings begin
        "--mode"
        help = "if true simulation is ran sequentially, distributed and threaded"
        arg_type = String
        default = "distributed"
        "--ncores"
        help = "number of cores"
        arg_type = Int
        default = length(Sys.cpu_info()) - 1
        "--nbits"
        help = "number of bits"
        arg_type = Int
        default = 10
        "--tbit"
        help = "bit duration"
        arg_type = Float64
        default = 50.0
        "--ntrials"
        help = "number of trials"
        arg_type = Int
        default = 16
        "--coupling-strength"
        help = "coupling strength between the nodes"
        arg_type = Float64
        default = 20.0
        "--high"
        help = "high level of control link weight"
        arg_type = Float64
        default = 100.0
        "--low"
        help = "low level of control link weight"
        arg_type = Float64
        default = eps()
        "--srcnode"
        help = "source node of control link"
        arg_type = Int
        default = 5
        "--dstnode"
        help = "destination node of control link"
        arg_type = Int
        default = 8
        "--time-scaling"
        help = "time scaling used to scale network dynamics"
        arg_type = Float64
        default = 1.0
        "--dt"
        help = "sampling period"
        arg_type = Float64
        default = 0.01
        "--minsnr"
        help = "minimum snr level"
        arg_type = Int
        default = 0
        "--maxsnr"
        help = "maximum snr level"
        arg_type = Int
        default = 12
        "--stepsnr"
        help = "number of snr level"
        arg_type = Int
        default = 2
        "--simdir"
        help = "simulation directory"
        arg_type = String
        default = "/tmp"
        "--maxiters"
        help = "maximum number of iterations of solver"
        arg_type = Int
        default = typemax(Int)
        "--simprefix"
        help = "simulation prefix"
        arg_type = String
        default = "Simulation-"
        "--savenoise"
        help = "if true noise is saved in data files"
        action = :store_true
        "--loglevel"
        help = "logging level"
        arg_type = String
        default = "Info"
    end

    return parse_args(settings)
end


text:
"""
    convective_adjust!(x)

Remove negative gradients from `x` (generally a temperature profile).
"""
function convective_adjust!(x)
    for i in length(x)-3:-1:2
        if x[i] > x[i+1]
            if x[i-1] > x[i]
                x[i] = x[i+1]
            else
                x[i] = (x[i-1] + x[i+1]) / 2
            end
        end
    end
end


text:
include("header.jl")

@testset "loss" begin
    for f in (logp, logsumexp)
        a = rand(10,10)
        @test gradcheck(f,a)
        @test gradcheck(f,a,kw=(:dims=>1,))
        @test gradcheck(f,a,kw=(:dims=>2,))
        if gpu() >= 0
            k = KnetArray(a)
            @test gradcheck(f,k)
            @test gradcheck(f,k,kw=(:dims=>1,))
            @test gradcheck(f,k,kw=(:dims=>2,))
            @test isapprox(f(a),f(k))
            @test isapprox(f(a,dims=1),f(k,dims=1))
            @test isapprox(f(a,dims=2),f(k,dims=2))
        end

        a = rand(10,10,10)
        @test gradcheck(f,a)
        @test gradcheck(f,a,kw=(:dims=>1,))
        @test gradcheck(f,a,kw=(:dims=>2,))
        @test gradcheck(f,a,kw=(:dims=>3,))
        @test gradcheck(f,a,kw=(:dims=>(1,2),))
        @test gradcheck(f,a,kw=(:dims=>(3,2),))
        @test gradcheck(f,a,kw=(:dims=>(1,3),))

        if gpu() >= 0
            k = KnetArray(a)
            @test gradcheck(f,k)
            @test isapprox(f(a),f(k))
            for d in [1,2,3]
                @test gradcheck(f,k,kw=(:dims=>d,))
                @test isapprox(f(a,dims=d),f(k,dims=d))
            end
            for dims in [(1,2), (1,3), (3,1), [1,3,2]]
                @test isapprox(f(a,dims=dims),f(k,dims=dims))
            end
        end

        a = rand(10,10, 10)
        for d in [1, 2, 3, (1,2), (1,3), [2,3], [1,2,3]]
            @test softmax(a, dims=d) ≈ exp.(logsoftmax(a, dims=d))
            @test all(sum(softmax(a, dims=d), dims=d) .≈ 1)
            if gpu() > 0
                k = KnetArray(a)
                @test softmax(k, dims=d) ≈ exp.(logsoftmax(k, dims=d))
                @test all(sum(softmax(k, dims=d), dims=d) .≈ 1)
            end
        end
    end

    a = rand(10,10)
    indices = rand(1:10,10)
    @test gradcheck(nll, a, indices, kw=(:dims=>1,), args=1)
    @test gradcheck(nll, a, indices, kw=(:dims=>2,), args=1)
    if gpu() >= 0
        k = KnetArray(a)
        @test gradcheck(nll, k, indices, kw=(:dims=>1,), args=1)
        @test gradcheck(nll, k, indices, kw=(:dims=>2,), args=1)
        @test isapprox(nll(k, indices, dims=1), nll(a, indices, dims=1))
        @test isapprox(nll(k, indices, dims=2), nll(a, indices, dims=2))
    end
    @test gradcheck(logistic,a[:],a[:])
    @test gradcheck(bce,a[:],a[:])
end


text:
# Make a given Tree monotone increasing

# First idea just bring children nodes to parent node level or above

#for i in parent
#    for j in children of i
#        diff = state[j] - state[i]
#        if diff < 0
#            state[j] = state[i]
#        end
#    end
#end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ fe6fc4e0-688c-11eb-3cba-fbadc8c45088
using Plots

# ╔═╡ b3150b20-68a7-11eb-0c1c-1389cd8227ec
html"""
CellularBase docs

"""

# ╔═╡ ec8999c0-68a7-11eb-02dd-dfb6f5672595
md"This is the documentation of [CellularBase.jl](https://github.com/DhruvaSambrani/idc621/blob/main/CellularBase.jl)"

# ╔═╡ 37e3fad0-688c-11eb-0760-991b2ada6775
md"# Boundary Conditions

Define the kinds of Boundary Conditions

- Periodic
- FixedMin - Fixed BC with all elements of the min of the `possible_states`
- FixedMax - Fixed BC with all elements of the max of the `possible_states`
- Clamp - Nearest cell within the region of interest

"

# ╔═╡ 29996140-688c-11eb-1317-6d36dc7fe265
@enum BoundaryCondition Periodic FixedMin FixedMax Clamp

# ╔═╡ a222c570-688c-11eb-0a24-bf1da8cab724
md"# Neighborhoods
Every neighborhood has a `radius` and a set of `CartesianIndex`es associated with it. The most general way of evaluating these is by referencing the appropriate property.
"

# ╔═╡ 2ddfe490-688c-11eb-11fe-7fe9cf9cc425
begin
	# Neighborhoods
	abstract type AbstractNeighborhood end
	radius(neighborhood::AbstractNeighborhood) = neighborhood.radius
	cartesians(neighborhood::AbstractNeighborhood) = neighborhood.cartesians
end

# ╔═╡ eabf0730-688c-11eb-25cd-0fdf3696c423
md"## VonNeumann Neighborhood

The VonNeumann Neighborhood is the simplest neighborhood which contains $k$ neighbors adjacent to it.

It is defined for multiple dimensions.

A 2D representation of radius(=$k$) 3 is shown below
"

# ╔═╡ a1003e70-688c-11eb-12c6-193712a896ff
struct VonNeumannNeighborhood <: AbstractNeighborhood
	cartesians::Array{CartesianIndex}
	radius::Int
	dimensions::Int
	function VonNeumannNeighborhood(radius::Int, dimensions::Int)::VonNeumannNeighborhood
		arr = CartesianIndex[CartesianIndex(Tuple(zeros(Int, dimensions)))]
		for i in 1:radius
			for j in 1:dimensions
				zs = zeros(Int, dimensions)
				zs[j] = i
				push!(arr, CartesianIndex(Tuple(zs)))
				zs[j] = -i
				push!(arr, CartesianIndex(Tuple(zs)))
			end
		end
		new(arr, radius, dimensions)
	end
end

# ╔═╡ af3a67d0-688d-11eb-1d4d-bff9ccfbd737
begin
	vline([-4, 4], color = :black)
	hline!([-4, 4], color = :black)
	scatter!(Tuple.(cartesians(VonNeumannNeighborhood(3, 2))), grid=false, legend=false)
	scatter!((0,0), size=(240, 160), axis = false, ticks = false, lims = (-4,4))
end

# ╔═╡ 4241b7e0-688e-11eb-3c30-f7c0f1730888
md"## Moore Neighborhood

The Moore Neighborhood is a simple neighborhood which contains all neighbours in the $ n$ d-cube centered at the cell and of side length $2k+1$

It is defined for multiple dimensions.

A 2D representation of radius(= $k$) 3 is shown below
"

# ╔═╡ 3752f380-688e-11eb-31bd-699522b02cbd
struct MooreNeighborhood <: AbstractNeighborhood
	cartesians::Array{CartesianIndex}
	radius::Int
	dimensions::Int
	function MooreNeighborhood(radius::Int, dimensions::Int)::MooreNeighborhood
		new(
			CartesianIndices(Tuple(-radius:radius for _ in 1:dimensions)) |>
			x -> reshape(x, :),
			radius, dimensions
		)
	end
end

# ╔═╡ d069c712-688e-11eb-38c1-47e1dec2e2aa
begin
	vline([-4, 4], color = :black)
	hline!([-4, 4], color = :black)
	scatter!(Tuple.(cartesians(MooreNeighborhood(3, 2))), grid=false, legend=false)
	scatter!((0,0), size=(240, 160), axis = false, ticks = false, lims = (-4,4))
end

# ╔═╡ 61003e30-688f-11eb-1732-3dc61a8abad2
md"# Grids
These are the actual cellular automata which are the aim of this package.

`AbstractGrid` is concretized to a working grid which has the real working mechanism
"

# ╔═╡ 43d5cff0-688f-11eb-0cbe-d32c52e694cf
abstract type AbstractGrid{T} end

# ╔═╡ a86140ce-688f-11eb-21eb-c3cd76fcdcd7
md"## `getindex`
`Base.getindex` is specialized to AbstractGrids to handle Boundary conditions
"

# ╔═╡ e27b1022-688f-11eb-0aee-85fb8c9d3834
md"""## State functions

These are functions related to the _state_ of the grid.

!!! warning "Optimization Possibility"
    `newstate` and `state!` can often be specialized and optimized to be more efficient and hence throw warnings
"""

# ╔═╡ f6fb04b0-688f-11eb-16a7-659171a260de
begin
	function state(grid::AbstractGrid{T})::Array{T} where {T}
		grid.state
	end

	function state!(grid::AbstractGrid{T}, newstate)::Nothing  where {T}
		@warn "Possible unspecialized call"
		grid.state = newstate
		return nothing
	end

	function newstate(grid::AbstractGrid{T})::Array{T} where {T}
		@warn "Possible unspecialized call"
		similar(state(grid))
	end
end;

# ╔═╡ 1f21cbe0-6890-11eb-309b-a7512538a923
md"## Accessor functions"

# ╔═╡ feac0330-688f-11eb-03a7-f7039760a343
begin
	boundary_conditions(grid::AbstractGrid)::BoundaryCondition = grid.bc
	neighborhood(grid::AbstractGrid)::Array{CartesianIndex} = cartesians(grid.neighborhood)
	function possible_states(grid::AbstractGrid{T})::Array{T} where {T}
		grid.possible_states
	end
	function neighbors(grid::AbstractGrid{T}, i::CartesianIndex)::Array{T} where {T}
		neighborindex = neighborhood(grid) .|> x -> x + i
		grid[neighborindex]
	end
end;

# ╔═╡ 5b49fe40-688f-11eb-333c-c722511ca6e0
function Base.getindex(grid::AbstractGrid{T}, index::CartesianIndex)::T where {T}
    try
        state(grid)[index]
    catch
        if boundary_conditions(grid) == Periodic
            grid[CartesianIndex(mod1.(Tuple(index), size(grid)))] # mod-1 arithmetic
        elseif boundary_conditions(grid) == FixedMax
            maximum(possible_states(grid))
        elseif boundary_conditions(grid) == FixedMin
            minimum(possible_states(grid))
        else
            grid[
                CartesianIndex(
                    Tuple(
                        clamp(Tuple(index)[i], 1, size(grid)[i])
                        for i in 1:ndims(grid)
                    )
                )
            ]
        end
    end
end

# ╔═╡ c6e11070-6890-11eb-3f5a-051919a603b2
md"## Useful Helpers"

# ╔═╡ ea9a4450-6890-11eb-334e-0d4d584222db
begin
	Base.ndims(grid::AbstractGrid) = ndims(state(grid))
	Base.size(grid::AbstractGrid) = size(state(grid))
	function Base.getindex(grid::AbstractGrid{T}, indexs::Union{Array{Int},Array{CartesianIndex{N}}})::Array{T} where {T,N}
		indexs .|> i -> grid[i]
	end

	function Base.getindex(grid::AbstractGrid{T}, is::Int...)::T where T
		grid[CartesianIndex(is)]
	end
	Base.CartesianIndices(grid::AbstractGrid)::CartesianIndices = CartesianIndices(state(grid))
end

# ╔═╡ 0c8a9f10-6891-11eb-044c-abe0b2b66629
md"# Evolutions and Simulations
This section defines the most abstract functions which evolve and simulate a Cellular Automaton
"

# ╔═╡ 4a177510-6891-11eb-12ca-35c2f8dad57d
md"""## Evolutions

`evolve!` is a mutator function that evolves the grid by one step. A new state is created by calling `newstate`. Each index is iterated through, and the values of the neighbours are sent as parameters to the `grid`. The grid then evaluates the next state of the present index, and hence the evolution rules must be coded into the `grid` as a functional object

!!! danger "Grid is a Function too!"
    Every AbstractGrid is also assumed to be a function-like object and the code will throw an error if it is not defined to be as such

The new state is then set back into the grid.

`tabular_evolution!` is a helper function that can be used for rules which can be encoded into a table, like Wolfram's CA.

!!! info "tabular_evolution!"
    This function can be used in a specialized `evolve!`. Send the table via the `kwargs` to `evolve!` and then unpack to `tabular_evolution!`. See [caWolfram.jl](https://github.com/DhruvaSambrani/idc621/blob/main/caWolfram.jl) for an implementation.

"""

# ╔═╡ 2ea1bbb0-6891-11eb-3c44-2706f3faa991
begin
	function evolve!(grid::AbstractGrid{T}; kwargs...)::Nothing where {T}
	    _newstate = newstate(grid)
	    Threads.@threads for i in CartesianIndices(grid)
	        _newstate[i] = grid(neighbors(grid, i); kwargs...)
		end
		state!(grid, _newstate)
		return nothing
	end

	function tabular_evolution!(grid::AbstractGrid, table::Dict)::Nothing
	    newstate = newstate(grid)
	    for i in CartesianIndices(grid)
	        newstate[i] = table[grid(neighbors(grid, i))]
		end
		state!(grid, newstate)
		return nothing
	end
end

# ╔═╡ bc977ad0-6892-11eb-065f-8561963d653e
md"""## Simulation
The outer most function that will almost always be called _unspecialized_ by the user.
The `simulate!` function takes a grid and a number of steps to simulate for. The `store_results` parameter allows users to not store results after every `evolve!` call. `postrunhook` is called after every step with parameters as `(step, grid, kwargs...)`. All `kwargs` are forwarded to `evolve!` and `postrunhook`.

!!! tips "postrunhook Special Return Values"
    While `postrunhook` can mutate the grid on specific steps, it can also act as a messenger to `simulate!` to change its behavior. See below for available messages

List of special messages-

1. `:shortcurcuit` ⟹ do not call `evolve!` on the grid in the next step. This can be useful if you want to stop evolution for some steps, or if you know that the grid will not change anymore, such as when it has reached a stable state.
2. `:interrupt` ⟹ return immediately.

!!! info "Implementation Example"
    See [Term Paper 1#Equilibrium infection distribution](https://dhruvasambrani.github.io/idc621/termpaper1) for example

"""

# ╔═╡ e450da2e-6892-11eb-21e4-c3e5d6141a6d
function simulate!(grid::AbstractGrid{T}, max_steps::Int; store_results=true, postrunhook=nothing, kwargs...)::Array{Array{T}}  where {T}
	if store_results results = Array{T}[] end
	push!(results, copy(state(grid)))
	if !isnothing(postrunhook)
		postrunhook(grid, step; kwargs...)
	end
	ss = false
	step = 1
	while step < max_steps
		if (!ss) evolve!(grid; kwargs...) end
		if store_results push!(results, copy(state(grid))) end
		if !isnothing(postrunhook)
			message = postrunhook(grid, step; kwargs...)
			ss = message == :shortcurcuit
			if message == :interrupt
				break
			end
		end
	end
	if store_results results else nothing end
end

# ╔═╡ Cell order:
# ╟─b3150b20-68a7-11eb-0c1c-1389cd8227ec
# ╟─ec8999c0-68a7-11eb-02dd-dfb6f5672595
# ╠═fe6fc4e0-688c-11eb-3cba-fbadc8c45088
# ╟─37e3fad0-688c-11eb-0760-991b2ada6775
# ╠═29996140-688c-11eb-1317-6d36dc7fe265
# ╟─a222c570-688c-11eb-0a24-bf1da8cab724
# ╠═2ddfe490-688c-11eb-11fe-7fe9cf9cc425
# ╟─eabf0730-688c-11eb-25cd-0fdf3696c423
# ╟─af3a67d0-688d-11eb-1d4d-bff9ccfbd737
# ╠═a1003e70-688c-11eb-12c6-193712a896ff
# ╟─4241b7e0-688e-11eb-3c30-f7c0f1730888
# ╟─d069c712-688e-11eb-38c1-47e1dec2e2aa
# ╠═3752f380-688e-11eb-31bd-699522b02cbd
# ╟─61003e30-688f-11eb-1732-3dc61a8abad2
# ╠═43d5cff0-688f-11eb-0cbe-d32c52e694cf
# ╟─a86140ce-688f-11eb-21eb-c3cd76fcdcd7
# ╠═5b49fe40-688f-11eb-333c-c722511ca6e0
# ╟─e27b1022-688f-11eb-0aee-85fb8c9d3834
# ╠═f6fb04b0-688f-11eb-16a7-659171a260de
# ╟─1f21cbe0-6890-11eb-309b-a7512538a923
# ╠═feac0330-688f-11eb-03a7-f7039760a343
# ╟─c6e11070-6890-11eb-3f5a-051919a603b2
# ╠═ea9a4450-6890-11eb-334e-0d4d584222db
# ╟─0c8a9f10-6891-11eb-044c-abe0b2b66629
# ╟─4a177510-6891-11eb-12ca-35c2f8dad57d
# ╠═2ea1bbb0-6891-11eb-3c44-2706f3faa991
# ╟─bc977ad0-6892-11eb-065f-8561963d653e
# ╠═e450da2e-6892-11eb-21e4-c3e5d6141a6d


text:
function count_move!(r, side :: HorizonSide, count)
    for i = 1:count
        move!(r,side)
    end
end

function move_count_mark!(r,side::HorizonSide, level, len)
    for i in 1:(len - level)
        putmarker!(r)
        move!(r, side)
    end
    putmarker!(r)
end

function movements_to_side!(r,side::HorizonSide)
    count = 0
    while isborder(r,side) ==  false
        move!(r,side)
        count = count + 1
    end
    return count
end

function move_mark_to_side!(r,side::HorizonSide)

    count = 0

    while isborder(r,side) ==  false
        putmarker!(r)
        move!(r,side)
        count = count + 1
    end
    putmarker!(r)
    return count
end

function corner(r)
    x = 0
    y = 0
    while (isborder(r, Sud) == false) || (isborder(r, West) == false)
        y = y + movements_to_side!(r, Sud)
        x = x + movements_to_side!(r, West)
    end
    return x, y
end

function go_home!(r, x, y)
    corner(r)
    count_move!(r, Nord, y)
    count_move!(r, Ost, x)
end

function ladder(r)

    level = 1
    len = 0
    x = 0
    y = 0

    x, y = corner(r)

    len = move_mark_to_side!(r, Ost)
    movements_to_side!(r, West)

    while (isborder(r, Nord) == false)
        move!(r, Nord)
        move_count_mark!(r, Ost, level, len)
        movements_to_side!(r, West)
        level = level + 1
    end

    go_home!(r, x, y)

end

ladder(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
export iauTpstv
"""
In the tangent plane projection, given the star's rectangular
coordinates and the direction cosines of the tangent point, solve
for the direction cosines of the star.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   xi,eta  double     rectangular coordinates of star image (Note 2)
   v0      double[3]  tangent point's direction cosines

Returned:
   v       double[3]  star's direction cosines

   1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

   2. The eta axis points due north in the adopted coordinate system.
   If the direction cosines represent observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  If the direction cosines are with
   respect to a right-handed triad, (xi,eta) are also right-handed.
   The units of (xi,eta) are, effectively, radians at the tangent
   point.

   3. The method used is to complete the star vector in the (xi,eta)
   based triad and normalize it, then rotate the triad to put the
   tangent point at the pole with the x-axis aligned to zero
   longitude.  Writing (a0,b0) for the celestial spherical
   coordinates of the tangent point, the sequence of rotations is
   (b-pi/2) around the x-axis followed by (-a-pi/2) around the
   z-axis.

   4. If vector v0 is not of unit length, the returned vector v will
   be wrong.

   5. If vector v0 points at a pole, the returned vector v will be
   based on the arbitrary assumption that the longitude coordinate
   of the tangent point is zero.

   6. This function is a member of the following set:

      spherical      vector         solve for

      iauTpxes      iauTpxev         xi,eta
      iauTpsts    > iauTpstv <        star
      iauTpors      iauTporv         origin

References:

   Calabretta M.R. & Greisen, E.W., 2002, "Representations of
   celestial coordinates in FITS", Astron.Astrophys. 395, 1077

   Green, R.M., "Spherical Astronomy", Cambridge University Press,
   1987, Chapter 13.

This revision:   2018 January 2

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauTpstv(double xi, double eta, double v0[3], double v[3])

function iauTpstv(xi::Real, eta::Real, v0::AbstractVector{<:Real})
   v = zeros(Float64, 3)

   ccall((:iauTpstv, libsofa_c), Cvoid,
       (Cdouble, Cdouble,
       Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Float64, xi),
       convert(Float64, eta),
       convert(Array{Float64, 1}, v0),
       v)

   return SVector{3}(v)
end


text: println("e")

max_stars_repo_name: kamilrakowski/e

max_issues_repo_name: kamilrakowski/e

max_forks_repo_name: kamilrakowski/e


text:
# -----------------------------------------------------------------------------
# Scaled-diagonally-dominant matrices [JuMP]
# -----------------------------------------------------------------------------

# Ahmadi, A., & Majumdar, A. (2017). Dsos and sdsos optimization: More
# tractable alternatives to sum of squares and semidefinite optimization.
# SIAM Journal on Applied Algebra and Geometry, 3. https://doi.org/
# 10.1137/18M118935X

# JuMP constraints for scaled-diagonal-dominance ------------------------------
struct sdd end
struct sddConstraint <: AbstractConstraint
        mat::AbstractArray
end

JuMP.build_constraint(::Function, m::AbstractArray, ::sdd) = sddConstraint(m)

function JuMP.add_constraint(model::Model, c::sddConstraint, ::String = "")
        X = Matrix{GenericAffExpr{Float64,VariableRef}}(c.mat)
        n = size(X, 1)
        for j in 1:n, i in 1:j-1
                r = @variable(model, [1:3])
                g = [r[1], r[3], √2 * r[2]]
                @constraint(model, g in MOI.RotatedSecondOrderCone(3))
                X[i, i] -= r[1]
                X[i, j] -= r[2]
                X[j, i] -= r[2]
                X[j, j] -= r[3]
        end
        @constraint(model, X .== zeros(size(X)))
end

# JuMP constraints for the SDD family of cones --------------------------------
# SDD(U) := {M ∈ Symmetric | M = U'QU for some sdd matrix Q}.
# One can think of SDD(U) as the set of matrices that are sdd after a change of
# coordinates via the matrix U.
struct SDD
        u::AbstractArray
end
struct SDDConstraint <: AbstractConstraint
        mat::AbstractArray
        u::AbstractArray
end

JuMP.build_constraint(::Function, m::AbstractArray, d::SDD) =
        SDDConstraint(m, d.u)

function JuMP.add_constraint(model::Model, c::SDDConstraint, ::String = "")
        n = size(c.mat, 1)
        Q = @variable(model, [1:n, 1:n], Symmetric)
        @constraint(model, Q in sdd())
        xs = JuMP.vectorize(c.mat - c.u' * Q * c.u, SymmetricMatrixShape(n))
        @constraint(model, xs in MOI.Zeros(length(xs)))
end


text:
"""
    struct SolutionVector{T} <: AbstractVector{T}
Solution vector stored as 3D vector evaluated at nodes/quadrature points.
`A[:,iS,iK]` is the solution evaluated for state iS and element iK
"""
struct SolutionVector{T} <: AbstractVector{T}
    data::Array{T,3}

    function SolutionVector(data::Array{T,3}) where {T}
        new{T}(data)
    end

    function SolutionVector{T}(s1,s2,s3) where {T}
        new{T}(Array{T,3}(undef, s1,s2,s3))
    end
end

Base.similar(u::SolutionVector, ::Type{T}) where {T} = SolutionVector(Base.similar(u.data, T))

"""
    size(A::SolutionVector)
Return a tuple containing the dimensions of underlying operator A
"""
function Base.size(A::SolutionVector)
    (n,ns,ne) = size(A.data)
    return (n*ns*ne,)
end

"""
    true_size(A::SolutionVector)
Return a tuple containing the dimensions of underlying data for A
"""
function true_size(A::SolutionVector)
    return size(A.data)
end
"""
    true_size(A::SolutionVector, i::Integer)
Return the size of underlying data for A in dimension i
"""
function true_size(A::SolutionVector, i::Integer)
    return size(A.data,i)
end

"""
    getindex(A::SolutionVector, i::Integer)
Linear scalar indexing into array: access A[i]
"""
@inline function Base.getindex(A::SolutionVector, i::Integer)
    @boundscheck checkbounds(A.data, i)
    return A.data[i]
end

"""
    setindex!(A::SolutionVector, v, i::Integer)
Linear scalar indexed assignment: set A[i] = v
"""
@inline function Base.setindex!(A::SolutionVector, v, i::Integer)
    @boundscheck checkbounds(A.data, i)
    setindex!(A.data, v, i)
end

(==)(A::SolutionVector, B::SolutionVector) = A.data == B.data
(-)(A::SolutionVector) = SolutionVector(-A.data)
(+)(A::SolutionVector, B::SolutionVector) = SolutionVector(A.data+B.data)
(-)(A::SolutionVector, B::SolutionVector) = SolutionVector(A.data-B.data)
(*)(c::Number, A::SolutionVector) = SolutionVector(c*A.data)
(*)(A::SolutionVector, c::Number) = SolutionVector(c*A.data)
(/)(A::SolutionVector, c::Number) = SolutionVector(A.data / c)

"""
    lp_error!(u::SolutionVector, u_true::SolutionVector, p::Real=2)
Compute the lp-norm error ||u-u_true||_p of a SolutionVector (see: lp_norm).
Reuses the storage in u_true to compute this error.
"""
function lp_error!(u::SolutionVector, u_true::SolutionVector, p::Real=2)
    u_true -= u
    return lp_norm(u_true, p)
end

"""
    lp_norm(u::SolutionVector, p::Real=2)
Compute the lp-norm ||u||_p of a SolutionVector (see: norm).
If u has multiple states, then return array of L2 norms for each state.
"""
function lp_norm(u::SolutionVector, p::Real=2)
    (n,ns,ne) = true_size(u)
    lp_norms = zeros(ns)
    for iS = 1:ns
        @views lp_norms[iS] += LinearAlgebra.norm(vec(u.data[:,iS,:]),p)
    end
    return lp_norms
end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
#-----------------------------------------------------------
#-----------------------------------------------------------
# Define function to access windows registry
#
#File: jWinReg.jl
#Author: Leonardo Cocco
#Creation Date: 11-04-2016
#Last update 21-11-2016
#ver: 0.12
#----------------------------------------------------------
#----------------------------------------------------------

module jWinReg

if is_windows()
	include("wintypes.jl")
	include("commons.jl")

	include("regfunc.jl")
	export regread,
		   regwrite,
		   regcreatekey,
		   regdelete,
		   regkeylist,
		   regvaluelist
end #windows_only

end # module


text:
include("sr.jl")

include("discrete.jl")
#include("continuous.jl")

include("veronese.jl")


text:
struct DenseSplitBatchedCache{Ta,Tb,Tc,Td,Te}
    σr::Tc
    σc::Tc
    out::Tb

    δℒℒ::Td
    δℒr::Te
    δℒc::Te

    θ::Ta
    θ_tmp::Ta
    valid::Bool
end

function cache(l::DenseSplit{Ta,Tb}, arr_T, in_T ,in_sz, batch_sz) where {Ta,Tb}
    c = DenseSplitBatchedCache(
            similar(l.Wr, size(l.Wr,2), batch_sz),
            similar(l.Wr, size(l.Wr,2), batch_sz),
            similar(l.b, size(l.b,1), batch_sz),

            similar(l.Wr, size(l.Wr,1), batch_sz),
            similar(l.Wr, size(l.Wr,2), batch_sz),
            similar(l.Wr, size(l.Wr,2), batch_sz),

            similar(l.b, size(l.b,1), batch_sz),
            similar(l.b, size(l.b,1), batch_sz),
            false)
    return c
end

batch_size(c::DenseSplitBatchedCache) = size(c.out, 2)

function (l::DenseSplit)(c::DenseSplitBatchedCache, (xr, xc))
    # The preallocated caches
    out  = c.out
    θ = c.θ
    θ₂ = c.θ_tmp

    # Store the input to this layer for the backpropagation
    σr = copyto!(c.σr, xr)
    σc = copyto!(c.σc, xc)

    #θ .= net.b .+ net.W * x
    mul!(θ, l.Wr, σr)
    mul!(θ₂, l.Wc, σc)
    θ .+= l.b .+ θ₂

    # Apply the nonlinear function
    out  .= l.σ.(θ)
    return out
end

function backprop(∇, l::DenseSplit, c::DenseSplitBatchedCache, δℒ)
    # The preallocated caches
    θ = c.θ
    δℒℒ = c.δℒℒ

    # Compute the actual sensitivity
    copyto!(δℒℒ, δℒ)
    δℒℒ .*= fwd_der.(l.σ, θ)

    _batched_outer_prod_noconj!(∇.Wr, δℒℒ, c.σr)
    _batched_outer_prod_noconj!(∇.Wc, δℒℒ, c.σc)
    ∇.b .= δℒℒ

    δℒr = mul!(c.δℒr, transpose(l.Wr), δℒℒ)
    δℒc = mul!(c.δℒc, transpose(l.Wc), δℒℒ)
    return (δℒr, δℒc)
end


text:
# Indefinite integral of pdist.pdf*f is clenshaw_asin(cdf(x),clenshaw_coefs)

function clenshaw_coefs{N}(pdist::Distribution, f::Function, degreetype::Val{N} = Val{50}())
    f1 = Fun(x-> 1/2*cos(x)*f(quantile(pdist, (sin(x)+1)/2)), -π/2..π/2, N-1)
    g1 = cumsum(f1)
    len = min(length(g1.coefficients),N)
    (len < length(g1.coefficients)) && warn("I'm truncating the coefficients!")
    u =  SVector{N, Float64}(vcat(g1.coefficients[1:len],zeros(N-len)))
    return u
end

#clenshaw_asin(x,c) = ApproxFun.clenshaw(2*(asin(2x-1))/π, c)
clenshaw_asin(x,c) = clenshaw_halved(4*(asin(2x-1))/π, c)
@generated function clenshaw_halved{N, R<:Real}(x::R, c::SVector{N})
    Expr(:block,
    :((a,b) = (zero(R), zero(R))),
    (:((a,b) = (muladd(x,a,c[$k]-b),a)) for k = N:-1:2)...,
    :(muladd(x/2,a,c[1]-b)))
end


text:
# this script does not run as a whole
# (they code snippets are demonstrative, and should be used on their own, and/or in the REPL
# on Julia prompt:
# julia> STDIN
# TTY(open, 0 bytes waiting)
# julia> STDOUT
# TTY(open, 0 bytes waiting)

# read(STDIN, Char)
write(STDOUT, "Julia") #> Julia5
# readbytes(STDIN,3)
# readline(STDIN)
# input:  # Julia
# output: # "Julia\r\n" (Windows) - "Julia\n" (Linux)

# stream can be STDIN or any other input stream:
stream = STDIN
for line in eachline(stream)
    print("Found $line")
end

while !eof(stream)
    x = read(stream, Char)
    println("Found: $x")
end

function process(item)
    # do the processing
end

# files:
names(IOStream)
# 4-element Array{Symbol,1}:
#  :handle
#  :ios
#  :name
#  :mark

IOStream.types
# (Ptr{None},Array{Uint8,1},String,Int64)

# reading:
fname = "example.dat"
f1 = open(fname)
# IOStream(<file example.dat>)

data = readlines(f1)
# 3-element Array{Union(ASCIIString,UTF8String),1}:
#  "this is line 1.\r\n"
#  "this is line 2.\r\n"
#  "this is line 3."

for line in data
	println(line)
	process(line)
end
close(f1)

open(fname) do file
    process(file)
end

open(fname) do file
    for line in eachline(file)
        print(line)
        process(line)
    end
end

# writing:
fname = "example2.dat"
f2 = open(fname, "w")
write(f2, "I write myself to a file\n") #> 24
println(f2, "even with println!")
close(f2)

# looping over all files in a folder:
for file in readdir()
	process(file)
end

max_stars_repo_path: Module 1/Chapter08/io.jl

max_issues_repo_path: Module 1/Chapter08/io.jl

max_forks_repo_path: Module 1/Chapter08/io.jl


text:
module Disorderness

	include("randboxcars.jl");
	include("SDE.jl");
	export randboxcars
	export SDE
end


text:
import HDF5
import PyPlot
import SparseArrays
import .MultipeakFunctions
import .MasslistFunctions
import .ResultFileFunctions
import .InterpolationFunctions

#include("masslistFunctions.jl")

function deconvolute(
  filepath;
  outputfilename="results/_result.hdf5",
  binWidth = 6,
  calcTransposed = false
  )
  #binWidth += 1 # Maybe crosstalk is calculated from "<" while summing is done from "<=" ?? Gives better crostalk removal.
  file = joinpath(filepath, outputfilename)

  peakShapesY = HDF5.h5read(file, "MassDepPeakshape")
  peakShapesY = peakShapesY ./ sum(peakShapesY,dims=1) # Normalize!

  peakShapesCenterMass = HDF5.h5read(file, "MassDepPeakshapeCenterMasses")
  totalAvgSpectrum = HDF5.h5read(file, "AvgSpectrum") - HDF5.h5read(file,"AvgBaseline")
  massAxis = HDF5.h5read(file, "MassAxis")
  masslistElements = HDF5.h5read(file,"ElementNames")
  compositionsOrig = HDF5.h5read(file, "ElementalCompositions")
  massesOrig = HDF5.h5read(file, "MassList")
  massBordersLowOrig = HDF5.h5read(file, "MassListIntegrationBordersLow")
  massBordersHighOrig = HDF5.h5read(file, "MassListIntegrationBordersHigh")
  massCenterIdxOrig = HDF5.h5read(file, "MassListIdx")
  massLowIdxOrig = HDF5.h5read(file, "MassListIntegrationBordersIdxLow")
  massHighIdxOrig = HDF5.h5read(file, "MassListIntegrationBordersIdxHigh")
  massScaleMode = HDF5.h5read(file, "MassCalibMode")
  massScaleParameters = HDF5.h5read(file, "MassCalibParameters")
  if unique(massesOrig) != massesOrig
    println("Multiple entries of the same mass --> will produce singular matrix!!!")
  end

  PyPlot.figure()
  ax = PyPlot.subplot(111)

    selector = (massesOrig .> 0)
    masses = massesOrig[selector]
    massBordersLow = massBordersLowOrig[selector]
    massBordersHigh = massBordersHighOrig[selector]
    massCenterIdx = massCenterIdxOrig[selector]
    massLowIdx = massLowIdxOrig[selector]
    massHighIdx = massHighIdxOrig[selector]

    compositions = compositionsOrig[:,selector]

    print("Populating matrix for inversion of linear system...")
    mtrx = MultipeakFunctions.calculateCrossTalkMatrix(masses, massCenterIdx, massLowIdx, massHighIdx, massScaleMode, massScaleParameters, compositions, peakShapesCenterMass, peakShapesY)
    stickRaw = [InterpolationFunctions.interpolatedSum(massLowIdx[i], massHighIdx[i], totalAvgSpectrum) for i=1:length(masses)]
    println(" DONE")

    print("Inverting Matrix...")
    deconvolutionMatrix = SparseArrays.sparse(inv(mtrx)) # sparse on one processor is as fast as dense on mp, --> choosing less power consumption.
    #deconvolutionMatrix = (inv(mtrx)) # there seems to be no multithread support for sparse matrix multiplication, still sparse is same speed
    println(" DONE")

    print("Applying deconvolution kernel...")
    counts = deconvolutionMatrix * stickRaw
    println(" DONE")

  print("Reconstructing Spectrum for visual check...")

  PyPlot.plot(massAxis,totalAvgSpectrum, "-o", label="Original", color="r")
  #reconstructedSpectrum = reconstructSpectrum(massAxis, masses[(masses.>158) & (masses.<162)], masslistElements, compositions[:,(masses.>158) & (masses.<162)], counts[(masses.>158) & (masses.<162)], peakShapesCenterMass, peakShapesY)
  reconstructedSpectrum = MultipeakFunctions.reconstructSpectrum(massAxis, massScaleMode, massScaleParameters, masses, compositions, counts, peakShapesCenterMass, peakShapesY)

  PyPlot.plot(massAxis, reconstructedSpectrum, label="Fit", color="b")
  PyPlot.plot(massAxis, totalAvgSpectrum-reconstructedSpectrum, label="Residual", color="g")
  assyErrorX = [(masses-massBordersLow)'; (massBordersHigh-masses)']
  y = InterpolationFunctions.interpolate(masses, massAxis, totalAvgSpectrum)
  PyPlot.errorbar(TOFFunctions.timebin2mass(massCenterIdx, massScaleMode, massScaleParameters),y,xerr=assyErrorX, fmt="o")
  PyPlot.errorbar(masses,y,xerr=assyErrorX, fmt="x")
  #=
  fittedPeaks = Array{Float64}(length(masses),2001)
  for i=1:length(masses)
    approxMassIndex = searchsortedfirst(massAxis,masses[i])
    fittedPeaks[i,:] = reconstructSpectrum(massAxis[approxMassIndex-300 : approxMassIndex+1700], masses[i], masslistElements, compositions[:,i], counts[i], peakShapesCenterMass, peakShapesY)
    plot(massAxis[approxMassIndex-300 : approxMassIndex+1700],fittedPeaks[i,:],"--", color="green")
  end
  =#

  PyPlot.legend()
  ax[:set_ylim]([minimum(totalAvgSpectrum),maximum(totalAvgSpectrum)])
  println(" DONE")

  ############ WRITE OUTPUT TO FILE ##############################################
  haveStickCps = false
  # Correct timetraces

  fh = HDF5.h5open(file,"r+")

  if Base.haskey(fh, "CorrStickCps")
  HDF5.o_delete(fh,"CorrStickCps")
  end

  if Base.haskey(fh, "CorrStickCpsErrors")
  HDF5.o_delete(fh,"CorrStickCpsErrors")
  end

  if Base.haskey(fh, "StickCps")
    haveStickCps = true

    # Create empty Dataspace
    nbrSpectra = ResultFileFunctions.getNbrTraceSamples(file)
    dset = HDF5.d_create(fh, "CorrStickCps", HDF5.datatype(Float32), HDF5.dataspace(nbrSpectra, length(masses)), "chunk", (1,length(masses)), "compress", 3)

    toProcessLow = 0
    toProcessHigh = 0
    while toProcessHigh < nbrSpectra
      toProcessLow = toProcessHigh + 1
      toProcessHigh = toProcessLow + 9999
      if toProcessHigh > nbrSpectra
        toProcessHigh = nbrSpectra
      end
      print("Correcting spectrum $toProcessLow to $toProcessHigh of $nbrSpectra: Loading...")
      #samplesSubRange = convert(SharedArray, ResultFileFunctions.getTraceSamples(file,toProcessLow:toProcessHigh, raw=true)[:,selector])
      samplesSubRange = ResultFileFunctions.getTraceSamples(file,toProcessLow:toProcessHigh, raw=true)[:,selector]
      print("Deconvoluting...")
      for i=1:(toProcessHigh - toProcessLow + 1)
        #traces[i,:] = deconvolutionMatrix * traces[i,:]
        #tracesErrors[i,:] = abs(deconvolutionMatrix) * sqrt(abs(traces[i,:])/5
        dset[toProcessLow - 1 + i,:] = deconvolutionMatrix *samplesSubRange[i,:]
      end
      println("DONE")
    end
    #HDF5.h5write(file, "CorrStickCps", convert(Array,traces))
    #HDF5.h5write(file, "CorrStickCpsErrors", tracesErrors)
  end
  HDF5.close(fh)

  if calcTransposed
      ResultFileFunctions.transposeStickCps(file)
  end

  fh = HDF5.h5open(file,"r+")
  if Base.haskey(fh, "AvgStickCps")

    traces = HDF5.h5read(file, "AvgStickCps")[:,selector]
    tracesErrors = similar(traces)
    for i=1:size(traces,1)
      #traces[i,:] = deconvolutionMatrix \ traces[i,:]
      traces[i,:] = deconvolutionMatrix * traces[i,:]
      #tracesErrors[i,:] = abs(deconvolutionMatrix) * sqrt(abs(traces[i,:])/3600)
  end

  if Base.haskey(fh, "CorrAvgStickCps")
  HDF5.o_delete(fh,"CorrAvgStickCps")
  end
  if Base.haskey(fh, "CorrAvgStickCpsErrors")
  HDF5.o_delete(fh,"CorrAvgStickCpsErrors")
  end
  HDF5.h5write(file, "CorrAvgStickCps", traces)
  HDF5.h5write(file, "CorrAvgStickCpsErrors", tracesErrors)
  end
  HDF5.close(fh)

  #if haveStickCps
  #  traces = HDF5.h5read(file, "CorrStickCps")[:,selector]
  #else
  #  traces = HDF5.h5read(file, "CorrAvgStickCps")[:,selector]
  #end
  return deconvolutionMatrix
end


text:
import Base: sin, cos, tan, cot, sec, csc
import Base: real, imag, reim, ==
import Base: asin, acos, atan
import Base: asinh, acosh, atanh
import Base: sinh, cosh, tanh, coth, sech, csch
import Base: abs, exp, log, zero
import Base: +, -, *, /, ^

## By constructing bicomplex, we can represent multi-complex numbers using recursion
struct bicomplex
    i1
    i2
end

"""
Using matrix to represent bicomplex.
"""
function mat(bi::bicomplex)
    i1, i2 = bi.i1, bi.i2
    return [i1 -i2; i2 i1]
end

function mat2bicomplex(matrix)
    out1, out2 = matrix[1, 1], matrix[2, 1]
    return bicomplex(out1, out2)
end

reim(bi::bicomplex) = (bi.i1, bi.i2)

==(a::bicomplex, b::bicomplex) = a.i1 == b.i1 && a.i2 == b.i2

## Basic operators: +, -, *, /, ^

+(bi1::bicomplex, bi2::bicomplex) = mat2bicomplex(mat(bi1)+mat(bi2))
+(bi::bicomplex, n) = mat2bicomplex(n.+mat(bi))
+(n, bi::bicomplex) = +(bi, n)

-(bi1::bicomplex, bi2::bicomplex) = mat2bicomplex(mat(bi1)-mat(bi2))
-(bi::bicomplex, n) = mat2bicomplex(mat(bi).-n)
-(n, bi::bicomplex) = -(bi, n)
-(bi::bicomplex) = mat2bicomplex(-mat(bi))

*(bi1::bicomplex, bi2::bicomplex) = mat2bicomplex(mat(bi1)*mat(bi2))
*(n, bi::bicomplex) = mat2bicomplex(n*mat(bi))
*(bi::bicomplex, n) = *(n, bi)

/(bi1::bicomplex, bi2::bicomplex) = mat2bicomplex(mat(bi1)/mat(bi2))
/(n, bi::bicomplex) = mat2bicomplex(n*inv(mat(bi)))
/(bi::bicomplex, n) = mat2bicomplex(mat(bi)/n)

^(bi::bicomplex, power) = mat2bicomplex(mat(bi)^power)

zero(bi::bicomplex) = bicomplex(zero(bi.i1), zero(bi.i2))

## Basic functions: abs, exp, log

abs(bi::bicomplex) = sqrt(bi.i1^2 + bi.i2^2)

function exp(bi::bicomplex)
    i1, i2 = bi.i1, bi.i2
    out1 = exp(i1)*cos(i2)
    out2 = exp(i1)*sin(i2)
    return bicomplex(out1, out2)
end

function log(bi::bicomplex)
    out1 = log(abs(bi))
    out2 = angle(bi)
    return bicomplex(out1, out2)
end

angle(bi::bicomplex) = atan(bi.i2/bi.i1)

## Basic math functions: sin, cos, tan, cot, sec, csc

function sin(bi::bicomplex)
    i1, i2 = bi.i1, bi.i2
    out1 = cosh(i2)*sin(i1)
    out2 = sinh(i2)*cos(i1)
    return bicomplex(out1, out2)
end

function cos(bi::bicomplex)
    i1, i2 = bi.i1, bi.i2
    out1 = cosh(i2)*cos(i1)
    out2 = -sinh(i2)*sin(i1)
    return bicomplex(out1, out2)
end

tan(bi::bicomplex) = sin(bi)/cos(bi)
cot(bi::bicomplex) = cos(bi)/sin(bi)
sec(bi::bicomplex) = 1/cos(bi)
csc(bi::bicomplex) = 1/sin(bi)

## Basic inverse hyperbolic trig function: asinh, acosh, atanh, acoth, asech, acsch

function atan(bi::bicomplex)
    i1, i2 = bi.i1, bi.i2
    temp = atanh(bicomplex(-i2, i1))
    return bicomplex(temp.i2, temp.i1)
end

## Basic hyperbolic functions: sinh, cosh, tanh, coth, sech, csch

function sinh(bi::bicomplex)
    i1, i2 = bi.i1, bi.i2
    out1 = cosh(i1)*cos(i2)
    out2 = sinh(i1)*sin(i2)
    return bicomplex(out1, out2)
end

function cosh(bi::bicomplex)
    i1, i2 = bi.i1, bi.i2
    out1 = sinh(i1)*cos(i2)
    out2 = cosh(i1)*sin(i2)
    return bicomplex(out1, out2)
end

tanh(bi::bicomplex) = sinh(bi)/cosh(bi)
coth(bi::bicomplex) = cosh(bi)/sinh(bi)
sech(bi::bicomplex) = 1/cosh(bi)
csch(bi::bicomplex) = 1/sinh(bi)

image2(bi::bicomplex) = imag(bi.i2)

biderivative(f, point, h) = image2(f(bicomplex(point+im*h, h)))/h^2


text:
using OptiMo
using LinearAlgebra: norm
using Printf

"""
    restore( prob, [x0, with_indicator, verbose, ...] )
"""
function restore(
    prob::TP;
    x0::Tx = prob.meta.x0,
    with_indicator::Bool = false,
    verbose::Bool = false,
    max_out_iter::Int = 100,
    max_sub_iter::Int = 1000,
    ifactor::Real = 0.1,
    wfactor::Real = 0.1,
    wmin::Real = 1e-16,
    tol_cviol::Real = 1e-8,
    tol_inf_cviol::Real = 1e-5,
    tol_infst::Real = 1e-10,
    tol_optim::Real = 1e-14,
    loss::Symbol = :l2sq,
    huber_rho::Real = 1.0,
    huber_mu::Real = 1.0,
) where {TP<:AbstractOptiModel,Tx<:AbstractVector}

    @assert 0 < tol_optim
    @assert 0 < tol_cviol
    @assert 0 < tol_inf_cviol
    @assert 0 < tol_infst
    @assert 0 < wmin
    @assert 0 < wfactor < 1
    @assert 0 < ifactor
    @assert 0 < max_sub_iter
    @assert 0 < max_out_iter

    solver = Bazinga.ZEROFPR( verbose = false,
                              max_iter = max_sub_iter,
                              tol_optim = tol_optim)

    # initial point = base point
    x = copy(x0)
    R = eltype(x)
    # allocation
    s = similar(x)
    infst = -1.0
    if with_indicator
        bigeps = 1e-10
    end
    # diagonal scaling
    d0 = min.(1.0, 1.0 ./ abs.(x0))
    # initial infeasibility
    fprob = FeasOptiModel(
        prob,
        x0 = x,
        with_indicator = with_indicator,
        loss = loss,
        huber_rho = huber_rho,
        huber_mu = huber_mu,
    )
    infea = infeasibility(fprob, x)
    cviol = cviolation(fprob, x)
    # initial weight estimate
    w = max(wmin, ifactor * infea)

    if verbose
        @printf(
            "%4s | %7s | %4s  %7s | %7s  %7s | %7s | \n",
            "o.it",
            "w",
            "i.it",
            "optim",
            "infea",
            "pdist",
            "cviol"
        )
    end

    niter = 0
    while true

        fprob = FeasOptiModel(
            prob,
            x0 = x,
            with_indicator = with_indicator,
            xprox = x0,
            wprox = w,
            dprox = d0,
            loss = loss,
            huber_rho = huber_rho,
            huber_mu = huber_mu,
        )
        out = solver(fprob)
        copyto!(x, out.x)
        infea = infeasibility(fprob, x)
        cviol = cviolation(fprob, x)
        w = max(wmin, min(wfactor * w, ifactor * infea))

        # infeasibility stationarity
        infeasibilitygrad!(fprob, x, s)
        if with_indicator
            s .= x .- bigeps .* s
            prox!(fprob, s, bigeps, s)
            s .= (x .- s) ./ bigeps
        end
        infst = norm(s, Inf)

        if verbose
            proxd = proxdistance(fprob, x)
            @printf(
                "%4d | %.1e | %4d  %.1e | %.1e  %.1e | %.1e | %.1e \n",
                niter,
                w,
                out.iterations,
                out.optimality,
                infea,
                proxd,
                cviol,
                infst
            )
        end

        niter += 1
        if cviol <= tol_cviol
            if verbose
                @printf "cviolation satisfied"
            end
            break
        elseif cviol >= tol_inf_cviol && infst <= tol_infst
            if verbose
                @printf "local infeasibility detected"
            end
            break
        elseif w <= wmin && out.iterations == 1
            if verbose
                @printf "tired restoration"
            end
            break
        elseif niter >= max_out_iter
            if verbose
                @printf "max iterations reached"
            end
            break
        end

    end
    return x, cviol, infst
end


text:
using Revise
# using Profile
# using ProfileView
using RLInterface
using POMDPs
using POMDPModels
using Test
using Random
# using BenchmarkTools
# using POMDPModelTools

function sim(env, nsteps=100, rng=MersenneTwister(0))
    o = reset!(env)
    step = 1
    done = false
    r_tot = 0.0
    na = length(actions(env))
    dims = obs_dimensions(env)
    while !done && step <= nsteps
        action = sample_action(env)
        obs, rew, done, info = step!(env, action)
        r_tot += rew
        step += 1
    end
    return r_tot
end

env = MDPEnvironment(SimpleGridWorld())
# env = POMDPEnvironment(TigerPOMDP())

sim(env)

@inferred sim(env)

@btime sim($env, $4)

@code_warntype sim(env)

@inferred reset!(env)

@inferred step!(env, :up)

@code_warntype step!(env, :up)

@code_warntype gen(DDNOut(:sp, :r), env.problem, env.state, :up, env.rng)

@inferred gen(DDNOut(:sp, :r), env.problem, env.state, :up, env.rng)

function teststep(env, a)
    s, r = gen(DDNOut(:sp, :r), env.problem, env.state, :up, env.rng)
    return (s, r)
end

@inferred teststep(env, :up)

@code_warntype teststep(env, :up)

using POMDPs, POMDPModels, Random, Test

mdp = SimpleGridWorld()
rng = MersenneTwister(1)
s = initialstate(mdp, rng)
a = :up

@inferred gen(DDNOut(:sp, :r), mdp, s, a, rng)

struct MDPEnv{M, S, R}
    prob::M
    s::S
    rng::R
end

function step(env::MDPEnv, a)
    return gen(DDNOut(:sp, :r), env.prob, env.s, a, env.rng)
end

env = MDPEnv(mdp, s, rng)

@inferred step(env, a)

@code_warntype step(env, a)

# Profile.clear()
# @profile for i=1:1000; sim(env); end

# ProfileView.view()


text:
macro tensors(x...)
    return :(nGraph.Tensor.(($(esc(x[1])),), $(esc(x[2]))))
end

@testset "Ops" begin
    # NOTE: Tests are not in alphabetical order because we need the existence of some
    # ops to bootstrap the testing of other ops.
    #
    # I've added tests here in the order I created them, so the ordering should reflect
    # the bootstrap requirement order.

    backend = nGraph.Backend{:CPU}()

    @testset "Add" begin
        x = randn(Float32, 10, 10)
        y = randn(Float32, 10, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)

        Z = X + Y

        # Compile and run.
        ex = nGraph.compile(backend, [X,Y], [Z])
        vX, vY, vZ = @tensors backend (x,y,Z)
        ex([vX, vY], [vZ])

        @test isapprox(parent(vZ), x .+ y)
    end

    @testset "AvgPool" begin
        x = randn(Float32, 30, 30, 10, 10)

        X = nGraph.Node(x)
        Z = nGraph.avgpool(X, (3,3); pad = 0, stride = (3,3))
        ex = nGraph.compile(backend, [X], [Z])
        vX, vZ = @tensors backend (x,Z)
        ex([vX], [vZ])

        # Do an equivalent Flux MeanPool
        m = Flux.MeanPool((3,3))
        @test isapprox(parent(vZ), m(x))
    end

    # Test that constants get slurped up by nGraph
    # TODO: Run over all element types exported by nGraph.
    @testset "Constant" begin
        x = randn(Float32, 10, 10)
        y = randn(Float32, 10, 10)

        X = nGraph.constant(x)
        Y = nGraph.Node(y)
        Z = X + Y
        ex = nGraph.compile(backend, [Y], [Z])
        vY, vZ = @tensors backend (y,Z)

        ex([vY], [vZ])
        @test isapprox(parent(vZ), x .+ y)
    end

    @testset "Broadcast" begin
        x = randn(Float32, 10)
        dims = (10, 8)

        X = nGraph.Node(x)
        Z = nGraph.broadcast(X, dims)
        @test isa(Z, nGraph.Node{Float32,2})
        @test size(Z) == dims

        ex = nGraph.compile(backend, [X], [Z])
        vX, vZ = @tensors backend (x,Z)
        ex([vX], [vZ])

        # Broadcast `x` into `z` for comparison.
        z = zeros(Float32, dims)
        z .= x

        @test isapprox(parent(vZ), z)
    end

    @testset "Convert Eltype" begin
        x = rand(Int32, 10)

        X = nGraph.Node(x)
        Z = nGraph.convert_eltype(Int64, X)
        ex = nGraph.compile(backend, [X], [Z])
        vX, vZ = @tensors backend (x,Z)
        ex([vX], [vZ])

        @test parent(vX) == parent(vZ)
    end

    @testset "Divide" begin
        x = randn(Float32, 10, 10)
        y = randn(Float32, 10, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)

        Z = X ./ Y
        ex = nGraph.compile(backend, [X,Y], [Z])
        vX, vY, vZ = @tensors backend (x,y,Z)
        ex([vX,vY], [vZ])

        @test isapprox(parent(vZ), x ./ y)
    end

    # @testset "Concat" begin
    #     x = randn(Float32, 10, 1)
    #     y = randn(Float32, 10, 1)
    #     z = randn(Float32, 10, 1)

    #     X,Y,Z = nGraph.Node.((x,y,z))

    #     # Concat along dimensions 1 and 2
    #     A = cat(X,Y,Z; dims = 1)
    #     @test size(A) == (30,1)
    #     ex = nGraph.compile(backend, [X,Y,Z], [A])
    #     vX,vY,vZ,vA = @tensors backend (x,y,z,A)
    #     ex([vX,vY,vZ], [vA])
    #     @test parent(vA) == cat(x,y,z; dims = 1)

    #     B = cat(X,Y,Z; dims = 2)
    #     @test size(B) == (10,3)
    #     ex = nGraph.compile(backend, [X,Y,Z], [B])
    #     vX,vY,vZ,vB = @tensors backend (x,y,z,B)
    #     ex([vX,vY,vZ], [vB])
    #     @test parent(vB) == cat(x,y,z; dims = 2)
    # end

    @testset "Dot" begin
        # 1x1
        x = randn(Float32, 10, 1)
        y = randn(Float32, 1, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)
        Z = X * Y
        ex = nGraph.compile(backend, [X,Y], [Z])
        tX,tY,tZ = @tensors backend (x,y,Z)
        ex([tX,tY], [tZ])

        @test isapprox(parent(tZ), x * y)

        # 2x2
        x = randn(Float32, 10, 10)
        y = randn(Float32, 10, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)
        Z = X * Y
        ex = nGraph.compile(backend, [X,Y], [Z])
        tX,tY,tZ = @tensors backend (x,y,Z)
        ex([tX,tY], [tZ])

        @test isapprox(parent(tZ), x * y)

        # Mixed
        x = randn(Float32, 10, 10)
        y = randn(Float32, 10, 1)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)
        Z = X * Y
        ex = nGraph.compile(backend, [X,Y], [Z])
        tX,tY,tZ = @tensors backend (x,y,Z)
        ex([tX,tY], [tZ])
        @test isapprox(parent(tZ), x * y)

        x = randn(Float32, 1, 10)
        y = randn(Float32, 10, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)
        Z = X * Y
        ex = nGraph.compile(backend, [X,Y], [Z])
        tX,tY,tZ = @tensors backend (x,y,Z)
        ex([tX,tY], [tZ])
        @test isapprox(parent(tZ), x * y)
    end

    @testset "GOE" begin
    end

    #####
    ##### Misc Element-Wise + Binary Ops
    #####

    @testset "Log" begin
        # Add 1 to make sure the log is more or less nicely behaved.
        x = rand(Float64, 10, 10) .+ 1
        X = nGraph.Node(x)

        Z = log.(X)
        ex = nGraph.compile(backend, [X], [Z])
        tX,tZ = @tensors backend (x,Z)
        ex([tX],[tZ])
        @test isapprox(parent(tZ), log.(x))
    end

    @testset "Negative" begin
        # Add 1 to make sure the log is more or less nicely behaved.
        x = randn(Float64, 10, 10)
        X = nGraph.Node(x)

        Z = (-).(X)
        ex = nGraph.compile(backend, [X], [Z])
        tX,tZ = @tensors backend (x,Z)
        ex([tX],[tZ])
        @test isapprox(parent(tZ), (-).(x))
    end

    @testset "Maximum" begin
        x = randn(Float64, 10, 10)
        y = randn(Float64, 10, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)

        Z = max.(X,Y)
        ex = nGraph.compile(backend, [X,Y], [Z])
        tX, tY, tZ = @tensors backend (x,y,Z)
        ex([tX, tY],[tZ])
        @test isapprox(parent(tZ), max.(x,y))
    end

    @testset "Minimum" begin
        x = randn(Float64, 10, 10)
        y = randn(Float64, 10, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)

        Z = min.(X,Y)
        ex = nGraph.compile(backend, [X,Y], [Z])
        tX, tY, tZ = @tensors backend (x,y,Z)
        ex([tX, tY],[tZ])
        @test isapprox(parent(tZ), min.(x,y))
    end

    @testset "Multiply" begin
        x = randn(Float64, 10, 10)
        y = randn(Float64, 10, 10)

        X = nGraph.Node(x)
        Y = nGraph.Node(y)

        Z = X .* Y
        ex = nGraph.compile(backend, [X,Y], [Z])
        tX, tY, tZ = @tensors backend (x,y,Z)
        ex([tX, tY],[tZ])
        @test isapprox(parent(tZ), x .* y)
    end
end

# @testset "Reshape" begin
#     backend = nGraph.Backend()
#
#     x = rand(Float32, 100)
#     f = nGraph.compile(backend, x -> reshape(x, 1, :), x)
#     Z = f()
#     @test reshape(x, 1, :) == Z.base
#
#     x = rand(Float32, 2, 2, 2)
#     g = x -> reshape(x, :)
#     f = nGraph.compile(backend, g, x)
#     @test g(x) == f().base
#
#     x = rand(Float32, 3, 2, 1)
#     g = x -> reshape(x, 1, 2, 3)
#     f = nGraph.compile(backend, g, x)
#     @test g(x) == f().base
#
#     # More extravagent reshape
#     x = rand(Float32, 1, 2, 3, 4, 5, 6)
#     g = x -> reshape(x, 6, 5, :, 3, 2)
#
#     N = nGraph.Node(x)
#     M = g(N)
#     @test size(M) == (6, 5, 4, 3, 2)
#     f = nGraph.compile(backend, g, x)
#
#     @test g(x) == f().base
# end
#
# @testset "Softmax" begin
#     backend = nGraph.Backend()
#
#     # 1D case
#     x = rand(Float32, 100)
#     z = softmax(x)
#     f = nGraph.compile(backend, softmax, x)
#     @test isapprox(z, f().base)
#
#     # 2D case
#     x = rand(Float32, 100, 100)
#     z = softmax(x)
#     f = nGraph.compile(backend, softmax, x)
#     @test isapprox(z, f().base)
# end


text:
include("plotting.jl")
include("algConnectivity.jl")
include("growOptimalACgraph.jl")
include("utilities.jl")
fldrs = "empirical_results_2016_04_16_08_57_32__1_"
fname = "final"
num_fldrs = 5
empirical_opt_connectivities = zeros(num_fldrs)
random_connectivities = similar(empirical_opt_connectivities)
opt_connectivity=0
n = 0
edg = 0
for i = 1:num_fldrs
    out_fldr = string(fldrs,i)
    f = open(string(joinpath(out_fldr,fname),".dot"),"r")

    g = SimpleGraphs.SimpleGraph()

    for j in enumerate(eachline(f))
        if j[1] == 1 || eof(f)
            continue
        else
            inew = [parse(Int,x) for x in split(j[2],"--")]
            if length(inew) == 1
                add!(g,inew[1])
            else
                add!(g,inew[1],inew[2])
            end
        end
    end
    if i == 1
        n = NV(g)
        edg = NE(g)
        opt_connectivity = growOptimalA(n,edg)[2]
    end
    close(f)
    empirical_opt_connectivities[i] = AlgConnectivity(g,false)
end

# create random graphs and store connectivites
for i = 1:num_fldrs
    g = simple_NE_sgraph(n,edg)
    random_connectivities[i] = AlgConnectivity(g,false)
end
vals = cat(2,empirical_opt_connectivities,random_connectivities)
plot_c_box(vals,n,opt_connectivity,"","test")


text:
function reset_ensemble(e::Bayes_IPM_ensemble)
    new_e=deepcopy(e)
    for i in 1:length(e.models)
        if string(i) in [basename(record.path) for record in e.models]
            new_e.models[i]=e.models[findfirst(isequal(string(i)), [basename(record.path) for record in e.models])]
        else
            new_e.models[i]=e.retained_posterior_samples[findfirst(isequal(string(i)), [basename(record.path) for record in e.retained_posterior_samples])]
        end
    end

    new_e.log_Li=[new_e.log_Li[1]]
    new_e.log_Xi=[new_e.log_Xi[1]]
    new_e.log_wi=[new_e.log_wi[1]]
    new_e.log_Liwi=[new_e.log_Liwi[1]]
    new_e.log_Zi=[new_e.log_Zi[1]]
    new_e.Hi=[new_e.Hi[1]]

    new_e.sample_posterior && (new_e.retained_posterior_samples=Vector{Model_Record}())

    model_counter=length(new_e.models)+1

    return new_e
end


text:
################################################################################
#
# CKKS scheme setup.
#
################################################################################

using SEAL

function setup_ckks(poly_modulus_degree, primes)
    # Set the encryption scheme to CKKS.
    params = EncryptionParameters(SchemeType.ckks)

    # Set the degree of the cyclotomic polynomial to be a power of 2.
    set_poly_modulus_degree!(params, poly_modulus_degree)

    # Set the ciphertext coefficient modulus (a large integer) as the product of
    # large distinct prime numbers.
    set_coeff_modulus!(params, coeff_modulus_create(poly_modulus_degree, primes))

    # Now that all the parameters are set, we are ready to construct a SEALContext
    # object. This is a heavy class that checks the validity and properties of the
    # parameters we just set.
    context = SEALContext(params)

    # Check context parameters:
    print_parameters(context)
    println("Parameter validation: ", parameter_error_message(context))

    # Generate the public and private keys by using the KeyGenerator, which
    # automatically generates the public and private keys for us.
    keygen = KeyGenerator(context)

    pub_key = PublicKey()
    create_public_key!(pub_key, keygen)

    pri_key = secret_key(keygen)

    relin_keys = RelinKeys()
    create_relin_keys!(relin_keys, keygen)

    galois_keys = GaloisKeys()
    create_galois_keys!(galois_keys, keygen)

    # Create an encryptor, evaluator and decryptor.
    encryptor = Encryptor(context, pub_key)
    evaluator = Evaluator(context)
    decryptor = Decryptor(context, pri_key)

    # Create a CKKSEncoder for encoding values and vectors.
    encoder = CKKSEncoder(context)

    (
        context = context,
        public_key = pub_key,
        private_key = pri_key,
        relinearization_keys = relin_keys,
        galois_keys = galois_keys,
        encryptor = encryptor,
        evaluator = evaluator,
        decryptor = decryptor,
        encoder = encoder
    )
end

function print_parameters(context::SEALContext)
    context_data = key_context_data(context)
    encryption_parms = parms(context_data)
    scheme_type = scheme(encryption_parms)
    if scheme_type == SchemeType.ckks
      scheme_name = "CKKS"
    elseif scheme_type == SchemeType.bfv
      scheme_name = "BFV"
    else
      error("unsupported scheme")
    end

    println("/")
    println("| Encryption parameters:")
    println("|   scheme: ", scheme_name)
    println("|   poly_modulus_degree: ", poly_modulus_degree(encryption_parms))

    print("|   coeff_modulus size: ", total_coeff_modulus_bit_count(context_data), " (")
    bit_counts = [bit_count(modulus) for modulus in coeff_modulus(encryption_parms)]
    print(join(bit_counts, " + "))
    println(") bits")

    if scheme_type == SchemeType.ckks
      println("|   plain_modulus: ", value(plain_modulus(encryption_parms)))
    end

    println("\\")
end


text:
f1 = DataFrame(
    A = [false, true, false, true],
    B = [false, false, true, true],
    potential = [0.75, 0.60, 0.25, 0.40]
    ) |> Table

f2 = DataFrame(
    A = [false, true],
    potential = [0.9, 0.1]
    ) |> Table

let
    # factor multiplication
    f12 = f1 * f2
    @test size(f12) == (4,3)
    @test elementwise_isapprox(partialsort(f12, Assignment(:A=>false, :B=>false))[:potential], [0.75*0.9])
    @test elementwise_isapprox(partialsort(f12, Assignment(:A=>true,  :B=>false))[:potential], [0.60*0.1])
    @test elementwise_isapprox(partialsort(f12, Assignment(:A=>false, :B=>true))[:potential], [0.25*0.9])
    @test elementwise_isapprox(partialsort(f12, Assignment(:A=>true,  :B=>true))[:potential], [0.40*0.1])
end

let
    # factor marginalization
    f1_sans_B = sumout(f1, :B)
    @test size(f1_sans_B) == (2,2)
    @test elementwise_isapprox(partialsort(f1_sans_B, Assignment(:A=>false))[:potential], [0.75 + 0.25])
    @test elementwise_isapprox(partialsort(f1_sans_B, Assignment(:A=>true))[:potential], [0.60 + 0.40])

    f1_sans_A = sumout(f1, :A)
    @test size(f1_sans_A) == (2,2)
    @test elementwise_isapprox(partialsort(f1_sans_A, Assignment(:B=>false))[:potential], [0.75 + 0.60])
    @test elementwise_isapprox(partialsort(f1_sans_A, Assignment(:B=>true))[:potential], [0.25 + 0.40])
end

let
    # factor normalization
    f3 = BayesNets.normalize!(Table(DataFrame(
        A = [false, true],
        potential = [1.0, 3.0]
    )))

    @test elementwise_isapprox(f3[:potential], [0.25, 0.75])
end

let
    # estimation
    table = fit(Table, DataFrame(
        A = [false, false, true, true, true]
    ))
    @test elementwise_isapprox(table[:potential], [2/5, 3/5])

    # TODO: properly test this
    # estimate_convergence(df, Assignment(:A=>true))
end


text:
 #module MultilevelGenesis

#############################################
###### Dependencies #########################
#############################################
using Random
#using Plots
using DataFrames
using JSON
using CSV
using Combinatorics

include("Data.jl")
include("Model_CRS.jl")
include("MoBlue.jl")
#############################################
###### Rection Functions ####################
#############################################

function CalculatePropensities(concentrations::Array{Float64,1}, CRS::CRS)
    ########## NOT UNIT TESTED! ##########
    n_rxn = length(CRS.reaction_list)::Int64
    PropensityVec = zeros(n_rxn)
    rID = 0
    reaction_list = CRS.reaction_list::Array{Reaction,1}
    for rxn in CRS.reaction_list
        Ap = 0.0
        if rxn.propensity == "STD"
            Ap = STD_propensity(rxn, concentrations)::Float64
        end
        rID += 1
        #println("reaction_list ID: ", rID, "    Propensity: ", Ap)
        PropensityVec[rID] = Ap
    end

    return PropensityVec
end

function PickReactionID(PropensityVec::Array{Float64,1})
    Ap_tot = sum(PropensityVec)
    dice_roll = Ap_tot*rand()
    checkpoint = 0.0

    rID = 0
    #println("Dice Roll: ", dice_roll)
    while checkpoint < dice_roll
        rID += 1
        checkpoint += PropensityVec[rID]
        #println("rID: ", rID, "     Checkpoint: ", checkpoint)
    end
    #println("Returning: ", rID, "   Rxn Propensity: ", PropensityVec[rID])
   return rID
end

function ExecuteReaction(concentrations::Array{Float64,1}, CRS::CRS, rID::Int64)
    ########## NOT UNIT TESTED! ##########
    rxn = CRS.reaction_list[rID]
    for (i, r) in enumerate(rxn.reactants)
        concentrations[r] -= rxn.react_coef[i]
    end
    for (j, p) in enumerate(rxn.products)
        concentrations[p] += rxn.prod_coef[j]
    end
    if any(x->x<0, concentrations)
        println(CRS.reaction_list[rID])
        for (i, r) in enumerate(rxn.reactants)
            concentrations[r] += rxn.react_coef[i]
        end
        for (j, p) in enumerate(rxn.products)
            concentrations[p] -= rxn.prod_coef[j]
        end
        println("Concentration of Reactants: ",concentrations[rxn.reactants])
        println("Rxn Propensity based on concentrations: ", STD_propensity(rxn, concentrations))
        error("Concentrations went negative on reaction $rID")

    end

    return concentrations
end

function STD_propensity(rxn::Reaction, concentrations::Array{Float64,1})
    ########## NOT UNIT TESTED! ##########
    Ap  = rxn.rate_constant
    for (ri,r) in enumerate(rxn.reactants)
        rconc = concentrations[r]
        r_order = rxn.react_coef[ri]
        Ap = Ap*(prod([(rconc-i) for i in 0:(r_order-1)]))
    end
    cat_enhance = 0.0
    for (ci,c) in enumerate(rxn.catalysts)
        cconc = concentrations[c]
        cat_enhance += cconc*rxn.cat_coef[ci]
    end
    Ap = Ap*(1 + cat_enhance)
    return Ap
end

# function sort_reactions_by_propensity(concentrations, myCRS)
#     println(myCRS.reaction_list[1:15])
#     propensities = CalculatePropensities(concentrations, myCRS)
#     println(propensities)
#     sorted_indices = sortperm(propensities, rev = true)
#     println(sorted_indices)
#     myCRS.reaction_list[:] = myCRS.reaction_list[sorted_indices]
#     println(myCRS.reaction_list[1:15])
#     return myCRS
# end

#############################################
##### Stochastic Evolution Functions ########
#############################################

function Run_Binary_Ligation(nA::Int64, nB::Int64, max_L::Int64, kf::Float64, kb::Float64, volume::Float64, tau_max::Float64, tau_freq::Float64 = 0.1, T::Float64 = 300.0, R::Float64= 8.3144598, mA::Int64 = 100, mB::Int64 = 100)
    ### Unique Experiment number required

    ## Check for the correct output files and directories
    param_file = "../data/Binary_Ligation_Run_Parameters.csv"
    save_dir = "../data/Binary_Ligation_timeseries"

    if !isdir(save_dir)
        mkdir(save_dir)
    end

    ## Initialize Time values
    tau = 0.0
    freq_count = 0.0

    binary_CRS = generate_binary_ligation_CRS(max_L, kf, kb, volume, T)
    ### Redo output initialization using DataFrames
    output_DF = InitializeOutput(binary_CRS)

    nmolecules = length(binary_CRS.molecule_list)
    concentrations= zeros(nmolecules)

    Aindex = binary_CRS.molecule_dict["A"]
    Bindex = binary_CRS.molecule_dict["B"]
    concentrations[Aindex] = nA
    concentrations[Bindex] = nB

    #Initialize Propensities
    propensities = CalculatePropensities(concentrations, binary_CRS)

    ####### Main LOOP #######
    while tau < tau_max
        # Pick Reaction
        rID = PickReactionID(propensities)
        # Execute Reaction
        concentrations = ExecuteReaction(concentrations, binary_CRS, rID)
        # Calculate Propensities
        propensities = CalculatePropensities(concentrations, binary_CRS)
        # Record Data
        if tau >= freq_count
            output_DF[Symbol(freq_count)] = concentrations[:]
            freq_count += tau_freq
            println(tau)
        end

        #Update Time
        Ap_tot = sum(propensities)
        tau -= (log(rand())/Ap_tot)
    end

    # Save time series and params
    save_name, parameter_df = generate_output_data(binary_CRS, save_dir)#
    println(parameter_df)
    if !isfile(param_file)
        CSV.write(param_file, parameter_df)
    else
        runs_df = CSV.read(param_file)
        runs_df = [runs_df; parameter_df]
    end

    CSV.write(save_name, output_DF)
    return concentrations
end

function Run_MoBlue(Mo1_mass::Int64, total_time::Float64, k_f::Float64, stable_backward::Float64, dimerization_ratio::Float64, mo36_enhance::Float64, wheel_enhancement_multiplier::Float64= 1.0,ball_growth_multipler::Float64 = 1.0,out_count::Float64 = 100.0, volume::Float64 = 1.0, T::Float64 = 300.0, R::Float64= 8.3144598)::Array{Float64,1}

    ### Unique Experiment number required
    tau_freq = total_time/out_count
    ## Check for the correct output files and directories
    param_file = "../data/MoBlue_Run_Parameters.csv"
    save_dir = "../data/MoBlue_timeseries"

    if !isdir(save_dir)
        mkdir(save_dir)
    end

    ## Initialize Time values
    tau = 0.0
    freq_count = 0.0

    MoBlueCRS = make_MoBlueCRS(k_f, stable_backward, dimerization_ratio, mo36_enhance,wheel_enhancement_multiplier,ball_growth_multipler, volume, T, R)::CRS
    ### Redo output initialization using DataFrames
    output_DF = InitializeOutput(MoBlueCRS)

    nmolecules = length(MoBlueCRS.molecule_list)::Int64
    concentrations= zeros(nmolecules)::Array{Float64,1}

    Monomer_index = MoBlueCRS.molecule_dict["Mo1"]
    concentrations[Monomer_index] = Mo1_mass

    #Initialize Propensities
    propensities = CalculatePropensities(concentrations, MoBlueCRS)

    ####### Main LOOP #######
    while tau < total_time
        # Pick Reaction
        rID = PickReactionID(propensities)::Int64
        #println("rID is : ", rID, "     PropensityVec: ", propensities[rID] )
        #println("Picked reaction is: ", MoBlueCRS.reaction_list[rID])
        #println("Reactants: ", MoBlueCRS.molecule_list[MoBlueCRS.reaction_list[rID].reactants])
        # println("Products: ", MoBlueCRS.molecule_list[MoBlueCRS.reaction_list[rID].products])
        #println("Concentrations: ", concentrations[MoBlueCRS.reaction_list[rID].reactants])
        #println("Propensity: ", STD_propensity(MoBlueCRS.reaction_list[rID],concentrations))

        # Execute Reaction
        concentrations = ExecuteReaction(concentrations, MoBlueCRS, rID)::Array{Float64,1}

        # Calculate Propensities
        propensities = CalculatePropensities(concentrations, MoBlueCRS)::Array{Float64,1}
        # Record Data
        if tau >= freq_count
            # Calculate mass and check conservation...
            current_mass = calculate_mass(concentrations, MoBlueCRS)::Float64
            if current_mass != Mo1_mass
                println(MoBlueCRS.reaction_list[rID])
                println(MoBlueCRS.molecule_dict[MoBlueCRS.reaction_list[rID].reactants],MoBlueCRS.molecule_list[MoBlueCRS.reaction_list[rID].products] )
                error("Conservation of Mass not maintained.. Exiting Simulation")
            end
            output_DF[Symbol(freq_count)] = concentrations[:]
            freq_count += tau_freq
            #MoBlueCRS = sort_reactions_by_propensity(concentrations, MoBlueCRS)

            sort!(MoBlueCRS.reaction_list, by=x -> STD_propensity(x, concentrations), rev = true)
            propensities = CalculatePropensities(concentrations, MoBlueCRS)::Array{Float64,1}
            #println(tau)
        end

        #Update Time
        Ap_tot = sum(propensities)
        tau -= (log(rand())/Ap_tot)
    end

    # Save time series and params
    save_name, parameter_df = generate_output_data(MoBlueCRS, save_dir)#
    #println(parameter_df)
    if !isfile(param_file)
        #println("Not Finding File")
        CSV.write(param_file, parameter_df)
    else
        CSV.write(param_file, parameter_df; append = true)
    end
    CSV.write(save_name, output_DF)
    return concentrations
end

function calculate_mass(concentrations, CRS)
    mass = 0
    nmolecules= length(CRS.molecule_list)
    for i in 1:nmolecules
        M = count_mo_num(CRS.molecule_list[i])
        mass += M*concentrations[i]
    end
    return mass
end

#############################################
#end #### module #############################
#############################################


text:
module TestModels

  using SearchLight

  ######## Model from Genie-Searchligth-example-app extracted ############
  export Book, BookWithInterns
  using SearchLight, Dates

  ######## Model from Genie-Searchligth-example-app extracted ############
  export Callback
  export seed, fields_to_store

  mutable struct Book <: AbstractModel

    ### FIELDS
    id::DbId
    title::String
    author::String
    cover::String

    ### VALIDATION
    # validator::ModelValidator

    ### CALLBACKS
    # before_save::Function
    # after_save::Function
    # on_save::Function
    # on_find::Function
    # after_find::Function

    ### SCOPES
    # scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

    ### constructor
    Book(;
      ### FIELDS
      id = DbId(),
      title = "",
      author = "",
      cover = "",

      ### VALIDATION
      # validator = ModelValidator([
      #   ValidationRule(:title, BooksValidator.not_empty)
      # ]),

      ### CALLBACKS
      # before_save = (m::Todo) -> begin
      #   @info "Before save"
      # end,
      # after_save = (m::Todo) -> begin
      #   @info "After save"
      # end,
      # on_save = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On save"
      # end,
      # on_find = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On find"
      # end,
      # after_find = (m::Todo) -> begin
      #   @info "After find"
      # end,

      ### SCOPES
      # scopes = Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}()

    ) = new(id, title, author, cover                                                             ### FIELDS
            # validator,                                                                  ### VALIDATION
            # before_save, after_save, on_save, on_find, after_find                       ### CALLBACKS
            # scopes                                                                      ### SCOPES
            )
  end

  mutable struct BookWithInterns <: AbstractModel

    ### FIELDS
    id::DbId
    title::String
    author::String
    cover::String

    ### VALIDATION
    # validator::ModelValidator

    ### CALLBACKS
    # before_save::Function
    # after_save::Function
    # on_save::Function
    # on_find::Function
    # after_find::Function

    ### SCOPES
    # scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

    ### constructor
    BookWithInterns(;
      ### FIELDS
      id = DbId(),
      title = "",
      author = "",
      cover = "",

      ### VALIDATION
      # validator = ModelValidator([
      #   ValidationRule(:title, BooksValidator.not_empty)
      # ]),

      ### CALLBACKS
      # before_save = (m::Todo) -> begin
      #   @info "Before save"
      # end,
      # after_save = (m::Todo) -> begin
      #   @info "After save"
      # end,
      # on_save = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On save"
      # end,
      # on_find = (m::Todo, field::Symbol, value::Any) -> begin
      #   @info "On find"
      # end,
      # after_find = (m::Todo) -> begin
      #   @info "After find"
      # end,

      ### SCOPES
      # scopes = Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}()

    ) = new("bookwithinterns", "id", Symbol[],                                                ### INTERNALS
            id, title, author, cover                                                             ### FIELDS
            # validator,                                                                  ### VALIDATION
            # before_save, after_save, on_save, on_find, after_find                       ### CALLBACKS
            # scopes                                                                      ### SCOPES
            )
  end

  Base.@kwdef mutable struct Callback <: AbstractModel
    id::DbId = DbId()
    title::String = ""
    indicator::Bool = true
    created_at::String = string(Dates.now())
    # callbacks
    before_save::Function = (m::Callback) -> begin
      @info "Do something before saving"
    end
    after_save::Function = (m::Callback) -> begin
      @info "Do something after saving"
    end
  end

  function seed()
    BillGatesBooks = [
      ("The Best We Could Do", "Thi Bui"),
      ("Evicted: Poverty and Profit in the American City", "Matthew Desmond"),
      ("Believe Me: A Memoir of Love, Death, and Jazz Chickens", "Eddie Izzard"),
      ("The Sympathizer!", "Viet Thanh Nguyen"),
      ("Energy and Civilization, A History", "Vaclav Smil")
    ]
  end

end ### End Module

max_stars_repo_name: GenieFramework/SearchLightMySQL

max_issues_repo_name: GenieFramework/SearchLightMySQL

max_forks_repo_name: GenieFramework/SearchLightMySQL


text:
using Pkg
Pkg.activate(".")

using Genie, Genie.Router, Genie.Renderer, Genie.Renderer.Html

form = """
<form action="/" method="POST" enctype="multipart/form-data">
  <input type="text" name="greeting" value="hello genie" />
  <input type="submit" value="post" />
</form>
"""

layout = "<? @yield ?>"

route("/") do
  html(form)
end

route("/", method = POST) do
  @show @params

  @show @params(:greeting)

  @params(:greeting)
end

Genie.AppServer.startup(; open_browser = false, async = false)


text:
polyEval(coeffs, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)])

#Importing parameters
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport parameters_julia

out_dir = parameters_julia.out_dir
files_ID = parameters_julia.files_ID
month_w = parameters_julia.month_w
year = parameters_julia.year
instances = parameters_julia.instances_ID
key = "(6, 1.5, 1000.0, 1)"

cross_validation_best_key_Apr_AM = key
cross_validation_best_key_Apr_MD = key
cross_validation_best_key_Apr_PM = key
cross_validation_best_key_Apr_NT = key

coeffs_dict_Apr_AM = readstring( out_dir * "coeffs_dict_"* month_w * "_" * instances[1] * ".json")
coeffs_dict_Apr_AM = JSON.parse(coeffs_dict_Apr_AM)

coeffs_dict_Apr_MD = readstring( out_dir * "coeffs_dict_"* month_w * "_" * instances[2] * ".json")
coeffs_dict_Apr_MD = JSON.parse(coeffs_dict_Apr_MD)

coeffs_dict_Apr_PM = readstring( out_dir * "coeffs_dict_"* month_w * "_" * instances[3] * ".json")
coeffs_dict_Apr_PM = JSON.parse(coeffs_dict_Apr_PM)

coeffs_dict_Apr_NT = readstring( out_dir  * "coeffs_dict_"* month_w * "_" * instances[4] * ".json")
coeffs_dict_Apr_NT = JSON.parse(coeffs_dict_Apr_NT)

#coeffs_dict_Apr_weekend = readall("../temp_files/coeffs_dict_Apr_weekend.json")
#coeffs_dict_Apr_weekend = JSON.parse(coeffs_dict_Apr_weekend)

coeffs_best_Apr_AM = coeffs_dict_Apr_AM[cross_validation_best_key_Apr_AM]
coeffs_best_Apr_MD = coeffs_dict_Apr_MD[cross_validation_best_key_Apr_MD]
coeffs_best_Apr_PM = coeffs_dict_Apr_PM[cross_validation_best_key_Apr_PM]
coeffs_best_Apr_NT = coeffs_dict_Apr_NT[cross_validation_best_key_Apr_NT]
#coeffs_best_Apr_weekend = coeffs_dict_Apr_weekend[cross_validation_best_key_Apr_weekend]

using PyPlot

xs = linspace(0, 3.2, 20)

#true_coeffs = [1, 0, 0, 0, .15]
#ys_true = map(x->polyEval(true_coeffs, x), xs)
#plot(xs, ys_true, "k", label="True")

ys_Apr_AM = map(x->polyEval(coeffs_best_Apr_AM, x), xs)
ys_Apr_MD = map(x->polyEval(coeffs_best_Apr_MD, x), xs)
ys_Apr_PM = map(x->polyEval(coeffs_best_Apr_PM, x), xs)
ys_Apr_NT = map(x->polyEval(coeffs_best_Apr_NT, x), xs)
#ys_Apr_weekend = map(x->polyEval(coeffs_best_Apr_weekend, x), xs)

plot(xs, ys_Apr_AM, ">-r", label="AM Peak", linewidth=1.5)
plot(xs, ys_Apr_MD, "s-c", label="Middle Day", linewidth=1.5)
plot(xs, ys_Apr_PM, "<-b", label="PM Peak", linewidth=1.5)
plot(xs, ys_Apr_NT, "^-m", label="Night", linewidth=1.5)
#plot(xs, ys_Apr_weekend, "v-g", label="Weekend", linewidth=1.5)

legend(loc="upper left")

grid("on")
#=
xlim(-0.1, 1.25);
ylim(0.98, 1.5);

font1 = ["family"=>"serif",
    "color"=>"black",
    "weight"=>"normal",
    "size"=>12]
xlabel("Scaled Flow", fontdict=font1)
ylabel("Scaled Cost", fontdict=font1)
=#
savefig(out_dir * "fitting_comparison_"* month_w * ".pdf")

max_stars_repo_path: 03. Plot_comparison_results.jl

max_issues_repo_path: 03. Plot_comparison_results.jl

max_forks_repo_path: 03. Plot_comparison_results.jl


text:
# DELETE THIS IT IS OLD !!!
type class_newton_solver <: abstract_newton_solver
	direction::class_variables

	update_system::Function
	compute_direction::Function
	maximum_delta_decrease::Function

	function class_newton_solver(qp::class_quadratic_program, settings::class_settings)
		this = new();

		linear_system_solver = settings.linear_system_solver;
		linear_system_solver.initialize(qp, settings);

		this.direction = class_variables(qp.n,qp.m);
		rhs = zeros(qp.n + qp.m);
		h = zeros(qp.n + qp.m);
		p = zeros(qp.n + qp.m);

		function compute_h_vector(vars::class_variables)
			try
				rhs[1:(qp.n)] = qp.g(vars) - qp.H(vars) * vars.x()/vars.tau();
				rhs[(qp.n + 1):(qp.m + qp.n)] = -qp.b;

				GLOBAL_timer::class_algorithm_timer
				GLOBAL_timer.start("Solve")
				linear_system_solver.ls_solve!(rhs, h, qp);
				GLOBAL_timer.stop("Solve")
			catch e
				println("ERROR in class_newton_solver.compute_h_vector")
				throw(e)
			end
		end

		this.update_system = function(vars::class_variables)
			try
				GLOBAL_timer::class_algorithm_timer
				GLOBAL_timer.start("Factor")
				inertia = linear_system_solver.ls_factor(vars, qp);
				GLOBAL_timer.stop("Factor")

				compute_h_vector(vars) # move ???

				return inertia # is inertia is correct?
			catch e
				println("ERROR in class_newton_solver.update_system")
				throw(e)
			end
		end

		this.compute_direction = function(vars::class_variables, residuals::class_residuals, gamma::Float64)
			try
				eta = 1 - gamma;

				r_G = residuals.r_G;
				r_D = residuals.r_D;
				r_P = residuals.r_P;
				x = vars.x();
				s = vars.s();
				tau = vars.tau();
				kappa = vars.kappa();
				mu = vars.mu();
				g = qp.g(vars);

				xs = -x .* s + gamma * mu * ones(qp.n);
				tk = -tau*kappa + gamma * mu;

				rhs[1:(qp.n)] = eta * r_D + xs ./ x;
				rhs[(qp.n + 1):(qp.m + qp.n)] = eta * r_P;

				GLOBAL_timer::class_algorithm_timer
				GLOBAL_timer.start("Solve")
				linear_system_solver.ls_solve!(rhs, p, qp);
				GLOBAL_timer.stop("Solve")

				v_tmp = [g' + x'/tau * qp.H(vars) qp.b'];
				numerator_dtau = eta*r_G + 1/tau * tk + (v_tmp * p)[1];
				denominator_dtau = kappa/tau + (1/tau)^2 * (x' * qp.H(vars) * x)[1] + (v_tmp * h)[1];
				dtau = numerator_dtau/denominator_dtau;

				dir = this.direction;
				dir.v[dir.tau_ind] = dtau;
				dir.v[1:(qp.n + qp.m)] = p - h * dtau; # (dx,dy)

				dir.v[dir.s_ind] = (xs - s.*dir.x()) ./ x;
				dir.v[dir.kappa_ind] = (tk - kappa*dir.tau()) / tau;

				# is this direction valid ?
				dir.check_direction();
			catch e
				println("ERROR in class_newton_solver.compute_direction")
				throw(e)
			end
		end

		return this;
	end
end

max_forks_repo_name: ohinder/Everest