text:
function ParAmulB!(nthreads::Int64, yt::Array{Float64,1}, y::Array{Float64,1}, A::SparseMatrixCSC{Float64,Int64}, x::Array{Float64,1})
    (m,n) = size(A)
    length(yt) == nthreads*m || throw(DimensionMismatch("yt"))
    if length(x) != n || length(y) != m
       throw(DimensionMismatch("length(x) != n || length(y) != m"))
    end
    p = ccall((:a_mul_b_prr_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                                                    &nthreads,  &m,         &n,         A.nzval,      A.rowval,   A.colptr,   x,            y,            yt           )
    return nothing
end

function AmulB!(y::Array{Float64,1}, A::SparseMatrixCSC{Float64,Int64}, x::Array{Float64,1})
    (m,n) = size(A)
    if length(x) != n || length(y) != m
       throw(DimensionMismatch("length(x) != n || length(y) != m"))
    end
    p = ccall((:a_mul_b_rr_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}),
                                                   &m,         &n,         A.nzval,      A.rowval,   A.colptr,   x,            y             )
    return nothing
end

function AmulB!(y::Array{Complex128,1}, A::SparseMatrixCSC{Float64,Int64}, x::Array{Complex128,1})
    (m,n) = size(A)
    if length(x) != n || length(y) != m
       throw(DimensionMismatch("length(x) != n || length(y) != m"))
    end
    p = ccall((:a_mul_b_rc_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}),
                                                   &m,         &n,         A.nzval,      A.rowval,   A.colptr,   x,               y                )
    return nothing
end

function AmulB!(y::Array{Complex128,1}, A::SparseMatrixCSC{Complex128,Int64}, x::Array{Complex128,1})
    (m,n) = size(A)
    if length(x) != n || length(y) != m
       throw(DimensionMismatch("length(x) != n || length(y) != m"))
    end
    p = ccall((:a_mul_b_cc_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}),
                                                      &m,         &n,         A.nzval,         A.rowval,   A.colptr,   x,               y                )
    return nothing
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
module OptimalPlasticityRatios

using StaticArrays
using Flux
using Random
using LinearAlgebra
using Statistics
using RecipesBase
using Tracker

include("miscellaneous.jl")
include("updateInfo.jl")
include("trainingInfo.jl")
include("lossInfo.jl")
include("recordingInfo.jl")
include("updateMethods.jl")
include("trainingLoop.jl")
include("dataParsing.jl")

export training_loop!, get_model, recorder, extract_endmean

end # module

ext: jl


text:
using ArrayFire
using Test
using Libdl,Random,SparseArrays,LinearAlgebra,Statistics
using FFTW

@testset "Main" begin
    @testset "Bugs" begin
        include("bugs.jl")
    end

    allowslow(AFArray, false)

    @testset "FFT" begin
        include("fft.jl")
    end

    allowslow(AFArray) do
        include("indexing.jl")
    end
    include("sparse.jl")
    include("math.jl")
    include("blackscholes.jl")
    include("array.jl")
end


text:
using Base.Test
using Sigma
import Sigma: measure, logmeasure

# What should a domain be
b = HyperBox([0.0 0.0
              2.0 1.0])
@test ndims(b) == 2
@test measure(b) == 2
@test logmeasure(b) == log(2) + 0


text:
module CIS3223

include("types.jl")
export
    Graph

include("utils.jl")

include("chapter0.jl")
include("chapter1.jl")
include("chapter2.jl")
include("chapter3.jl")
include("chapter4.jl")
include("chapter5.jl")
include("chapter6.jl")
include("chapter7.jl")
include("chapter8.jl")
include("chapter9.jl")
include("chapter10.jl")

using .Chapter0
export
    fib1,
    fib2

using .Chapter1
export
    multiply,
    divide,
    modexp,
    euclid,
    extendedeuclid,
    primality,
    primality2

using .Chapter2
export
    karatsuba,
    karatsuba2,
    mergesort,
    iterativemergesort,
    polynomialmultiplication

using .Chapter3
export
    explore,
    dfs

using .Chapter4
export
    dijkstra

using .Chapter5
using .Chapter6
export
    lcs

using .Chapter7
using .Chapter8
using .Chapter9
using .Chapter10

end # module


text:
include("MakeData.jl")
include("Transform.jl")
include("Train.jl")
include("Analyze.jl")

function RunProject()

# generate the raw training data
MakeData()

# transform the raw statistics, and split out params and stats
info = Transform()
# when this is done, can delete raw_data.bson

# train the net using the transformed training/testing data
Train()
include("Importance.jl")

end
RunProject()


text:
using MyWeb
using Documenter
using DocumenterTools: Themes

# Generate Sass files.
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "style.scss"), String)
    theme = read(joinpath(@__DIR__, "$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "$(w).scss"), header*"\n"*theme)
end

# Compile Sass files.
Themes.compile(joinpath(@__DIR__, "light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

makedocs(;
    modules=[MyWeb],
    repo="https://github.com/Shoichiro-Tsutsui/MyWeb.jl/blob/{commit}{path}#L{line}",
    sitename = "Shoichiro Tsutsui's Web Page",
    authors = "Shoichiro Tsutsui",
    doctest = false,
    linkcheck = false,
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Shoichiro-Tsutsui.github.io/MyWeb.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "CV" => "cv.md",
        "Publications" => "publications.md",
        "Presentations" => "presentations.md",
        "Links" => "links.md"
    ],
)

deploydocs(;
    repo="github.com/Shoichiro-Tsutsui/MyWeb.jl",
)

max_stars_repo_name: Shoichiro-Tsutsui/MyWeb

max_issues_repo_name: Shoichiro-Tsutsui/MyWeb

max_forks_repo_name: Shoichiro-Tsutsui/MyWeb


text:
#-----------------------------------------------------------------------
#		PLOTS AND RESULTS
#-----------------------------------------------------------------------

using Plots, LinearAlgebra, Serialization, BenchmarkTools

include("../src/squirrel.jl")
include("../src/metric.jl")
include("../src/type.jl")
include("plotfunc.jl")

#-----------------------------------------------------------------------
#       PLOT PARAMETERS
#-----------------------------------------------------------------------

Neval	= 1000                  # Evaluated test cases (for filenames)
N	= 1000                  # Number of test cases to plot

pfx     = "td"                  # Filename string prefix
sufx 	= ""                    # Filename string suffix
dir 	= "../res/"             # Results directory
Nes	= string(Neval)         # Turn Neval into a string
Ns      = string(N)             # Turn Ns into a string

tpfl    = Float64               # Floating point precision to use

hisfac  = Int(round(N/100))     # Histogram factor

#-----------------------------------------------------------------------
#       CONVERSION FACTORS
#-----------------------------------------------------------------------

mfact   = 0.004435              # conversion factor to meters
cmfact  = 0.4435                # conversion factor to centimeters
mmfact  = 4.435                 # conversion factor to millimeters

#-----------------------------------------------------------------------
ne	= 5
#-----------------------------------------------------------------------

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

δ1	= 0.001                 # Atmospheric perturbation parameter
δ2	= 0.10                  # Ionospheric perturbation parameter
sfx	= "n"*string(ne)*"p"*string(Int(round(δ2*100))) # Filename str.

tdLtup	= Serialization.deserialize(dir*pfx*"-"*Ns*"-"*sfx*sufx)
tdL     = squirrel.seval.tup2td( tdLtup )

lfact   = mfact                         # Choose conversion factors
unitstr = "m"                           # Unit string

thresh = (20,20,20)                             # Threshold tuples

h   = lfact*sort(tdL.erh,rev=true)              # horizontal errors
v   = lfact*sort(tdL.erv,rev=true)              # vertical errors
e   = lfact*sort(tdL.err,rev=true)              # total error

pf = dir*pfx*"-"*Ns*"-"*sfx*sufx                # Filename string

ers = PrintErr( h , v , e , thresh , pf , unitstr )     # Error string

binwidth = 0.05                                 # Histogram bin width
histL    = 0.0                                  # Histogram lower limit
histU    = 50.0                                 # Histogram upper limit

xrange   = (0,10)               # x range for histogram plot
yrange   = (0,15*hisfac)        # y range for histogram plot

histpar  = (binwidth,histL,histU)
labelsh  = ("Horizontal error ("*unitstr*")","Number of test cases")
labelsv  = ("Vertical error ("*unitstr*")","Number of test cases")

PlotErr( h , (ers[1],ers[4]) , histpar , xrange , yrange , labelsh ,
         dir*"plotHerr-"*sfx*".pdf" )
PlotErr( v , (ers[2],ers[5]) , histpar , xrange , yrange , labelsv ,
         dir*"plotVerr-"*sfx*".pdf" )

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

δ1	= 0.001                 # Atmospheric perturbation parameter
δ2	= 0.01                  # Ionospheric perturbation parameter
sfx	= "n"*string(ne)*"p"*string(Int(round(δ2*100))) # Filename str.

tdStup	= Serialization.deserialize(dir*pfx*"-"*Ns*"-"*sfx*sufx)
tdS     = squirrel.seval.tup2td( tdStup )

lfact   = cmfact                        # Choose conversion factors
unitstr = "cm"                          # Unit string

thresh = (20,20,20)                          # Threshold tuples

h   = lfact*sort(tdS.erh,rev=true)              # horizontal errors
v   = lfact*sort(tdS.erv,rev=true)              # vertical errors
e   = lfact*sort(tdS.err,rev=true)              # total error

pf = dir*pfx*"-"*Ns*"-"*sfx*sufx                # Filename string

ers = PrintErr( h , v , e , thresh , pf , unitstr )     # Error string

binwidth = 0.5                                  # Histogram bin width
histL    = 0.0                                  # Histogram lower limit
histU    = 500.0                                # Histogram upper limit

xrange   = (0,100)              # x range for histogram plot
yrange   = (0,15*hisfac)        # y range for histogram plot

histpar  = (binwidth,histL,histU)
labelsh  = ("Horizontal error ("*unitstr*")","Number of test cases")
labelsv  = ("Vertical error ("*unitstr*")","Number of test cases")

PlotErr( h , (ers[1],ers[4]) , histpar , xrange , yrange , labelsh ,
         dir*"plotHerr-"*sfx*".pdf" )
PlotErr( v , (ers[2],ers[5]) , histpar , xrange , yrange , labelsv ,
         dir*"plotVerr-"*sfx*".pdf" )

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

sfx	= "n"*string(ne)*"p0"

td0tup 	= Serialization.deserialize(dir*pfx*"-"*Ns*"-"*sfx*sufx)
td0     = squirrel.seval.tup2td( td0tup )

lfact   = mmfact                        # Choose conversion factors
unitstr = "mm"                          # Unit string

thresh = (20,20,20)                             # Threshold tuples

h   = lfact*sort(td0.erh,rev=true)              # horizontal errors
v   = lfact*sort(td0.erv,rev=true)              # vertical errors
e   = lfact*sort(td0.err,rev=true)              # total error

pf = dir*pfx*"-"*Ns*"-"*sfx*sufx                # Filename string

ers = PrintErr( h , v , e , thresh , pf , unitstr )     # Error string

binwidth = 0.01                                 # Histogram bin width
histL    = 0.0                                  # Histogram lower limit
histU    = 10.0                                 # Histogram upper limit

xrange   = (0,2)                # x range for histogram plot
yrange   = (0,20*hisfac)        # y range for histogram plot

histpar  = (binwidth,histL,histU)
labelsh  = ("Horizontal error ("*unitstr*")","Number of test cases")
labelsv  = ("Vertical error ("*unitstr*")","Number of test cases")

PlotErr( h , (ers[1],ers[4]) , histpar , xrange , yrange , labelsh ,
         dir*"plotHerr-"*sfx*".pdf" )
PlotErr( v , (ers[2],ers[5]) , histpar , xrange , yrange , labelsv ,
         dir*"plotVerr-"*sfx*".pdf" )

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

sfx	= "n"*string(ne)*"k"

tdktup  = Serialization.deserialize(dir*pfx*"-"*Ns*"-"*sfx*sufx)
tdk     = squirrel.seval.tup2td( tdktup )

lfact   = mmfact                        # Choose conversion factors
unitstr = "mm"                          # Unit string

thresh = (20,20,20)                             # Threshold tuples

h   = lfact*sort(tdk.erh,rev=true)              # horizontal errors
v   = lfact*sort(tdk.erv,rev=true)              # vertical errors
e   = lfact*sort(tdk.err,rev=true)              # total error

pf = dir*pfx*"-"*Ns*"-"*sfx*sufx                # Filename string

ers = PrintErr( h , v , e , thresh , pf , unitstr )     # Error string

binwidth = 0.0025                               # Histogram bin width
histL    = 0.0                                  # Histogram lower limit
histU    = 5.0                                # Histogram upper limit

xrange   = (0,0.5)              # x range for histogram plot
yrange   = (0,30*hisfac)        # y range for histogram plot

histpar  = (binwidth,histL,histU)
labelsh  = ("Horizontal error ("*unitstr*")","Number of test cases")
labelsv  = ("Vertical error ("*unitstr*")","Number of test cases")

PlotErr( h , (ers[1],ers[4]) , histpar , xrange , yrange , labelsh
        , dir*"plotHerr-"*sfx*".pdf" )
PlotErr( v , (ers[2],ers[5]) , histpar , xrange , yrange , labelsv
        , dir*"plotVerr-"*sfx*".pdf" )

pfC = pf*"C"

lfactC   = cmfact                       # Choose conversion factors
unitstrC = "cm"                         # Unit string

threshC = (2,5,20)                            # Threshold tuples

hC  = lfactC*sort(tdk.erhC,rev=true)            # horizontal errors
vC  = lfactC*sort(tdk.ervC,rev=true)            # vertical errors
eC  = lfactC*sort(tdk.errC,rev=true)            # total errors

ersC = PrintErr( hC , vC , eC , threshC , pfC , unitstrC )  # Error str

labelsh  = ("Horizontal error ("*unitstrC*")","Number of test cases")
labelsv  = ("Vertical error ("*unitstrC*")","Number of test cases")

binwidth = 0.00025                              # Histogram bin width
histpar  = (binwidth,histL,histU)
xrange   = (0,0.05)             # x range for histogram plot
yrange   = (0,20*hisfac)        # y range for histogram plot
PlotErr( hC , (ersC[1],ersC[4]) , histpar , xrange , yrange , labelsh ,
         dir*"plotHerr-"*sfx*"C.pdf" )

binwidth = 0.02                              # Histogram bin width
histpar  = (binwidth,histL,histU)
xrange   = (0,4)               # x range for histogram plot
yrange   = (0,4*hisfac)        # y range for histogram plot
PlotErr( vC , (ersC[2],ersC[5]) , histpar , xrange , yrange , labelsv ,
         dir*"plotVerr-"*sfx*"C.pdf" )

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


text:
@testset "thin_edges" begin
    @info "Test: thin_edges"

    @testset "API" begin
        img_ref = Gray{N0f8}.(testimage("lake_gray.tif"))
        img = copy(img_ref)

        f = NonmaximaSuppression()
        g₁, g₂ = imgradients(img, KernelFactors.scharr)
        mag = hypot.(g₁, g₂)

        g₁_cpy = copy(g₁)
        g₂_cpy = copy(g₂)
        mag_cpy = copy(mag)

        thinned_edges₀ = zeros(eltype(mag), axes(mag))
        thin_edges!(thinned_edges₀, mag, g₁, g₂, f)
        thinned_edges₁ = thin_edges(mag, g₁, g₂, f)

        @test mag == mag_cpy # mag unchanged
        @test g₁ == g₁_cpy # g₁ unchanged
        @test g₂ == g₂_cpy # g₂ unchanged

        @test thinned_edges₀ == thinned_edges₁
        @test eltype(thinned_edges₀) == eltype(mag)
        @test eltype(thinned_edges₁) == eltype(mag)

        thin_edges!(mag, g₁, g₂, f)
        @test mag != mag_cpy
        @test thinned_edges₀ == mag

        mag = copy(mag_cpy)
        for T in generate_test_types([Float16, Float32, Float64], [Gray])
            thinned_edges₂ = zeros(T, axes(mag))
            thin_edges!(thinned_edges₂, mag, g₁, g₂, f)
            @test eltype(thinned_edges₂) == T
            @test eltype(thin_edges(T, mag, g₁, g₂, f)) == T
        end
    end

    @testset "Numerical" begin
        img = Gray{N0f8}.(load("algorithms/References/circle.png"))
        nms_ref = Gray{N0f8}.(load("algorithms/References/circle_nms.png"))

        f = NonmaximaSuppression()

        g₁, g₂ = imgradients(img, KernelFactors.scharr)
        mag = hypot.(g₁, g₂)

        thinned_edges₀ = zeros(eltype(mag), axes(mag))
        thin_edges!(thinned_edges₀, mag, g₁, g₂, f)

        thinned_edges₁ = thin_edges(mag, g₁, g₂, f)

        for i in CartesianIndices(nms_ref)
            nms_ref[i] == 0.0 ? (@test thinned_edges₀[i] == 0) : (@test thinned_edges₀[i] == mag[i])
            nms_ref[i] == 0.0 ? (@test thinned_edges₁[i] == 0) : (@test thinned_edges₁[i] == mag[i])
        end
    end
end


text:
# # Partitions
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/sketches/Partitions.ipynb)
#
# Partitions are a categorical construction that we derive from sets and functions.
# Given a set A, you can think of all of the ways to partition A into parts.
# These ways of partitioning are isomorphic to equivalence relations R ⊆ A × A.

# The first step is our Catlab imports
using Core: GeneratedFunctionStub
using Test

using Catlab, Catlab.Theories, Catlab.CategoricalAlgebra, Catlab.CategoricalAlgebra.FinSets
import Catlab.Theories: compose
using DataStructures
using PrettyTables
PrettyTables.pretty_table(f::FinFunction, name::Symbol=:f) = pretty_table(OrderedDict(:x=>1:dom(f).set, Symbol("$(name)(x)")=>collect(f)))
using LaTeXStrings

Quiversty = L"""
% contents of quiver.sty
% `tikz-cd` is necessary to draw commutative diagrams.
\RequirePackage{tikz-cd}
% `amssymb` is necessary for `\lrcorner` and `\ulcorner`.
\RequirePackage{amssymb}
% `calc` is necessary to draw curved arrows.
\usetikzlibrary{calc}
% `pathmorphing` is necessary to draw squiggly arrows.
\usetikzlibrary{decorations.pathmorphing}

% A TikZ style for curved arrows of a fixed height, due to AndréC.
\tikzset{curve/.style={settings={#1},to path={(\tikztostart)
    .. controls ($(\tikztostart)!\pv{pos}!(\tikztotarget)!\pv{height}!270:(\tikztotarget)$)
    and ($(\tikztostart)!1-\pv{pos}!(\tikztotarget)!\pv{height}!270:(\tikztotarget)$)
    .. (\tikztotarget)\tikztonodes}},
    settings/.code={\tikzset{quiver/.cd,#1}
        \def\pv##1{\pgfkeysvalueof{/tikz/quiver/##1}}},
    quiver/.cd,pos/.initial=0.35,height/.initial=0}

% TikZ arrowhead/tail styles.
\tikzset{tail reversed/.code={\pgfsetarrowsstart{tikzcd to}}}
\tikzset{2tail/.code={\pgfsetarrowsstart{Implies[reversed]}}}
\tikzset{2tail reversed/.code={\pgfsetarrowsstart{Implies}}}
% TikZ arrow styles.
\tikzset{no body/.style={/tikz/dash pattern=on 0 off 1mm}}
""";

# ## FinSet: the category of Finite Sets
# In FinSet the objects are sets n = {1...n} and the morphisms are functions between finite sets.
# You can wrap a plain old Int into a finite set with the FinSet(n::Int) function. These sets will
# serve as the domain and codomains of our functions.
n = FinSet(3)
m = FinSet(4)

# once you have some sets, you can define functions between them. A FinFunction from n to m, f:n→m, can be specified as
# an array of length n with elements from m.

f = FinFunction([2,4,3], n, m)

pretty_table(f)

# ## Surjective maps
# In order to use a map to represent a partition, we have to make sure that it is surjective.
# Given a FinFunction, we can compute the preimage of any element in its codomain.

preimage(f, 2)

preimage(f, 1)

# If the preimage is empty, then there is no element in the domain that maps to that element of the codomain.
# This gives us a definition of surjective functions by asserting that all the preimages are nonempty.
# Julia note: !p is the predicate x ↦ ¬p(x), f.(A) applies f to all of the elements in A.

is_surjective(f::FinFunction) = all((!isempty).(preimage(f,i) for i in codom(f)))
is_surjective(f)

# Our function f, wasn't surjective so it can't be used to induce a partition via its preimages.
# Let's try again,

g = FinFunction([1,2,3,3], m, n)
pretty_table(g, :g)
is_surjective(g)

# # Refinements of a Partition
# When defining partitions classically as A = ∪ₚ Aₚ with p ≠ r ⟹ Aₚ ≠ Aᵣ,
# it is not immediately obvious how to define comparisons between partitions.
# With the "a partition of A is a surjective map out of A" definition, the comparisons
# are obvious. The composition of surjective maps is surjective, so we can define
# the refinement order in terms of a diagram in Set.
#
# You can see a graphical definition in [quiver](https://q.uiver.app/?q=WzAsMyxbMCwwLCJBIl0sWzMsMCwiUSJdLFszLDIsIlAiXSxbMSwyLCJoIiwwLHsic3R5bGUiOnsiaGVhZCI6eyJuYW1lIjoiZXBpIn19fV0sWzAsMSwiZiIsMCx7InN0eWxlIjp7ImhlYWQiOnsibmFtZSI6ImVwaSJ9fX1dLFswLDIsImciLDIseyJzdHlsZSI6eyJoZWFkIjp7Im5hbWUiOiJlcGkifX19XV0=)

using TikzCDs

triangle = L"""
A &&& Q \\
\\
&&& P
\arrow["h", two heads, from=1-4, to=3-4]
\arrow["f", two heads, from=1-1, to=1-4]
\arrow["g"', two heads, from=1-1, to=3-4]
""";

TikzCD(triangle, preamble=Quiversty)

# Let's take a look at an example:

A = FinSet(4)
Q = FinSet(3)
P = FinSet(2)

f = FinFunction([1,2,3,3], A, Q)
g = FinFunction([1,1,2,2], A, P)
h = FinFunction([1,1,2], Q, P)

@test_throws ErrorException compose(g,h) #Catlab checks the domains match

pretty_table(compose(f,h), Symbol("(f⋅h)"))

compose(f,h) == g

# This triangle commutes, so f is a refinement of g equivalently g is coarser than f.

h′ = FinFunction([1,1], P, FinSet(1))

pretty_table(f⋅h⋅h′, Symbol("f⋅h⋅h′"))

# ### Properties of refinements
# We can show that refinement gives us a preorder on partitions directly from the
# nice properties of surjective maps.
# 1. Reflexive: Any partition is a refinement of itself.
# 2. Transitive: If f ≤ g ≤ h as partitions, then f ≤ h
# You can read these directly off the definition of refinements as a commutative
# triangle in the category of (Set, Surjections).
# You can edit this diagram in [quiver](https://q.uiver.app/?q=WzAsNCxbMCwwLCJBIl0sWzMsMCwiUSJdLFszLDIsIlAiXSxbMyw0LCJRXlxccHJpbWUiXSxbMSwyLCJoIl0sWzAsMSwiZiIsMCx7InN0eWxlIjp7ImhlYWQiOnsibmFtZSI6ImVwaSJ9fX1dLFswLDIsImciLDIseyJzdHlsZSI6eyJoZWFkIjp7Im5hbWUiOiJlcGkifX19XSxbMiwzLCJoXlxccHJpbWUiXSxbMCwzLCJmXFxjZG90IGhcXGNkb3QgaF5cXHByaW1lID0gZ1xcY2RvdCBoXlxccHJpbWUiLDIseyJzdHlsZSI6eyJoZWFkIjp7Im5hbWUiOiJlcGkifX19XV0=)
refinement = L"""
A &&& Q \\
\\
&&& P \\
\\
&&& {Q^\prime}
\arrow["h", from=1-4, to=3-4]
\arrow["f", two heads, from=1-1, to=1-4]
\arrow["g"', two heads, from=1-1, to=3-4]
\arrow["{h^\prime}", from=3-4, to=5-4]
\arrow["{f\cdot h\cdot h^\prime = g\cdot h^\prime}"', two heads, from=1-1, to=5-4]
""";

TikzCD(refinement, preamble=Quiversty)


text:
# stemplots

#JL Plots
Using Plots

x = sin.(0:0.3:2pi)
plot(x, line=:stem, marker=:star, markersize=20)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
#=
Filtering
=#

function _make_filter_gui()

    #Bandwidth Adjustment

    band_hw_frame=Frame("Hardware Filters")
    band_grid=Grid()
    band_hw_grid=Grid()
    band_grid[1,1]=band_hw_frame
    push!(band_hw_frame,band_hw_grid)

    band_sb1=SpinButton(0:1000)
    Gtk.GAccessor.value(band_sb1,300)
    band_hw_grid[1,1]=band_sb1
    band_hw_grid[2,1]=Label("Lower Bandwidth")

    band_sb2=SpinButton(1000:10000)
    Gtk.GAccessor.value(band_sb2,5000)
    band_hw_grid[1,2]=band_sb2
    band_hw_grid[2,2]=Label("Higher BandWidth")

    band_sb3=SpinButton(0:1000)
    Gtk.GAccessor.value(band_sb3,300)
    band_hw_grid[1,3]=band_sb3
    band_hw_grid[2,3]=Label("DSP High Pass")

    band_b1=Button("Update")
    band_hw_grid[1,4]=band_b1

    band_sw_frame=Frame("Software Filters")
    band_sw_grid=Grid()
    band_grid[1,2]=band_sw_frame
    push!(band_sw_frame,band_sw_grid)

    band_sw_sb3=SpinButton(1:1)
    Gtk.GAccessor.value(band_sw_sb3,1)
    band_sw_grid[1,1]=band_sw_sb3
    band_sw_grid[2,1]=Label("Channel Number")

    band_sw_sb4=SpinButton(1:1)
    Gtk.GAccessor.value(band_sw_sb4,1)
    band_sw_grid[1,2]=band_sw_sb4
    band_sw_grid[2,2]=Label("Filter Number")

    band_sw_check=CheckButton()
    band_sw_grid[1,3]=band_sw_check
    band_sw_grid[2,3]=Label("All Channels")

    filter_combo = ComboBoxText()
    for choice in ["High Pass"; "Low Pass"; "BandPass"; "BandStop"]
        push!(filter_combo,choice)
    end
    setproperty!(filter_combo,:active,0)

    band_sw_grid[1,4]=filter_combo
    band_sw_grid[2,4]=Label("Filter Type")

    band_sw_sb1=SpinButton(0:10000)
    Gtk.GAccessor.value(band_sw_sb1,10)
    band_sw_grid[1,5]=band_sw_sb1
    band_sw_sb1_l=Label("High Pass Cutoff")
    band_sw_grid[2,5]=band_sw_sb1_l

    band_sw_sb2=SpinButton(0:10000)
    Gtk.GAccessor.value(band_sw_sb2,10)
    band_sw_grid[1,6]=band_sw_sb2
    band_sw_sb2_l=Label("")
    band_sw_grid[2,6]=band_sw_sb2_l

    filter_combo_output = ComboBoxText()
    for choice in ["Spikes"; "LFP"]
        push!(filter_combo_output,choice)
    end
    setproperty!(filter_combo_output,:active,0)
    band_sw_grid[1,7]=filter_combo_output
    band_sw_grid[2,7]=Label("Output of Filter")

    band_sw_b1=Button("Add New")
    band_sw_grid[1,8]=band_sw_b1

    band_sw_b2=Button("Replace")
    band_sw_grid[2,8]=band_sw_b2

    band_sw_b_delete=Button("Delete")
    band_sw_grid[1,9]=band_sw_b_delete

    band_sw_c=Canvas(200,300)
    band_sw_grid[2,10]=band_sw_c

    filt_list = ListStore(String,Int32,Int32)

    filt_tv = TreeView(TreeModel(filt_list))
    filt_rtext2=CellRendererText()
    filt_rtext3=CellRendererText()
    filt_rtext4=CellRendererText()

    filt_c2 = TreeViewColumn("Filter Type",filt_rtext2,Dict([("text",1)]))
    filt_c3 = TreeViewColumn("Wn1",filt_rtext3,Dict([("text",2)]))
    filt_c4 = TreeViewColumn("Wn2",filt_rtext4,Dict([("text",3)]))

    filt_tv_s = Gtk.GAccessor.selection(filt_tv)
    #Gtk.GAccessor.mode(filt_tv_s,Gtk.GConstants.GtkSelectionMode.MULTIPLE)

    push!(filt_tv,filt_c2)
    push!(filt_tv,filt_c3)
    push!(filt_tv,filt_c4)

    filt_scroll=ScrolledWindow()
    Gtk.GAccessor.min_content_height(filt_scroll,500)
    Gtk.GAccessor.min_content_width(filt_scroll,250)
    push!(filt_scroll,filt_tv)

    band_grid[2,2]=filt_scroll

    band_win=Window(band_grid)
    setproperty!(band_win, :title, "Filtering")

    Gtk.showall(band_win)
    visible(band_win,false)
    Gtk.visible(band_sw_sb2,false)

    band_widgets=Band_Widgets(band_win,band_sb1,band_sb2,band_sb3,band_b1,filter_combo,band_sw_sb1,band_sw_sb2,band_sw_sb3,band_sw_b1,band_sw_b2,band_sw_b_delete,band_sw_check,band_sw_sb1_l,band_sw_sb2_l,filter_combo_output,band_sw_sb4,band_sw_c,10,10,1,1,0,1,falses(0),filt_tv,filt_list)
end

function add_filter_callbacks(w,handles,r,fpga)

    #Hardware Filters
    signal_connect(band_b1_cb,w.b1,"clicked",Void,(),false,(handles,fpga))

    signal_connect(w.win, :delete_event) do widget, event
        visible(w.win, false)
        true
    end

    signal_connect(filter_type_cb,w.sw_box, "changed",Void,(),false,(handles,r))
    signal_connect(change_channel_cb,w.sw_chan_sb,"value-changed",Void,(),false,(handles,r))

    #Software Filter bandwidth
    signal_connect(change_wn1_cb,w.wn_sb1,"value-changed",Void,(),false,(handles,r))
    signal_connect(change_wn2_cb,w.wn_sb2,"value-changed",Void,(),false,(handles,r))

    #Add Remove software Filter
    signal_connect(add_filter_cb,w.add_b,"clicked",Void,(),false,(handles,r))
    signal_connect(replace_filter_cb,w.re_b,"clicked",Void,(),false,(handles,r))
    signal_connect(delete_filter_cb,w.delete_b,"clicked",Void,(),false,(handles,r))

    signal_connect(change_filt_output_cb,w.output_box,"changed",Void,(),false,(handles,r))
    signal_connect(change_pos_cb,w.filt_num_sb,"value-changed",Void,(),false,(r,handles))

end

function make_filter(rhd::RHD2000,filter_type::Int64,wn1,wn2=0.0)

    if filter_type == 1
        responsetype = Highpass(wn1; fs = rhd.sr)
    elseif filter_type == 2
        responsetype = Lowpass(wn1; fs = rhd.sr)
    elseif filter_type == 3
        responsetype = Bandpass(wn1,wn2; fs=rhd.sr)
    elseif filter_type == 4
        responsetype = Bandstop(wn1,wn2; fs=rhd.sr)
    end
    designmethod=Butterworth(4)
    df1=digitalfilter(responsetype,designmethod)
    DF2TFilter(df1)
end

function apply_filter(rhd::RHD2000,ff::Intan_Filter,chan_num::Int64)

    temp=zeros(Float64,SAMPLES_PER_DATA_BLOCK)
    for i=1:SAMPLES_PER_DATA_BLOCK
        temp[i]=convert(Float64,rhd.v[i,chan_num])
    end

    filt!(temp,ff.filt,temp)

    if ff.output==0
        for i=1:SAMPLES_PER_DATA_BLOCK
            rhd.v[i,chan_num] = round(Int16,temp[i])
        end
    else
        for i=1:SAMPLES_PER_DATA_BLOCK
            rhd.lfps[i,chan_num] = round(Int16,temp[i])
        end
    end
    nothing
end

#=
Callbacks
=#

function filter_type_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    filter_type = getproperty(han.band_widgets.sw_box,:active,Int64)+1
    han.band_widgets.f_type = filter_type

    if filter_type == 1
        Gtk.visible(han.band_widgets.wn_sb1,true)
        Gtk.visible(han.band_widgets.wn_sb2,false)
        setproperty!(han.band_widgets.wn_sb1_l,:label,"High Pass Cutoff")
        setproperty!(han.band_widgets.wn_sb2_l,:label,"")
    elseif filter_type == 2
        Gtk.visible(han.band_widgets.wn_sb1,true)
        Gtk.visible(han.band_widgets.wn_sb2,false)
        setproperty!(han.band_widgets.wn_sb1_l,:label,"Low Pass Cutoff")
        setproperty!(han.band_widgets.wn_sb2_l,:label,"")
    elseif filter_type == 3
        Gtk.visible(han.band_widgets.wn_sb1,true)
        Gtk.visible(han.band_widgets.wn_sb2,true)
        setproperty!(han.band_widgets.wn_sb1_l,:label,"Low Cutoff")
        setproperty!(han.band_widgets.wn_sb2_l,:label,"High Cutoff")
    elseif filter_type == 4
        Gtk.visible(han.band_widgets.wn_sb1,true)
        Gtk.visible(han.band_widgets.wn_sb2,true)
        setproperty!(han.band_widgets.wn_sb1_l,:label,"Low Cutoff")
        setproperty!(han.band_widgets.wn_sb2_l,:label,"High Cutoff")
    end

    nothing
end

function add_filter_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    #Create Filter
    han,rhd = user_data

    filt_type = han.band_widgets.f_type

    chan_num = han.band_widgets.chan

    wn1 = han.band_widgets.wn1
    wn2 = han.band_widgets.wn2

    pos = han.band_widgets.f_pos

    output = han.band_widgets.f_out

    if (getproperty(han.band_widgets.sw_check,:active,Bool))
        for i=1:size(rhd.v,2)
            myfilt=make_filter(rhd,filt_type,wn1,wn2)
            if pos > length(rhd.filts[i])
                push!(rhd.filts[i],Intan_Filter(i,output,wn1,wn2,filt_type,myfilt))
            else
                insert!(rhd.filts[i],pos,Intan_Filter(i,output,wn1,wn2,filt_type,myfilt))
            end

            if output==1
                han.band_widgets.lfp_en[i]=true
            end
        end
    else
        #add new filter
        myfilt=make_filter(rhd,filt_type,wn1,wn2)
        if pos>length(rhd.filts[chan_num])
            push!(rhd.filts[chan_num],Intan_Filter(chan_num,output,wn1,wn2,filt_type,myfilt))
        else
            insert!(rhd.filts[chan_num],pos,Intan_Filter(chan_num,output,wn1,wn2,filt_type,myfilt))
        end
        if output==1
            han.band_widgets.lfp_en[chan_num]=true
        end
    end

    draw_filter_canvas(han,rhd,chan_num)
    Gtk.GAccessor.range(han.band_widgets.filt_num_sb,1,length(rhd.filts[chan_num])+1)

    nothing
end

function replace_filter_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    filt_type = han.band_widgets.f_type

    chan_num = han.band_widgets.chan

    wn1 = han.band_widgets.wn1
    wn2 = han.band_widgets.wn2

    pos = han.band_widgets.f_pos

    output = han.band_widgets.f_out

    myfilt=make_filter(rhd,filt_type,wn1,wn2)

    if rhd.filts[chan_num][pos].f_out == 1
        if output == 0
            han.band_widgets.lfp_en[chan_num]=false
        end
    end
    if output == 1
        han.band_widgets.lfp_en[chan_num]=true
    end

    rhd.filts[chan_num][pos]=Intan_Filter(rhd.filts[chan_num][pos].chan,output,wn1,wn2,filt_type,myfilt)

    draw_filter_canvas(han,rhd,pos)

    nothing
end

function band_adj_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})
    han, rhd = user_data

    visible(han.band_widgets.win,true)
end

function band_b1_cb(widget::Ptr,user_data::Tuple{Gui_Handles,Array{I,1}}) where I<:IC

    han, myic = user_data

    lower=getproperty(han.band_widgets.sb1,:value,Int64)
    upper=getproperty(han.band_widgets.sb2,:value,Int64)
    dsp_lower=getproperty(han.band_widgets.sb3,:value,Int64)

    change_bandwidth(myic,lower,upper,dsp_lower)

    nothing
end

function change_wn1_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    han.band_widgets.wn1 = getproperty(han.band_widgets.wn_sb1,:value,Int64)

    draw_filter_canvas(han,rhd,han.band_widgets.f_pos)

    nothing
end

function change_wn2_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    han.band_widgets.wn2 = getproperty(han.band_widgets.wn_sb2,:value,Int64)

    draw_filter_canvas(han,rhd,han.band_widgets.f_pos)

    nothing
end

function change_pos_cb(widget::Ptr,user_data::Tuple{RHD2000,Gui_Handles})

    rhd, han = user_data

    han.band_widgets.f_pos=getproperty(han.band_widgets.filt_num_sb,:value,Int64)

    draw_filter_canvas(han,rhd,han.band_widgets.f_pos)

    nothing
end

function change_filt_output_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    han.band_widgets.f_out = getproperty(han.band_widgets.output_box,:active,Int64)

    draw_filter_canvas(han,rhd,han.band_widgets.f_pos)

    nothing
end

function delete_filter_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    if han.band_widgets.f_pos > length(rhd.filts[han.band_widgets.chan])
    else

        if rhd.filts[han.band_widgets.chan][han.band_widgets.f_pos].output == 1
            han.band_widgets.lfp_en[han.band_widgets.chan] = false
        end

        deleteat!(rhd.filts[han.band_widgets.chan],han.band_widgets.f_pos)

        setproperty!(han.band_widgets.filt_num_sb,:value,1)

        Gtk.GAccessor.range(han.band_widgets.filt_num_sb,1,length(rhd.filts[han.band_widgets.chan])+1)

        draw_filter_canvas(han,rhd,han.band_widgets.f_pos)
    end

    nothing
end

function change_channel_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    han.band_widgets.chan=getproperty(han.band_widgets.sw_chan_sb,:value,Int64)

    Gtk.GAccessor.range(han.band_widgets.filt_num_sb,1,length(rhd.filts[han.band_widgets.chan])+1)

    setproperty!(han.band_widgets.filt_num_sb,:value,1)

    draw_filter_canvas(han,rhd,1)

    nothing
end

function draw_filter_canvas(han::Gui_Handles,rhd,pos)

    ctx = Gtk.getgc(han.band_widgets.c)

    set_source_rgb(ctx,1.0,1.0,1.0)
    paint(ctx)
    set_source_rgb(ctx,0.0,0.0,0.0)

    num=han.band_widgets.chan

    draw_hardware_filter(han,num,ctx)

    split=false
    lr=1

    for i=1:length(rhd.filts[num])

        if rhd.filts[num][i].output==1
            lr=0
            split=true
        elseif split
            lr=2
        else
            lr=1
        end

        if i != pos
            draw_software_filter(rhd.filts[num][i],ctx,i,lr,(0.0,0.0,0.0),false)
        else
            draw_software_filter(Intan_Filter(num,han.band_widgets.f_out,han.band_widgets.wn1,han.band_widgets.wn2,han.band_widgets.f_type,make_filter(rhd,1,10,10)),ctx,i,lr,(1.0,0.0,0.0),false)
        end
    end

    if pos>length(rhd.filts[num])
        if han.band_widgets.f_out==1
            lr=0
            split=true
        elseif split
            lr=2
        else
            lr=1
        end
        draw_software_filter(Intan_Filter(num,han.band_widgets.f_out,han.band_widgets.wn1,han.band_widgets.wn2,han.band_widgets.f_type,make_filter(rhd,1,10,10)),ctx,pos,lr,(1.0,0.0,0.0),true)
    end

    reveal(han.band_widgets.c)
    nothing
end

function draw_hardware_filter(han,num,ctx)

    move_to(ctx,30,15)
    show_text(ctx,"RHD Hardware")
    move_to(ctx,30,25)
    show_text(ctx,string("Analog Bandpass: ",300," - ",3000))
    move_to(ctx,30,35)
    show_text(ctx,string("Hardware Highpass DSP: ",300))

    draw_box(25,5,175,45,(0.0,0.0,0.0),1.0,ctx)

    nothing
end

function draw_software_filter(myfilt,ctx,i,lr,mycolor,myselect)

    yinit=85+(i-1)*80

    if lr==0 #LFP left
        x1=25
        x2=100
        line(ctx,x1,x1,yinit-40,yinit+5)
        line(ctx,175,175,yinit-40,yinit+85)
        line(ctx,x1,50,yinit-40,yinit-40)
        line(ctx,150,175,yinit-40,yinit-40)
    elseif lr==1
        x1=50
        x2=150
        line(ctx,100,100,yinit-40,yinit+5)
    else
        x1=100
        x2=175
        line(ctx,x2,x2,yinit-40,yinit+5)
    end

    if myselect
        dashes=[5.0,5.0,5.0]
        set_dash(ctx,dashes,0.0)
    end

    draw_box(x1,yinit+5,x2,yinit+45,mycolor,1.0,ctx)

    if myselect
        set_dash(ctx,Float64[])
    end

    wn1 = myfilt.wn1
    wn2 = myfilt.wn2

    move_to(ctx,x1+5,yinit+15)
    if (myfilt.f_type==1)|(myfilt.f_type==2)

        if myfilt.f_type==1
            show_text(ctx,"High Pass")
        else
            show_text(ctx,"Low Pass")
        end
        move_to(ctx,x1+5,yinit+25)
        show_text(ctx,string("Cutoff: ", wn1))
    else

        if myfilt.f_type==3
            show_text(ctx,"BandPass")
        else
            show_text(ctx,"BandStop")
        end
        move_to(ctx,x1+5,yinit+25)
        show_text(ctx,string("Band: ", wn1, " - ", wn2))
    end

    nothing
end


text:
using .Networks

using Printf
using DataStructures

mutable struct UncertainValue{T}
    definitive_value::Union{T, Missing}
    anticipated_value::SortedDict{String, Union{T, Missing}}
end

function UncertainValue{T}() where T
    return UncertainValue{T}( missing, SortedDict{String, T}() )
end
function UncertainValue{T}(scenarios::Vector{String}) where T
    antcipated_values = SortedDict(zip(scenarios, fill(missing, length(scenarios))))
    return UncertainValue{T}(missing, antcipated_values)
end

function get_scenarios(uncertain_value::UncertainValue{T}) where T
    return keys(uncertain_value.anticipated_value)
end

function is_missing_values(uncertain_value::UncertainValue{T})::Bool where T
    for (_,value) in uncertain_value.anticipated_value
        if ismissing(value)
            return true
        end
    end
    return false
end

function make_non_definitive(uncertain_value::UncertainValue{T}) where T
    uncertain_value.definitive_value = missing
end

function set_value!(uncertain_value::UncertainValue{T}, scenario::String, value::T)::Union{T, Missing} where T
    if is_definitive(uncertain_value)
        make_non_definitive(uncertain_value)
    end
    uncertain_value.anticipated_value[scenario] = value
end

function is_definitive(uncertain_value::UncertainValue{T})::Bool where T
    existing_value = get_value(uncertain_value)
    return !ismissing(existing_value)
end

function safeset_definitive_value!(uncertain_value::UncertainValue{T}, value::T)::Union{T, Missing} where T
    if is_definitive(uncertain_value)
        existing_value = get_value(uncertain_value)
        if is_different(existing_value, value)
            msg = @sprintf("Unable to set definitive value to %s : The definitive value `%s` was already set to the UncertainValue.",
                            value, existing_value)
            throw( error(msg) )
        end
    else
        uncertain_value.definitive_value = value
        for (scenario, _) in uncertain_value.anticipated_value
            uncertain_value.anticipated_value[scenario] = value
        end
    end
    return uncertain_value.definitive_value
end
function set_definitive_value!(uncertain_value::UncertainValue{T}, value::T)::Union{T, Missing} where T
    if is_definitive(uncertain_value)
        existing_value = get_value(uncertain_value)
        if existing_value != value
            msg = @sprintf("changed value to %s : The definitive value was already set to `%s` for the UncertainValue.",
                    value, existing_value)
            @debug msg
        end
    end

    uncertain_value.definitive_value = value
    for (scenario, _) in uncertain_value.anticipated_value
        uncertain_value.anticipated_value[scenario] = value
    end
    return uncertain_value.definitive_value
end

function get_value(uncertain_value::UncertainValue{T})::Union{T, Missing} where T
    return uncertain_value.definitive_value
end

function safeget_value(uncertain_value::UncertainValue{T})::T where T
    if is_definitive(uncertain_value)
        return get_value(uncertain_value)
    else
        throw( error("UncertainValue is not definitive") )
    end
end

function get_value(uncertain_value::UncertainValue{T}, scenario::String)::Union{Missing,T} where T
    if is_definitive(uncertain_value)
        return get_value(uncertain_value)
    else
        if haskey(uncertain_value.anticipated_value, scenario)
            return uncertain_value.anticipated_value[scenario]
        else
            return missing
        end
    end
end

function safeget_value(uncertain_value::UncertainValue{T}, scenario::String)::T where T
    value = get_value(uncertain_value, scenario)
    if !ismissing( value )
        return value
    else
        msg = @sprintf("No decision was made for scenario %s yet", scenario)
        throw( error(msg) )
    end
end

function Base.show(io::IO, uncertain_value::UncertainValue{T}) where T
    scenarios = get_scenarios(uncertain_value)
    if is_definitive(uncertain_value)
        @printf(io, "definitive value %s for scenarios %s",
                get_value(uncertain_value), scenarios)
    else
        for s in scenarios
            @printf(io, "%s:%s,", s, get_value(uncertain_value, s))
        end
    end
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
#
# Data
#

survo_puzzle_problems  = Dict(

# http://en.wikipedia.org/wiki/Survo_Puzzle, first example
#
# Solution:
#  12 6 2 10
#  8 1 5 4
#  7 9 3 11
#
:1 => Dict(
        :row_sums => [30,18,30],
        :col_sums => [27,16,10,25],
        :problem => [[0, 6, 0, 0],
                   [8, 0, 0, 0],
                   [0, 0, 3, 0]]
      ),

# http://en.wikipedia.org/wiki/Survo_Puzzle, second example
# difficulty 0
:2  => Dict(
        :row_sums => [9, 12],       # :row_sums
        :col_sums => [9, 7, 5],     # :col_sums
        :problem => [[0, 0, 3],  # :problem
                   [0, 6, 0]]
       ),

# http://en.wikipedia.org/wiki/Survo_Puzzle, third example
# difficulty 150 ("open puzzle", i.e. no hints]
# It's an unique solution.
# (817 propagations with Gecode/fz, and 33 failures, 88 commits]
# r = 3;
# c = 4;
# :row_sums = [24,15,39];
# :col_sums = [21,10,18,29];
# matrix = array2d(1..r, 1..c,
#   [
#     0, 0, 0, 0,
#     0, 0, 0, 0,
#     0, 0, 0, 0
#   ]];
# Note: this version has no hints
:3  => Dict(
        :row_sums => [24,15,39],      # :row_sums
        :col_sums => [21,10,18,29],   # :col_sums
        :problem => [[0, 0, 0, 0], # :problem
                   [0, 0, 0, 0],
                   [0, 0, 0, 0]]
        ),

# same as above but with hints: difficulty 0
# (15 propagations with Gecode/fz, no failures, no commits]
# matrix = array2d(1..r, 1..c,
#    [
#      7, 0, 5, 0,
#      0, 1, 0, 8,
#      0, 0, 11, 0
#    ]];
:4  => Dict(
       :row_sums => [24,15,39],      # :row_sums
       :col_sums => [21,10,18,29],   # :col_sums
       :problem => [[7, 0, 5, 0], # :problem
                  [0, 1, 0, 8],
                  [0, 0, 11, 0]]
       ),

# http://www.survo.fi/puzzles/280708.txt, third puzzle
# Survo puzzle 128/2008 (1700] #364-35846
#
#    A  B  C  D  E  F
# 1  *  *  *  *  *  * 30
# 2  *  * 18  *  *  * 86
# 3  *  *  *  *  *  * 55
#   22 11 42 32 27 37
:5  => Dict(
       :row_sums => [30, 86, 55],
       :col_sums => [22, 11, 42, 32, 27, 37],
       :problem => [[0, 0,  0, 0, 0, 0],
                  [0, 0, 18, 0, 0, 0],
                  [0, 0,  0, 0, 0, 0]]
        ),
#
# http://en.wikipedia.org/wiki/Survo_Puzzle, under "Swapping method"
# (open puzzle]
#
:6  => Dict(
       :row_sums => [51,36,32,17],
       :col_sums => [51,42,26,17],
       :problem => [[0, 0, 0, 0],
                  [0, 0, 0, 0],
                  [0, 0, 0, 0],
                  [0, 0, 0, 0]]
      ),
)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
export HermiteSMF, FixedHermiteSMF

struct HermiteSMF<:SeqFilter
		"Filter function"
		G::Function

		"Standard deviations of the measurement noise distribution"
		ϵy::AdditiveInflation

        "HermiteMap"
        M::HermiteMap

        "Observation dimension"
        Ny::Int64

        "State dimension"
        Nx::Int64

        "Time step dynamic"
        Δtdyn::Float64

        "Time step observation"
        Δtobs::Float64

		"Time step between reparametrization of the map"
		Δtfresh::Float64

		"Boolean: is state vector filtered"
		isfiltered::Bool
end

function Base.show(io::IO, smf::HermiteSMF)
	println(io,"Hermite Stochastic Map Filter with filtered = $(smf.isfiltered)")
end

function (smf::HermiteSMF)(X, ystar::Array{Float64,1}, t::Float64)
	Ny = smf.Ny
	Nx = smf.Nx

	Nystar = size(ystar, 1)
	Nypx, Ne = size(X)

	@assert Nystar == Ny "Size of ystar is not consistent with Ny"
	@assert Nypx == Ny + Nx "Size of X is not consistent with Ny or Nx"

	# Perturbation of the measurements
	smf.ϵy(X, 1, Ny)

	# if abs(round(Int64,  t / smf.Δtfresh) - t / smf.Δtfresh)<1e-6
	M = HermiteMap(30, X; diag = true, b = "CstLinProHermite")
		# Perform a kfold optimization of the map
	optimize(M, X, "kfolds"; withconstant = false, withqr = true,
		     verbose = false, start = Ny+1, P = serial, hessprecond = true)
	# else
	# 	L = LinearTransform(X; diag = true)
	# 	M = HermiteMap(40, Ny+Nx, L, smf.M.C)
	# 	# M = HermiteMap(40, X; diag = true)
	# 	# Only optimize the existing coefficients of the basis
	# 	optimize(M, X, nothing; withconstant = false, withqr = true,
	# 		   verbose = false, start = Ny+1, P = serial, hessprecond = true)
	# end

	# Evaluate the transport map
	F = evaluate(M, X; apply_rescaling = true, start = Ny+1, P = serial)

	# # Rescale ystar
	# ystar .-= view(M.L.μ,1:Ny)
	# ystar ./= M.L.L.diag[1:Ny]

	# Generate the posterior samples by partial inversion of the map
	hybridinverse!(X, F, M, ystar; start = Ny+1, P = serial)
	# @show getcoeff(M[Nypx])
	# @show "after inversion"
	# @show norm(X)
	return X
end

struct FixedHermiteSMF<:SeqFilter
		"Filter function"
		G::Function

		"Standard deviations of the measurement noise distribution"
		ϵy::AdditiveInflation

        "HermiteMap"
        M::HermiteMap

        "Observation dimension"
        Ny::Int64

        "State dimension"
        Nx::Int64

        "Time step dynamic"
        Δtdyn::Float64

        "Time step observation"
        Δtobs::Float64

		"Boolean: is state vector filtered"
		isfiltered::Bool
end

function Base.show(io::IO, smf::FixedHermiteSMF)
	println(io,"Fixed Order Hermite Stochastic Map Filter with filtered = $(smf.isfiltered)")
end

function (smf::FixedHermiteSMF)(X, ystar::Array{Float64,1}, t::Float64)
	Ny = smf.Ny
	Nx = smf.Nx

	Nystar = size(ystar, 1)
	Nypx, Ne = size(X)

	@assert Nystar == Ny "Size of ystar is not consistent with Ny"
	@assert Nypx == Ny + Nx "Size of X is not consistent with Ny or Nx"

	# Perturbation of the measurements
	smf.ϵy(X, 1, Ny)

	L = LinearTransform(X; diag = true)
	M = HermiteMap(smf.M.m, smf.M.Nx, L, smf.M.C)
	clearcoeff!(M)

	M = totalordermap(X, 2; b = "CstLinProHermite")

	optimize(M, X, nothing; withconstant = false, withqr = true,
			 verbose = false, start = Ny+1, P = serial, hessprecond = true)

	# Evaluate the transport map
	F = evaluate(M, X; apply_rescaling = true, start = Ny+1, P = serial)

	# Generate the posterior samples by partial inversion of the map
	hybridinverse!(X, F, M, ystar; start = Ny+1, P = serial)

	return X
end

#
# function assimilate_obs(M::HermiteMap, X, ystar, Ny, Nx; withconstant::Bool = false,
#                         withqr::Bool = false, verbose::Bool = false, P::Parallel = serial)
#
#         Nystar = size(ystar, 1)
#         Nypx, Ne = size(X)
#
#         @assert Nystar == Ny "Size of ystar is not consistent with Ny"
#         @assert Nypx == Ny + Nx "Size of X is not consistent with Ny or Nx"
#         @assert Ne == Neystar "Size of X and Ystar are not consistent"
#
#         # Optimize the transport map
#         M = optimize(M, X, "kfold"; withconstant = withconstant, withqr = withqr,
#                                verbose = verbose, start = Ny+1, P = P)
#
#         # Evaluate the transport map
#         F = evaluate(M, X; start = Ny+1, P = P)
#
#         inverse!(X, F, M, ystar; start = Ny+1, P = P)
#
#         return X
# end


text:
using Plots

# number of points to randomly generate in the square [(0,0), (0,1), (1,1), (1,0)]
N = 100_000

"function checks weather a given point is in the circle of radius equal to 1"
is_in( (x,y) ) = if x^2 + y^2 < 1
    1
else
    0
end

all_points = [ (rand(), rand()) for _ in 1:N]
points_in = [ p for p in all_points if is_in(p) == 1]
points_out = [ p for p in all_points if is_in(p) == 0]

myPI = 4 * length(points_in) / N

plot( [(0,0), (0,1), (1,1), (1,0), (0,0)],
    xlims = (0,1), ylims = (0,1),
    aspect_ratio=:equal,
    leg=nothing)
plot!( [(cos(x), sin(x)) for x in 0:0.01:2*pi])
scatter!(points_in)
scatter!(points_out)

print("PI = $(myPI)")

max_stars_repo_name: appblue/monte-carlo-pi-calculation

max_issues_repo_name: appblue/monte-carlo-pi-calculation

max_forks_repo_name: appblue/monte-carlo-pi-calculation


text:
# Code for 1D example with highly oscillating coefficient,
# using soft thresholding solver
using TensorTrains, TensorTrainFEM
using TensorTrains.Solvers

export multiscaleRefSol, stsolveOscillatory1D

function ttExMultiscaleD(L::Int64, δ::Float64)
	n = [1; 2*ones(Int64,L); 1]

	r = [1; 6*ones(Int64,L+1); 1]
	V = [ zeros(r[l],n[l],r[l+1]) for l ∈ 1:L+2 ]

	V[1][1,1,1] = 1/4
	V[1][1,1,2] = -1/4
	V[1][1,1,3] = -sin(2*π/δ)/8
	V[1][1,1,4] = cos(2*π/δ)/8
	V[1][1,1,5] = sin(2*π/δ)/8
	V[1][1,1,6] = -cos(2*π/δ)/8

	for l ∈ 1:L

		V[l+1][1,1,1] = 1
		V[l+1][1,2,1] = 1

		V[l+1][2,1,1] = -1/2
		V[l+1][2,1,2] = 1/2
		V[l+1][2,2,1] = 1/2
		V[l+1][2,2,2] = 1/2

		V[l+1][3,1,3] = cos(π/δ/2)
		V[l+1][3,1,4] = -sin(π/δ/2)
		V[l+1][4,1,3] = sin(π/δ/2)
		V[l+1][4,1,4] = cos(π/δ/2)

		V[l+1][3,1,3] = cos(π/δ/2)
		V[l+1][3,1,4] = -sin(π/δ/2)
		V[l+1][4,1,3] = sin(π/δ/2)
		V[l+1][4,1,4] = cos(π/δ/2)

		V[l+1][3,2,3] = cos(π/δ/2)
		V[l+1][3,2,4] = sin(π/δ/2)
		V[l+1][4,2,3] = -sin(π/δ/2)
		V[l+1][4,2,4] = cos(π/δ/2)

		V[l+1][5,1,3] = -cos(π/δ/2)/2
		V[l+1][5,1,4] = sin(π/δ/2)/2
		V[l+1][5,1,5] = cos(π/δ/2)/2
		V[l+1][5,1,6] = -sin(π/δ/2)/2
		V[l+1][6,1,3] = -sin(π/δ/2)/2
		V[l+1][6,1,4] = -cos(π/δ/2)/2
		V[l+1][6,1,5] = sin(π/δ/2)/2
		V[l+1][6,1,6] = cos(π/δ/2)/2

		V[l+1][5,2,3] = cos(π/δ/2)/2
		V[l+1][5,2,4] = sin(π/δ/2)/2
		V[l+1][5,2,5] = cos(π/δ/2)/2
		V[l+1][5,2,6] = sin(π/δ/2)/2
		V[l+1][6,2,3] = -sin(π/δ/2)/2
		V[l+1][6,2,4] = cos(π/δ/2)/2
		V[l+1][6,2,5] = -sin(π/δ/2)/2
		V[l+1][6,2,6] = cos(π/δ/2)/2

		δ = 2*δ;
	end

	V[L+2][1,1,1] = 1
	V[L+2][2,1,1] = 0
	V[L+2][3,1,1] = 0
	V[L+2][4,1,1] = 1
	V[L+2][5,1,1] = 0
	V[L+2][6,1,1] = 0

	return V
end

function ttExMultiscale(L::Int64, δ::Float64)
	n = [1; 2*ones(Int64,L); 1]

	r = [1; 7*ones(Int64,L+1); 1]
	U = [ zeros(r[l],n[l],r[l+1]) for l ∈ 1:L+2 ]

	U[1][1,1,1] = 1/6 + δ^2/(4*π)^2
	U[1][1,1,2] = 1/8
	U[1][1,1,3] = -1/24
	U[1][1,1,4] = (sin(π/δ)-δ*cos(π/δ)/π)/π/16
	U[1][1,1,5] = (cos(π/δ)-δ*sin(π/δ)/π )/π/16
	U[1][1,1,6] = -sin(π/δ)/π/16
	U[1][1,1,7] = -cos(π/δ)/π/16

	K = δ^(1/L)

	for l ∈ 1:L
		δ = 2*δ;

		U[l+1][1,1,1] = 1
		U[l+1][1,2,1] = 1

		U[l+1][2,1,1] = -1/2
		U[l+1][2,1,2] = 1/2
		U[l+1][2,2,1] = 1/2
		U[l+1][2,2,2] = 1/2

		U[l+1][3,1,2] = -3/4
		U[l+1][3,1,3] = 1/4
		U[l+1][3,2,2] = 3/4
		U[l+1][3,2,3] = 1/4

		U[l+1][4,1,4] = cos(π/δ)*K
		U[l+1][4,1,5] = sin(π/δ)*K
		U[l+1][5,1,4] = -sin(π/δ)*K
		U[l+1][5,1,5] = cos(π/δ)*K

		U[l+1][4,2,4] = cos(π/δ)*K
		U[l+1][4,2,5] = -sin(π/δ)*K
		U[l+1][5,2,4] = sin(π/δ)*K
		U[l+1][5,2,5] = cos(π/δ)*K

		U[l+1][6,1,4] = -cos(π/δ)/2*K
		U[l+1][6,1,5] = -sin(π/δ)/2*K
		U[l+1][6,1,6] = cos(π/δ)/2*K
		U[l+1][6,1,7] = sin(π/δ)/2*K
		U[l+1][7,1,4] = sin(π/δ)/2*K
		U[l+1][7,1,5] = -cos(π/δ)/2*K
		U[l+1][7,1,6] = -sin(π/δ)/2*K
		U[l+1][7,1,7] = cos(π/δ)/2*K

		U[l+1][6,2,4] = cos(π/δ)/2*K
		U[l+1][6,2,5] = -sin(π/δ)/2*K
		U[l+1][6,2,6] = cos(π/δ)/2*K
		U[l+1][6,2,7] = -sin(π/δ)/2*K
		U[l+1][7,2,4] = sin(π/δ)/2*K
		U[l+1][7,2,5] = cos(π/δ)/2*K
		U[l+1][7,2,6] = sin(π/δ)/2*K
		U[l+1][7,2,7] = cos(π/δ)/2*K

	end

	U[L+2][1,1,1] = 1
	U[L+2][2,1,1] = 1
	U[L+2][3,1,1] = 1
	U[L+2][4,1,1] = cos(π/δ)
	U[L+2][5,1,1] = sin(π/δ)
	U[L+2][6,1,1] = cos(π/δ)
	U[L+2][7,1,1] = sin(π/δ)

	return U
end

function multiscaleRefSol(L::Integer, K::Number)
    ttuex0 = ttExMultiscale(L, 2/K)
    dttuex0 = ttExMultiscaleD(L, 2/K)
    ttuex = Tensor(L)
    dttuex = Tensor(L)
    if L > 2
        ttuex[2:(L-1)] = ttuex0[3:L]
        dttuex[2:(L-1)] = dttuex0[3:L]
    end

    ttuex[1] = 2.0^(-L/2)*reshape(reshape(ttuex0[1], (1,7))*reshape(ttuex0[2], (7,2*7)), (1,2,7))
    ttuex[L] = reshape(reshape(ttuex0[L+1], (7*2,7))*reshape(ttuex0[L+2], (7,1)), (7,2,1))
    scale!(ttuex, 4.)

    dttuex[1] = 2.0^(-L/2)*reshape(reshape(dttuex0[1], (1,6))*reshape(dttuex0[2], (6,2*6)), (1,2,6))
    dttuex[L] = reshape(reshape(dttuex0[L+1], (6*2,6))*reshape(dttuex0[L+2], (6,1)), (6,2,1))
    scale!(dttuex, 4.)

    return ttuex, dttuex
end

# iteration with 1D bpx preconditioning,
# case with oscillatory coefficient
function stsolveOscillatory1D(L::Int64, ɛ::Float64, K::Number = 1)
  normest = 25.
  invnormest = .5 * 3. # taking into account min of coeff

  C = ttBPXdnC(L, 1)

  #F = ttBPXdnQ(L)
  F = ttBPXdnQLaplace(L, 1, 1)
  p = size(F[L],1)
  q = size(F[L],4)
  display(F)
  display((p,q))
  T = reshape(F[L],(p*4,q))
  S = reshape(F[L+1],(q,1))
  F[L] = reshape(T*S,(p,2,2,1))
  pop!(F)

  Ft = transpose(F)

  Dinv = ttdiagm(add(2., ttones(L), 1., ttcos(L, K*π, .5)))
  Drhs = ttones(L)
  d, _ = stsolve(L, 1, x->add(1., matvec(Dinv,x), -1., Drhs), 3., 1., .5*ɛ, .75)

  D = ttdiagm(d)

  println("operator ranks: ", ranks(F)," / ", ranks(D),
    " / ", ranks(Ft))

  f = ttrhsdn(L)
  g = matvec(C, f)
  svdtrunc!(g, svd!(g), .5*ɛ);
  println("rhs ranks ", ranks(g))

  Θ = 0.5

  function residual(u::Tensor, δ::Float64)
    gδ = deepcopy(g)
    svdtrunc!(gδ, svd!(gδ), .1*δ);
    t = matvec(F, u)
    maxrank = [maximum(ranks(t))]
    svdtrunc!(t, svd!(t), .4*δ / sqrt(normest));
    push!(maxrank, maximum(ranks(t)))
    s = matvec(D, t)
    svdtrunc!(s, svd!(s), .4*δ / sqrt(normest));
    push!(maxrank, maximum(ranks(s)))
    r = matvec(Ft, s)
    push!(maxrank, maximum(ranks(r)))
    add!(r, -1., g)
    push!(maxrank, maximum(ranks(r)))
    svdtrunc!(r, svd!(r), .1*δ);
    push!(maxrank, maximum(ranks(r)))
    return r, maxrank
  end

  u, cdata = inexstsolve(L, 1, residual, normest, invnormest, ɛ, Θ)

  return u, matvec(C, u), cdata
end


text:
function CGStep(m1,m2;a=1.,b=1.)
	# m1 = a*m1 + b*m2
	return a*m1 + b*m2
end

function CGStep(m1,m2,h1::Array{Header,1},h2::Array{Header,1};a=1.,b=1.)
	# m1 = a*m1 + b*m2
	return a*m1 + b*m2,h1
end

function CGStep(m1::AbstractString,m2::AbstractString;a=1.,b=1.)
	# m1 = a*m1 + b*m2

	#tmp = join(["tmp_CGStep_",string(round(Int,rand()*100000))])
	#@compat params = Dict(:a=>a,:b=>b)
	#SeisProcess(m1,m2,tmp,[CGStep],[params],key=["imx";"imy"])
	#SeisCopy(tmp,m1)
	#SeisRemove(tmp)
	d1,h,ext = SeisRead(m1)
	d2,h,ext = SeisRead(m2)
	SeisWrite(m1,a*d1[:,:] + b*d2[:,:],h,ext)
end

function CGStep(m1::Array{AbstractString,1},m2::Array{AbstractString,1};a=[1. 1. 1.],b=[1. 1. 1.])
	# m1 = a*m1 + b*m2
	for j = 1 : length(m1)
		CGStep(m1[j],m2[j],a=a[j],b=b[j])
	end

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
# Contents: Functions for creating/handling secure cookies.

################################################################################
# Create session cookie
#
# The scheme:
#     const_key, const_iv     = global constants, output from a cryptographically secure random number generator
#                               (used to encrypt session-specific secret keys)
#     timestamp               = milliseconds since epoch, represented as a string
#     session_key, session_iv = output from a cryptographic random number generator, unique for each session
#     encrypted_session_key   = AES CBC encrypt(const_key, const_iv, session_key)
#     data blob               = AES CBC encrypt(session_key, session_iv, arbitrary data)
#     hmac signature          = HMAC(session_key, timestamp * data_blob)
#     unencoded cookie_value  = session_iv * encrypted_secret_key * hmac signature * timestamp * data blob
#     cookie_value            = base64encode(unencoded cookie value)
#                               (the encoding is for transport in an http header)
#
################################################################################

"""
Create a secure session cookie for the response.
The cookie value includes the encryption of the supplied data.
The Secure and HttpOnly attributes are set according to global variables.
"""
function create_secure_session_cookie(data, res, cookie_name = "sessionid")
    cookie_value = create_secure_session_cookievalue(data)
    attr = Dict("Max-Age" => timeout_str)
    if encrypted_sessions_only
	attr["Secure"] = ""
    end
    if http_only
	attr["HttpOnly"] = ""
    end
    setcookie!(res, cookie_name, string(cookie_value), attr)
end

"""
Create the value of the secure session cookie.

Input:  Data (ASCIIString) to be embedded in the encrypted cookie value.
Output: Cookie value (ASCIIString)

Note: Binary data is base64 encoded for transport in http headers (base64 is more space efficient than hex encoding).
"""
function create_secure_session_cookievalue(plaintext)
    # Encrypt data
    session_key = csrng(key_length)
    session_iv  = csrng(block_size)
    data_blob   = encrypt(CIPHER_AES, session_key, plaintext, session_iv)    # Encryption is done in CBC mode

    # Compute HMAC signature
    timestamp      = string(get_timestamp())    # Millieconds since epoch
    ts_uint8       = convert(Array{UInt8, 1}, timestamp)
    hmac_signature = digest(MD_SHA256, vcat(ts_uint8, data_blob), session_key)

    # Compute cookie value
    encrypted_session_key = encrypt(CIPHER_AES, const_key, session_key, const_iv)
    cookie_value          = base64encode(vcat(session_iv, encrypted_session_key, hmac_signature, ts_uint8, data_blob))
    cookie_value
end

################################################################################
# Validate session cookie
################################################################################

# The application validates a session cookie as follows:
# 1. Decode hmac signature.
# 2. Compute HMAC(secret key, timestamp * data blob) and compare it to hmac signature. Fail if they differ.
# 3. Decode timestamp.
# 4. Verify that the current time in milliseconds since the epoch is not greater than timestamp + session timeout.
# TODO: If the cookie is not valid, the application must refuse the requested action and redirect the user to the login page.
#
# If the cookie is valid, the application can
# 1. Decrypt data blob
# 2. Parse or deserialize data blob as appropriate.
# At this point, the application has a valid state object for the user's session and can proceed with processing the requested action.

"""
Returns the decrypted cookie data.

Returns "" if the cookie doesn't exist.
"""
function get_session_cookie_data(req, cookie_name)
    result = ""
    if haskey(req.headers, "Cookie")
	cookie_value = get_cookie_value(req, cookie_name)
	if cookie_value != ""
	    cookie_is_valid, data_blob, session_key, session_iv = session_cookie_is_valid(cookie_value)
	    if cookie_is_valid
		result = decrypt(CIPHER_AES, session_key, data_blob, session_iv)
		result = String(result)
	    end
	end
    end
    result
end

"""
Returns the cookie value, which is encrypted.

Returns "" if the cookie doesn't exist.
"""
function get_cookie_value(req, cookie_name)
    cookie_value = ""
    ckie         = req.headers["Cookie"]    # ASCIIString: "name1=value1; name2=value2"
    names_values = split(ckie, ";")         # "name=value"
    for nv in names_values
	r = search(nv, cookie_name)         # first_idx:last_idx
	if length(r) > 0                    # cookie_name is in nv
	    r2           = search(nv, "=")
	    cookie_value = nv[(r2[1] + 1):end]
	    break
	end
    end
    String(cookie_value)      # Convert SubString to string for base64 decoding
end

"""
Returns: cookie_is_valid (Bool) and session data.

cookie_is_valid is true if session cookie:
1) Has not expired, and
2) hmac_signature == HMAC(secret key, timestamp * data_blob)
"""
function session_cookie_is_valid(cookie_value)
    # Extract cookie data
    cookie_value   = base64decode(cookie_value)
    session_iv     = cookie_value[1:block_size]
    offset         = block_size
    encrypted_session_key = cookie_value[(offset + 1):(offset + key_length + block_size)]
    session_key    = decrypt(CIPHER_AES, const_key, encrypted_session_key, const_iv)
    offset += key_length + block_size
    hmac_signature = view(cookie_value, (offset + 1):(offset + key_length))
    offset += key_length
    ts_uint8       = cookie_value[(offset + 1):(offset + 13)]
    timestamp      = parse(Int,  String(ts_uint8))    # Seconds since epoch
    offset += 13
    data_blob      = cookie_value[(offset + 1):end]

    # Determine conditions
    current_time = get_timestamp()
    expired      = current_time > timestamp + session_timeout
    hmac_sig2    = digest(MD_SHA256, vcat(ts_uint8, data_blob), session_key)
    hmac_ok      = hmac_sig2 == hmac_signature

    # Prepare results
    cookie_is_valid = false
    if !expired && hmac_ok
        cookie_is_valid = true
    end
    cookie_is_valid, data_blob, session_key, session_iv
end

"""
Invalidates the cookie with name == cookie_name.
Curently this works by setting the Max-Age to 0.
"""
function invalidate_cookie!(res, cookie_name)
    setcookie!(res, cookie_name, "", Dict("Max-Age" => "0"))
end

# EOF


text:
using BHAPtfem
using Test

#=
Compare formulas at:
http://www.awc.org/pdf/codes-standards/publications/design-aids/AWC-DA6-BeamFormulas-0710.pdf
=#

println("Test when editing in Atom\n")

N = 20
Np1 = N + 1
Nhp1 = Int(N/2) + 1

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(N, Np1, 3, 1, 1, Line(2, 3)),
  :properties => [1.0e6 1.0e6 1.0e6 3.0e5;],
  :x_coords => range(0, stop=4, length=Np1),
  :y_coords => zeros(Np1),
  :z_coords => zeros(Np1),
  :g_num => [
    collect(1:N)';
    collect(2:Np1)'],
  :support => [
    (1, [1 0 0 0 0 0]),
    (Np1, [0 0 0 0 0 0]),
    ],
  :loaded_nodes => [
    (Nhp1, [0.0 -10000.0 0.0 0.0 0.0 0.0])]
)

@time m, dis_df, fm_df = PtFEM.p44(data)

# See figure 24 in above reference (Δmax):
@test round.(m.displacements[2,Nhp1], digits=11) ≈ -0.00333333333

# See figure 24 in above reference (Mmax):
@test m.actions[12,Nhp1] ≈ 3999.9999999975844


text:
using cascadeddecay
using PyCall
@pyimport matplotlib.pyplot as plt

const A = [0.01:0.01:5]
const γ = 1.
#const Θ = pi/2

plt.figure(1, figsize=(6,4))
plt.figure(2, figsize=(6,4))

for Θ=range(0, stop=pi/2, length=11)
    Ω = [Omega(a, Θ, γ) for a=A]
    Γ = [Gamma(a, Θ, γ) for a=A]

    plt.figure(1)
    plt.plot(A, Ω, "b")
    plt.ylim(-0.5,0.5)

    plt.figure(2)
    plt.plot(A, Γ, "g")
    plt.ylim(-0.5,0.5)
end

plt.figure(1)
plt.ylim(-0.5,0.5)
plt.xlabel("a")
plt.ylabel("\$\\Omega(a,\\Theta)\$")
plt.savefig("images/collective_interaction_omega.pdf")

plt.figure(2)
plt.ylim(-0.5,0.5)
plt.xlabel("a")
plt.ylabel("\$\\Gamma(a,\\Theta)\$")
plt.savefig("images/collective_interaction_gamma.pdf")

plt.show()


text:
eps(::Type{Posit8}) = Posit8(0x02)
eps(::Type{Posit16}) = Posit16(0x0100)
eps(::Type{Posit32}) = Posit32(0x00a0_0000)

#TODO remove subtraction once pX1_sub is available in the C version
eps(::Type{Posit8_1}) = Posit8_1(0x41) - Posit8_1(0x40)
eps(::Type{Posit16_1}) = Posit16_1(0x4001) - Posit16_1(0x4000)
eps(::Type{Posit24_1}) = Posit24_1(0x4000_01) - Posit24_1(0x4000_00)

eps(::Type{Posit8_2}) = Posit8_2(0x2800_0000)
eps(::Type{Posit16_2}) = Posit16_2(0x0a00_0000)
eps(::Type{Posit24_2}) = Posit24_2(0x0280_0000)

eps(p::AbstractPosit) = max(p-prevfloat(p),nextfloat(p)-p)


text:
#########################################################################################################################################
######### nurbs.jl contains all basic functionality in order to compute with n-variate NURBS and their derivatives. ############
#########################################################################################################################################
#
# Types:
#       Point                  -
#       BezierPatch            -
#       NURBS                  -
#
# Functions:
#       nsdims()               -  determine the space dimensions of the BezierPatch or NURBS object
#       boundary()             -  return boundary Bezier or NURBS objects
#       ∂()                    -  short for boundary()
#       dimsplinespace()       -  number of degrees of freedom of NURBS object
#       nbezierpatches()       -  number of Bezier elements
#       bezierdecompoperator() -  return the Bezier decomposition operators for each spatial direction of the NURBS object
#       ienarray()             -  return IEN array of NURBS object
#       refine!()              -  knot insertion of a specific knot in a specific direction
#       globalrefine!()        -  repeated knot insertion in each direction
#       degreeelevate!()       -  degree elevate NURBS object in a specific spatial direction
#
#########################################################################################################################################
#
# File part of the Julia IsoGeometric Analysis toolbox
# R.R.Hiemstra
# 19-04-2016
#
#########################################################################################################################################

###################################################################################################################
############################################ Type definition BezierPatch ##########################################
###################################################################################################################

const Point = Vector{Float64}

"type definition of Bezier patch"
struct BezierPatch{N}
    deg    :: NTuple{N,Degree}
    cpts   :: Array{Float64}
    wpts   :: Vector{Float64}
    dec    :: NTuple{N,Matrix{Float64}}
end

"contruct a NURBS patch in parametric space"
function BezierPatch(p::NTuple{N,Degree}) where {N}
    cpts = allcomb(ntuple(k -> grevillepoints(p[k],buildvector([-1.0;1.0], [p[k]+1,p[k]+1])), N)...)
    wpts = ones(Float64,size(cpts,1))
    C    = ntuple(k -> eye(p[k]+1), N)
    return BezierPatch(p, cpts, wpts, C)
end

############################################ Elementary operations ################################################

nsdims(S::BezierPatch) = length(S.deg)

function boundary(S::BezierPatch{1}, k::Int)
    if k==1
        return vec(S.cpts[1,:])
    elseif k==2
        return vec(S.cpts[end,:])
    else
        error("Out of bounds")
    end
end

"boundary of a 2d Bezier Patch"
function boundary(S::BezierPatch{N}, dir::Int, k::Int) where {N}

    # initialize
    dims = S.deg+1

    # local numbering
    # I = dropdims(slicedim(reshape(1:prod(dims),dims), dir, k), dims=dir)[:]
    I = slicedim(reshape(1:prod(dims),dims), dir, k)
    J = 1:N.!=dir

    # return boundary k
    return BezierPatch(S.deg[J], S.cpts[I,:], S.wpts[I], S.dec[J])
end

# entire boundary of
boundary(S::BezierPatch{1}) = [∂(S,1) ∂(S,2)]
boundary(S::BezierPatch{N}) where {N} = BezierPatch{N-1}[∂(S,dir,k) for k in 1:2, dir in 1:N][:]

# short name
∂(S::BezierPatch, dir::Int, k::Int) = boundary(S, dir, k)
∂(S::BezierPatch, k::Int) = boundary(S, k)
∂(S::BezierPatch) = boundary(S)

###################################################################################################################
############################################ Type definition NURBS ###############################################
###################################################################################################################

"type definition of NURBS patch"
mutable struct NURBS{N}
    deg    :: Vector{Degree}
    kts    :: Vector{KnotVector}
    cpts   :: Matrix{Float64}
    wpts   :: Vector{Float64}
    dec    :: NTuple{N,Array{Float64}}
    ien    :: Array{Int64}
    id     :: Array{Int64}
end
NURBS(p, kts, cpts, wpts) = NURBS(p, kts, cpts, wpts, ntuple(k->Float64[], length(p)),Int64[], Int64[])

"contruct a NURBS patch in parametric space"
function NURBS(p, kts)
    cpts = allcomb(ntuple(k -> grevillepoints(p[k],kts[k]), length(p))...)
    wpts = ones(Float64,size(cpts,1))
    return NURBS(p, kts, cpts, wpts)
end

# dimension of patch
nsdim(S::NURBS) = length(S.deg)
dimsplinespace(S::NURBS) = ntuple(k -> dimsplinespace(S.deg[k],S.kts[k]), nsdim(S))
nbezierpatches(S::NURBS) = ntuple(k -> length(unique(S.kts[k]))-1, nsdim(S))

"bezier decomposition operators"
bezierdecompoperator(S::NURBS) = ntuple(k -> bezierdecompoperator(S.deg[k],S.kts[k]), nsdim(S))

"construct IEN array for single NURBS patch S"
function ienarray(S::NURBS)
    # initialize
    p     = S.deg
    nsd   = length(p)
    ndofs = dimsplinespace(S)
    ukts, umult, inda, indc = uniqueknots(S.kts)
    nelms = ntuple(k -> length(ukts[k])-1, nsd)

    # allocate space for IEN array
    m = prod(nelms)     # determine number of elements
    n = prod(S.deg.+1)  # number of Bezier degrees of freedom

    # construct IEN-array
    IEN = zeros(Int64,n,m)
    for k in 1:m
        subs = CartesianIndices(nelms)[k]
        span = ntuple(i -> indc[i][subs[i]+1]-1, nsd)
        ind  = allcomb(ntuple(i -> collect(span[i]-p[i]:span[i]), nsd)...)
        IEN[:,k] = (LinearIndices(ndofs))[CartesianIndex.(ntuple(i->ind[:,i], nsd)...)]
        # IEN[:,k] = sub2ind(ndofs, ntuple(i->ind[:,i], nsd)...)
    end
    return IEN
end

"refine NURBS patch"
function refine!(S::NURBS, dir::Int, u::Vector{Float64})
    # initialize
    p     = S.deg
    nsd   = nsdim(S)
    dims  = dimsplinespace(S)
    perm  = [dir:nsd; 1:dir-1]

    # compute refinement operator
    C, nkts = knotinsertionoperator(S.deg[dir],S.kts[dir],u)

    # compute new control points
    m  = size(S.cpts,2)
    n  = prod(dims[1:nsd.!=dir])*dimsplinespace(p[dir],nkts)
    Pw = zeros(n,m); [Pw[:,k] = permutedims(C * permutedims(reshape(S.wpts.*S.cpts[:,k],dims),perm), invperm(perm))[:] for k in 1:m]
    W  = permutedims(C * permutedims(reshape(S.wpts,dims),perm), invperm(perm))[:]

    # update data
    S.kts[dir] = nkts
    S.cpts = Pw./W
    S.wpts = W;
end

"global knot insertion in direction dir"
function globalrefine!(S::NURBS, dir::Int, k::Int)
    u = linearspace(unique(S.kts[dir]),k+2)[2:end-1,:][:]
    refine!(S, dir, u)
end

"global knot insertion in direction dir"
function globalrefine!(S::NURBS, k::Int)
    for i in 1:nsdim(S)
        u = linearspace(unique(S.kts[i]),k+2)[2:end-1,:][:]
        refine!(S, i, u)
    end
end

"global knot insertion in direction dir"
function degreeelevate!(S::NURBS, dir::Int, t::Int)
    # initialize
    nsd   = nsdim(S)
    dims  = dimsplinespace(S)
    perm  = [dir:nsd; 1:dir-1]

    # compute refinement operator
    C, nkts, p = degreeelevationoperator(S.deg[dir],S.kts[dir],t)

    # compute new control points
    m  = size(S.cpts,2)
    n  = prod(dims[1:nsd.!=dir])*dimsplinespace(p,nkts)
    Pw = zeros(n,m); [Pw[:,k] = permutedims(C * permutedims(reshape(S.wpts.*S.cpts[:,k],dims),perm), invperm(perm))[:] for k in 1:m]
    W  = permutedims(C * permutedims(reshape(S.wpts,dims),perm), invperm(perm))[:]

    # update data
    S.deg[dir] = p
    S.kts[dir] = nkts
    S.cpts = Pw./W;
    S.wpts = W;

end
degreeelevate!(S::NURBS,t::Int) = [degreeelevate!(S,dir,t) for dir in 1:nsdim(S)]

############################################ Boundary routines ##############################################

"boundary of a 2d Bezier Patch"
function boundary(S::NURBS{N}, dir::Int, k::Int) where {N}
    # initialize
    dims = dimsplinespace(S)

    # local numbering
    I = selectdim(reshape(1:prod(dims),dims), dir, k==1 ? k : dims[dir])
    J = 1:length(S.deg).!=dir

    # set degree, controlpoints and weights of boundary (dir,k)
    p    = S.deg[J]
    kts  = S.kts[J]
    cpts = S.cpts[I,:]
    wpts = S.wpts[I]
    dec  = S.dec[J]

    # set IEN-array if non-empty
    ien = Int64[]
    if isempty(S.ien)==false

        # determine boundary dofs
        dims2 = S.deg .+ 1
        I     = selectdim(reshape(1:prod(dims2),dims2...), dir, k==1 ? k : dims2[dir])

        # determine boundary elements
        dims3 = nbezierpatches(S)
        J     = selectdim(reshape(1:prod(dims3),dims3...), dir, k==1 ? k : dims3[dir])

        ien = S.ien[I,J]
    end

    # return boundary k
    return NURBS(p, kts, cpts, wpts, dec, ien, Int64[])
end

function boundary(S::NURBS{1}, dir::Int, k::Int)
    if k==1
        return vec(S.cpts[1,:])
    elseif k==2
        return vec(S.cpts[end,:])
    end
end
boundary(S::NURBS{1}) = [∂(S,1,1) ∂(S,1,2)]
boundary(S::NURBS{N}) where {N} = NURBS{N-1}[∂(S,dir,k) for k in 1:2, dir in 1:N][:]

∂(S::NURBS, dir::Int, k::Int) = boundary(S, dir, k)
∂(S::NURBS, k::Int) = boundary(S, k)
∂(S::NURBS) = boundary(S)

max_stars_repo_name: michakraus/iga-demos

max_issues_repo_name: michakraus/iga-demos

max_forks_repo_name: michakraus/iga-demos


text:
using LinearAlgebra, CUDA, Flux, JLD2
import KitML

cd(@__DIR__)
data = KitML.load_data("2021-02-20_09:30:26.csv"; mode = :csv, dlm = ",")

X = hcat(data.u_0, data.u_1, data.u_2, data.u_3) |> permutedims
Y = data.h |> permutedims

data = nothing # deallocate

#nn = KitML.ICNNChain(4, 1, [16, 32, 32, 8], tanh)
@load "model.jld2" nn

loss(x, y) = sum(abs2, nn(x) - y) / 15968000

for iter = 1:1000
    KitML.sci_train!(nn, (X, Y), ADAM(); epoch = 2, batch = 4000, device = gpu)

    if iter % 2 == 0
        @show loss(X, Y)
        KitML.save_model(nn)
    end
end

#fnn = KitML.FastICNN(4, 1, [16, 64, 64, 8], tanh)
#res = KitML.sci_train(fnn, (X, Y); device = cpu, maxiters = 200)


text:
using JSONTypeProvider
using Test
import JSON3

@testset "JSONTypeProvider.jl" begin
    json = JSON3.read(read("menu.json", String)) # either a JSON.Array or JSON.Object

    # build a type for the JSON
    raw_json_type = JSONTypeProvider.build_type(json)
    @test raw_json_type <: NamedTuple

    # turn the type into struct expressions, including replacing sub types with references to a struct
    json_exprs = JSONTypeProvider.to_exprs(raw_json_type, :MyStruct)
    @test length(json_exprs) == 3

    # write the types to a file, then can be edited/included as needed
    path = mktempdir()
    file = joinpath(path, "test.jl")
    JSONTypeProvider.write_exprs(json_exprs, file)
    lines = readlines(file)
    @test length(lines) > 3

    json_arr = JSON3.read(read("menu.json", String)) # either a JSON.Array or JSON.Object

    # build a type for the JSON
    raw_json_arr_type = JSONTypeProvider.build_type(json_arr)
    @test raw_json_type <: NamedTuple

    # turn the type into struct expressions, including replacing sub types with references to a struct
    json_arr_exprs = JSONTypeProvider.to_exprs(raw_json_arr_type, :MyStruct)
    @test length(json_exprs) == 3

    # write the types to a file, then can be edited/included as needed
    file_arr = joinpath(path, "test_arr.jl")
    JSONTypeProvider.write_exprs(json_arr_exprs, file_arr)
    arr_lines = readlines(file_arr)
    @test length(lines) > 3
    @test arr_lines == lines

end


text:
@testset "Toffoli test" begin
	N = 3
	state = MPSState(N)
	apply!(state, X(1))
	apply!(state, X(2))
	state_e = toExactState(state)
	circuit = QCircuit(N)
	add!(circuit,     H(3),
					  CNOT(2,3),
					  TDAG(3),
					  CNOT(1,3),
					  T(3),
					  CNOT(2,3),
					  TDAG(3),
					  CNOT(1,3),
					  T(3),
					  H(3),
					  TDAG(2),
					  CNOT(1,2),
					  TDAG(2),
					  CNOT(1,2),
					  T(1),
					  S(2))
	apply!(state, circuit)
	apply!(state_e, TOFFOLI(1,2,3))
	RoundToZero # set rounding mode
	for i = 0:1
		for j = 0:1
			for k = 0:1
				@test isapprox(round(probability(state, [i,j,k])),probability(state_e, [i,j,k]))
			end
		end
	end
end


text:
#Using FLoops to parallelize findNextEquivClass!
using FLoops

#Generate an iterator equivalent to combinations(x,2) but with base packages
#The combinations are not in order, but that doesn't matter in this case
#benefit: packages have more methods defined for the base iterators
allpairs(v) = Iterators.filter(i -> isless(i...), Iterators.product(v,v))

function findNextEquivClass!(newStep, dataParsed, g, findBestOperation::Function, debug)

    #Parallelization with synced threads (to avoid race conditions)
    @sync begin
        #Loop through all possible node combinations
        for (x,y) in combinations(vertices(g),2)
            #Spawn a new thread
            Threads.@spawn begin
                #Check if there is an edge between two vertices
                hasEdge = isadjacent(g,x,y)

                #Skip over adjacent edges if we're trying to insert and vice versa
                if findBestOperation == findBestInsert ? !hasEdge : hasEdge
                    #For this pair of nodes, the following function will:
                        #(1) check if a valid operator exists
                        #(2) score all valid operators and return the one with the highest score
                        #findBestOperation is either "findBestInsert" or "findBestDelete"
                    (bestSubset, bestScore) = findBestOperation(dataParsed, g, x, y, debug)

                    #If the best valid operator was better than any previously found...
                    if bestScore > newStep.Δscore
                        #...update newStep
                        newStep.edge = Edge(x,y,true)
                        newStep.subset = bestSubset
                        newStep.Δscore = bestScore

                        debug && message(newStep)
                    end
                end
            end
        end
    end
end

ThreadsX.maximum((f(x,y),x,y) for (x,y) in allpairs(v))

maximum(allpairs(v)) do (x,y)
    (f(x,y),x,y)
end

function testPar!(newStep, dataParsed, g, findBestOperation::Function, debug)

    for (x,y) in allpairs(vertices(g))

    end

end


text:
include("scenes.jl")
include("path_planner.jl")
include("model.jl")

using Printf: @sprintf
import Random

function make_scene()
    scene = Scene(0, 1, 0, 1)
    add!(scene, Tree(Point(0.30, 0.20), size=0.1))
    add!(scene, Tree(Point(0.83, 0.80), size=0.1))
    add!(scene, Tree(Point(0.80, 0.40), size=0.1))
    horiz = 1
    vert = 2
    wall_height = 0.30
    wall_thickness = 0.02
    walls = [
        Wall(Point(0.20, 0.40), horiz, 0.40, wall_thickness, wall_height)
        Wall(Point(0.60, 0.40), vert, 0.40, wall_thickness, wall_height)
        Wall(Point(0.60 - 0.15, 0.80), horiz, 0.15 + wall_thickness, wall_thickness, wall_height)
        Wall(Point(0.20, 0.80), horiz, 0.15, wall_thickness, wall_height)
        Wall(Point(0.20, 0.40), vert, 0.40, wall_thickness, wall_height)]
    for wall in walls
        add!(scene, wall)
    end
    return scene
end

const scene = make_scene()
const times = collect(range(0, stop=1, length=20))

function show_prior_samples()
    Random.seed!(0)
    figure(figsize=(32, 32))
    for i=1:15
        subplot(4, 4, i)
        ax = gca()
        trace = simulate(model, (scene, times))
        render(scene, trace, ax)
    end
    savefig("demo.png")
end
show_prior_samples()

@gen function stop_proposal(prev_trace)
    @trace(uniform(0, 1), :stop_x)
    @trace(uniform(0, 1), :stop_y)
end

@gen function speed_proposal(prev_trace)
    @trace(uniform(0, 1), :speed)
end

@gen function noise_proposal(prev_trace)
    @trace(uniform(0, 1), :noise)
end

function inference(measurements::Vector{Point}, start::Point, iters::Int)
    t = length(measurements)

    constraints = choicemap()
    for (i, pt) in enumerate(measurements)
        constraints[i => :x] = pt.x
        constraints[i => :y] = pt.y
    end
    constraints[:start_x] = start.x
    constraints[:start_y] = start.y

    (trace, _) = generate(model, (scene, times[1:t]), constraints)

    for iter=1:iters
        (trace, _) = metropolis_hastings(trace, stop_proposal, ())
        (trace, _) = metropolis_hastings(trace, speed_proposal, ())
        (trace, _) = metropolis_hastings(trace, noise_proposal, ())
    end

    return trace
end

function experiment()

    # generate simulated ground truth
    Random.seed!(0)
    constraints = choicemap()
    constraints[:start_x] = 0.1
    constraints[:start_y] = 0.1
    constraints[:stop_x] = 0.5
    constraints[:stop_y] = 0.5
    constraints[:noise] = 0.1
    (trace, _) = generate(model, (scene, times), constraints)

    figure(figsize=(4, 4))
    ax = gca()
    render(scene, trace, ax)
    savefig("ground_truth.png")

    assignment = get_choices(trace)
    measurements = [Point(assignment[i => :x], assignment[i => :y]) for i=1:length(times)]
    start = Point(assignment[:start_x], assignment[:start_y])

    #for t=1:length(times)
    for t in [length(times)]
        println("t: $t")
        traces = []
        for i=1:100
            println(i)
            @time trace = inference(measurements[1:t], start, 100)
            push!(traces, trace)
        end
        figure(figsize=(4, 4))
        ax = gca()
        for (i, trace) in enumerate(traces)
            render(scene, trace, ax; show_measurements=i>1, show_start=i>1,
                   show_path=false, show_noise=false, stop_alpha=0.2, path_alpha=0.2)
        end
        fname = @sprintf("inferred_%03d.png", t)
        savefig(fname)
    end
end

experiment()


text:
module Pipelines

using DataFrames
using Random

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.Utils

import AutoMLPipeline.AbsTypes: fit!, transform!
export fit!, transform!
export LinearPipeline, ComboPipeline, @pipeline, @pipelinex

mutable struct LinearPipeline <: Workflow
  name::String
  model::Dict
  args::Dict

  function LinearPipeline(args::Dict = Dict())
    default_args = Dict(
			:name => "linearpipeline",
			# machines as list to chain in sequence.
			:machines => Vector{Machine}(),
			# Transformer args as list applied to same index transformer.
			:machine_args => Dict()
		       )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

function LinearPipeline(machs::Vector{T}) where {T<:Machine}
  LinearPipeline(Dict(:machines => machs))
end

function LinearPipeline(machs...)
  combo=nothing
  if eltype(machs) <: Machine
    v=[x for x in machs] # convert tuples to vector
    combo = LinearPipeline(v)
  else
    error("argument setup error")
  end
  return combo
end

function fit!(pipe::LinearPipeline, features::DataFrame, labels::Vector=[])
  instances=deepcopy(features)
  machines = pipe.args[:machines]
  machine_args = pipe.args[:machine_args]

  current_instances = instances
  new_machines = Machine[]

  # fit-transform all except last element
  # last element calls fit only
  trlength = length(machines)
  for t_index in 1:(trlength - 1)
    machine = createmachine(machines[t_index], machine_args)
    push!(new_machines, machine)
    fit!(machine, current_instances, labels)
    current_instances = transform!(machine, current_instances)
  end
  machine = createmachine(machines[trlength], machine_args)
  push!(new_machines, machine)
  fit!(machine, current_instances, labels)

  pipe.model = Dict(
      :machines => new_machines,
      :machine_args => machine_args
  )
end

function transform!(pipe::LinearPipeline, instances::DataFrame)
  machines = pipe.model[:machines]

  current_instances = deepcopy(instances)
  for t_index in 1:length(machines)
    machine = machines[t_index]
    current_instances = transform!(machine, current_instances)
  end

  return current_instances
end

mutable struct ComboPipeline <: Workflow
  name::String
  model::Dict
  args::Dict

  function ComboPipeline(args::Dict = Dict())
    default_args = Dict(
			:name => "combopipeline",
			# machines as list to chain in sequence.
			:machines => Vector{Machine}(),
			# Transformer args as list applied to same index transformer.
			:machine_args => Dict()
		       )
    cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

function ComboPipeline(machs::Vector{T}) where {T<:Machine}
  ComboPipeline(Dict(:machines => machs))
end

function ComboPipeline(machs...)
  combo=nothing
  if eltype(machs) <: Machine
    v=[eval(x) for x in machs] # convert tuples to vector
    combo = ComboPipeline(v)
  else
    error("argument setup error")
  end
  return combo
end

function fit!(pipe::ComboPipeline, features::DataFrame, labels::Vector=[])
  instances=deepcopy(features)
  machines = pipe.args[:machines]
  machine_args = pipe.args[:machine_args]

  new_machines = Machine[]
  new_instances = DataFrame()
  trlength = length(machines)
  for t_index in eachindex(machines)
    machine = createmachine(machines[t_index], machine_args)
    push!(new_machines, machine)
    fit!(machine, instances, labels)
  end

  pipe.model = Dict(
      :machines => new_machines,
      :machine_args => machine_args
  )
end

function transform!(pipe::ComboPipeline, features::DataFrame)
  machines = pipe.model[:machines]
  instances = deepcopy(features)
  new_instances = DataFrame()
  for t_index in eachindex(machines)
    machine = machines[t_index]
    current_instances = transform!(machine, instances)
    new_instances = hcat(new_instances,current_instances,makeunique=true)
  end

  return new_instances
end

function processexpr(args)
  for ndx in eachindex(args)
    if typeof(args[ndx]) == Expr
      processexpr(args[ndx].args)
    elseif args[ndx] == :+
      args[ndx] = :LinearPipeline
    elseif args[ndx] == :*
      args[ndx] = :ComboPipeline
    #else
    #  esc(:(args[ndx] = eval(args[ndx]))) # refer to local variable
    #  #args[ndx] = eval(:($(args[ndx])))
    #  #args[ndx] = eval(args[ndx])
    #  #println(args[ndx])
    end
  end
  return args
end

macro pipeline(expr)
  lexpr = :($(esc(expr)))
  res = processexpr(lexpr.args)
  lexpr.args = res
  lexpr
end

macro pipelinex(expr)
  lexpr = :($(esc(expr)))
  res = processexpr(lexpr.args)
  res
end

end


text:
module Laplacian

# Create a Kronecker Delta matrix with the given relationship and dimensions
function δ(fi::Function, fj::Function, n::Integer)
    r = zeros(n,n)
    for j in [1:n]
        for i in [1:n]
            fi(i) == fj(j) && (r[i,j] = 1.0)
        end
    end
    return r
end

kronecker = δ

# Create a laplacian matrix with specified grid size and grid divisions
function lap(N::Integer, L::Real)

    Np1sq = (N+1)*(N+1)                 # have (N+1)^2 grid points
    a = -4.0δ(i->i, j->j, Np1sq)
    a += δ(i->mod1(i+1, Np1sq), j->j, Np1sq)
    a += δ(i->mod1(i-1, Np1sq), j->j, Np1sq)
    a += δ(i->mod1(i+N+1, Np1sq), j->j, Np1sq)
    a += δ(i->mod1(i-N-1, Np1sq), j->j, Np1sq)
    return a * (N*N / (L*L))

end

lap(N::Integer) = lap(N,1)

# Split a laplacian matrix into free and boundary parts
function opsplit(A::Array, β::Array{Int,1})

    ndims(A) == 2 || error("operator matrix must be 2-d")
    length(β) == size(A)[1] || error("specify which indices are boundary conditions")
    size(A)[1] == size(A)[2] || error("operator matrix must be square")

    B = zeros(size(A))
    F = zeros(size(A))
    # Set the boundary condition columns equal to the corresponding columns of A
    for i in 1:length(β)
        β[i] == 1 ? B[:,i] = A[:,i] : F[:,i] = A[:,i]
    end
    return F, B

end

# β[i] = 1 iff ϕ[i] is a boundary condition
function rhoeff(ρ::Array, Δb::Array, ϕb::Array, β::Array{Int,1})
    return (ρ + Δb*ϕb) .* (1 .- β)
end

# Cancel out bound rows to get Δeffective
function deltaeff(Δf, β::Array{Int,1})
    return Δf .* (1 .- β)
end

end


text:
n = rand(2:10)

A = matrixdepot("wilkinson", n)

@test A == matrixdepot("wilkinson", Float64, n)

@test issym(full(A))

@test_approx_eq diag(full(A),1) ones(Float64, n - 1)

# symmetric along antidiagonal
@test issym(full(A)[n:-1:1,:])

θ = matrixdepot("wilkinson", 1)

println("'wilkinson' passed test...")


text:
using Documenter
using MendelImpute

makedocs(
    sitename = "MendelImpute",
    format = Documenter.HTML(),
    modules = [MendelImpute],
    pages = [
        "Home" => "index.md",
        "Phasing and Imputation" => "man/Phasing+and+Imputation.md",
        "Performance Gotchas" => "man/performance.md",
        "Estimating ancestry" => "man/painting.md",
        "Ultra compression" => "man/ultra+compress.md",
        "API" => "man/api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo   = "github.com/biona001/MendelImpute.git",
    target = "build"
)


text:
function loopback_arr(; kwargs...)
    return (; bigarrs = (;kwargs...))
end

fns = (; loopback_arr)


text:
#####
##### Intervals -> AlignedSpan
#####

is_start_exclusive(::Interval{T,L,R}) where {T,L,R} = L == Open
is_stop_exclusive(::Interval{T,L,R}) where {T,L,R} = R == Open

# Interface methods:
duration(interval::Interval{<:TimePeriod}) = last(interval) - first(interval)

function start_index_from_time(sample_rate, interval::Interval,
                               mode::Union{RoundingMode{:Up},RoundingMode{:Down}})
    first_index, error = index_and_error_from_time(sample_rate, first(interval), mode)
    if is_start_exclusive(interval) && mode == RoundUp && iszero(error)
        first_index += 1
    end
    return first_index
end

function stop_index_from_time(sample_rate, interval::Interval,
                              mode::Union{RoundingMode{:Up},RoundingMode{:Down}})
    last_index, error = index_and_error_from_time(sample_rate, last(interval), mode)
    if is_stop_exclusive(interval) && mode == RoundDown && iszero(error)
        last_index -= 1
    end
    return last_index
end

#####
##### Onda
#####

function Onda.column_arguments(samples::Samples, span::AlignedSpan)
    span.sample_rate == samples.info.sample_rate ||
        throw(ArgumentError("Sample rate of `samples` ($(samples.info.sample_rate)) does not match sample rate of `AlignedSpan` argument ($(span.sample_rate))."))
    return indices(span)
end

#####
##### TimeSpans <--> AlignedSpan
#####

TimeSpans.istimespan(::AlignedSpan) = true
TimeSpans.start(span::AlignedSpan) = time_from_index(span.sample_rate, span.first_index)
TimeSpans.stop(span::AlignedSpan) = time_from_index(span.sample_rate, span.last_index + 1)

# TimeSpan -> AlignedSpan is supported by passing to Intervals
to_interval(span) = Interval{Nanosecond,Closed,Open}(start(span), stop(span))
to_interval(span::Interval) = span
to_interval(span::AlignedSpan) = Interval(span)

# Interface methods:

function start_index_from_time(sample_rate, span, mode)
    return start_index_from_time(sample_rate, to_interval(span), mode)
end

function stop_index_from_time(sample_rate, span, mode)
    return stop_index_from_time(sample_rate, to_interval(span), mode)
end

#####
##### Serialization
#####

StructTypes.StructType(::Type{AlignedSpan}) = StructTypes.Struct()

const ARROW_ALIGNED_SPAN = Symbol("AlignedSpans.AlignedSpan")
ArrowTypes.arrowname(::Type{AlignedSpan}) = ARROW_ALIGNED_SPAN
ArrowTypes.JuliaType(::Val{ARROW_ALIGNED_SPAN}) = AlignedSpan