text:
include("helpers/listing.jl")
using Colors

function main(window)
    push!(window.assets, "widgets")
    push!(window.assets, "codemirror")
    push!(window.assets, "tex")

    colors = ["violet", "blue", "green", "yellow", "orange", "red"]
    colortile(c, w, h) = fillcolor(c, size(w, h, empty))
    vbox(
        title(3, "Layout"),
        md"""Escher UI's are built using Tiles. These Tiles can be put together in
        several ways to allow for several layouts.

        We see a few examples on how to create layouts using Escher here. For the
        purposes of these examples let us assume the use following tiles which are
        filled with 6 different colors - violet, blue, green, yellow, orange, red.
        These can be replaced by any Tile created using the Content API.""",

        h1("Sizes of Tiles"),
        md"""`width` sets the width and `height` sets the height of a tile.

        To create an empty tile of 4em x 4em that is filled with red""",
        listing("""
        widetile = Escher.width(4em, empty)
        squaretile = Escher.height(4em, widetile)
        fillcolor("red", squaretile)
        """,
            fillcolor("red", Escher.height(4em, width(4em, empty)))
        ),
        md"""We could use the `size` function, which sets the width and height of a tile.

        Making our code a bit more generic and using `size`""",
        listing("""
        colortile(c, w, h) = fillcolor(c, size(w, h, empty))
        colortile("red", 4em, 4em)
        """,
            colortile("red", 4em, 4em)
        ),
        md"We now have a way of creating tiles of our required dimensions. Now lets
        look at how we can lay these tiles out!",
        h1("Padding"),
        md"""Often we require some amount of padding around our tiles. To do this
        for an individual tile, we can use `pad`. You can a also pass an array of
        sides to selectively pad them. `left`, `right`, `top`, and `bottom` are
        possible choices.
        Extending our colortile example...""",
        listing("""
        colortile(c, w, h) = fillcolor(c, size(w, h, empty))
        redtile = colortile("red", 4em, 4em)
        paddedtile = redtile |> pad(2em)
        vbox(
            vbox("Normal tile", redtile),
            vbox("Tile with padding on allsides", paddedtile),
            hbox(
                redtile |> pad([left], 2em),
                "Only left padding",
                redtile |> pad([right], 2em),
                "Only right padding",
            ),
            vbox(
                "Only bottom padding",
                redtile |> pad([bottom], 2em),
                "Only top padding",
                redtile |> pad([top], 2em),
            ),
        )
        """,
        ),
        vskip(2em),
        h1("Creating Flows"),
        md"""Putting together flows is how UIs are built in Escher. Flows are created
        by putting together tiles using `flow`. `flow` creates a flow container with
        the components and decides what axis the components will
        be laid out on, i.e, `horizontal` or `vertical`.""",
        h3(md"Horizontal Flow - `hbox`"),
        md"""To obtain a horizontal layout of tiles, use `hbox`. `hbox` creates a
        horizontal flow of all the tiles passed to it, i.e, the main axis is the
        horizontal axis. `hbox` can also take arguments as an array of tiles.""",
        listing("""
        colors = ["violet", "blue", "green", "yellow", "orange", "red"]
        hbox(
            map(
                c -> colortile(c, 4em, 4em),
                colors
            )
        )
        """,
            hbox(
                map(
                c -> colortile(c, 4em, 4em),
                colors
                )
            )
        ),
        h3(md"Vertical Flow - `vbox`"),
        md"""To obtain a vertical layout of tiles, use `vbox`. `vbox` creates a
        vertical flow of all the tiles passed to it, i.e, the main axis is the
        vertical axis. `vbox` can also take arguments as an array of tiles.""",
        listing("""
        colors = ["violet", "blue", "green", "yellow", "orange", "red"]
        vbox(
            map(
                c -> colortile(c, 4em, 4em),
                colors
            )
        )
        """,
            vbox(
                map(
                c -> colortile(c, 4em, 4em),
                colors
                )
            )
        ),
        h3(md"Combine `vbox` and `hbox` to make complex layouts!"),
        listing("""
        colortile(c, w, h) = fillcolor(c, size(w, h, empty))
        colors = ["violet", "blue", "green", "yellow", "orange", "red"]
        colortiles = map(
            c -> colortile(c, 4em, 4em),
            colors
        )
        vcolorset = vbox(colortiles)
        hcolorset = hbox(colortiles[2:end])
        hbox(
            vcolorset,
            hcolorset
        )
        """
        ),
        h3(md"Growing tiles"),
        md"""
        Tiles can be grown to fill remaining space existing in a flow container.
        We use `grow` to create tiles that can be grown and specify a relative factor
        to decide on their rate of growth in comparison to other tiles that can
        grow. Tiles are grown only if the elements do not occupy the complete
        space along the main axis of the container.
        """,
        listing("""
        colors = ["violet", "blue", "green", "yellow", "orange", "red"]
        colortiles = map(
            c -> colortile(c, 4em, 4em),
            colors
        )
        hbox(
            colortiles[1:3]...,
            grow(2, colortiles[4]),
            colortiles[5],
            grow(1, colortiles[end]),
        )
        """
        ),
        h3(md"Shrinking tiles"),
        md"""
        Tiles can be shrinked to fit into the space existing in a flow container.
        We use `shrink` to create tiles that can be shrunk and specify a relative factor
        to decide on the rate at which they are shrunk in comparison to other tiles that can
        be shrunk. Tiles are shrunk only if the elements occupy more than the complete
        space along the main axis of the container.
        """,
        listing("""
        colors = ["violet", "blue", "green", "yellow", "orange", "red"]
        colortiles = map(
            c -> colortile(c, 8em, 4em),
            colors
        )
        hbox(
            colortiles[1:3]...,
            shrink(2, colortiles[4]),
            colortiles[5],
            shrink(1.5, colortiles[end])
        )
        """
        ),
        h3(md"Flexing tiles"),
        md"""To create a tile that will either shrink/grow to fill the remaining
        space in a flow container, use `flex`.""",
        listing("""
        colors = ["violet", "blue", "green", "yellow", "orange", "red"]
        colortiles = map(
            c -> colortile(c, 4em, 4em),
            colors
        )
        hbox(
            colortiles[1:3]...,
            flex(colortiles[4]),
            colortiles[5:end]...,
        )
        """
        ),
        h1(md"Packing spaces"),
        md"""
        While creating flows, we have several ways of packing the contents
        so as to manage space. They are :

        - `axisstart`
        - `axisend`
        - `center`
        - `spacebetween`
        - `spacearound`
        - `stretch`
        - `baseline`

        """,
        h3(md"Packing items along the main axis"),
        md"""
        We can use `packitems` to specify what part of the main axis we want to
        pack our items across. `baseline` and `stretch` cannot be used here.
        """,
        listing("""
        packings = [axisstart, axisend, center, spacebetween, spacearound]
        vbox(
            map(
                p -> vbox(
                string(typeof(p)),
                hbox(colortiles) |> packitems(p)
                ),
                packings
            )
        )
        """
        ),
        h3(md"Packing items along the cross axis"),
        md"""
        We can use `packacross` to specify what part of the cross axis we want to
        pack our items across. `spacearound` and `spacebetween` cannot be used here.
        """,
        listing("""
        packings = [axisstart, axisend, center, stretch]
        vbox(
            map(
                p -> vbox(
                string(typeof(p)),
                vbox(colortiles) |> packacross(p)
                ),
                packings
            )
        )
        """),
        h3(md"Packing items as lines"),
        md"""
        Pack wrapped lines of tiles across the cross axis.
        """,
        listing("""
        packings = [axisstart, axisend, center, spacearound, spacebetween, stretch]
        vbox(
            map(
                p -> vbox(
                string(typeof(p)),
                hbox(colortiles) |> packlines(p)
                ),
                packings
            )
        )
        """
        )
        ) |> pad(2em)
end


text:
include("applications.jl")

function train(nsite;depth::Int=2)
    symmetry = :twoqubit
    model = TFI(nsite; h=0.5, periodic=false)
    ansatz = simple_ansatz(nsite, symmetry, depth; load_params=false)

    run_train(ansatz, model; SAVE_ID=Symbol(symmetry,:_d,depth), niter=500, start_point=0)
end

function measure(task::String; symmetry::Symbol=:su2, depth::Int=5)
    model = simple_model_heis(6)
    ansatz = simple_ansatz(6, symmetry, depth; load_params=true)
    nbit = nbit_simulated(ansatz)

    if task == "energy"
        eng_sample = energy(ansatz, model)/nbit
        println("Sampled value of energy/site = $eng_sample")
    elseif task == "energy_exact"
        eng = energy_exact(ansatz, model)/nbit
        println("Exact value of energy/site = $eng")
    elseif task == "fidelity"
        gs = ground_state(model)
        fid = fidelity_exact(ansatz, gs)[]
        println("Exact value of fidelity = $fid")
    elseif task == "szsz"
        correlation_matrix(ansatz; SAVE_ID=symmetry)
    else
        throw(ArgumentError("$task is not defined!"))
    end
end

ext: jl

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
# build PyCall with the python on path
using Pkg
ENV["PYTHON"] = "python"
#ENV["PYTHON"] = raw"C:\Users\Yi-Kun.Yang\local\conda3\envs\qt\python.exe"
Pkg.build("PyCall")


text:
"""
    modelX_gaussian_knockoffs(X::Matrix{T})

Creates model-free multivariate normal knockoffs by sequentially sampling from
conditional multivariate normal distributions.

# Inputs
+ `X`: A `n × p` numeric matrix. Each row is a sample, and each column is standardized
to mean 0 variance 1.

# Reference:
"Panning for Gold: Model-X Knockoffs for High-dimensional Controlled
Variable Selection" by Candes, Fan, Janson, and Lv (2018)
"""
function modelX_gaussian_knockoffs(X::Matrix{T}, method::Symbol, μ::Vector{T}) where T <: AbstractFloat
    n, p = size(X)
    full_svd = n > p ? true : false
    # compute gram matrix using full svd
    U, σ, V = svd(X, full=full_svd)
    Σ = V * Diagonal(σ)^2 * V'
    Σinv = V * inv(Diagonal(σ)^2) * V'
    # compute s vector using the specified method
    if method == :equi
        λmin = typemax(T)
        for σi in σ
            σi^2 < λmin && (λmin = σi^2)
        end
        s = min(1, 2λmin) .* ones(size(Σ, 1))
    elseif method == :sdp
        svar = Variable(p)
        problem = maximize(sum(svar), svar ≥ 0, 1 ≥ svar, 2Σ - Diagonal(svar) in :SDP)
        solve!(problem, () -> SCS.Optimizer(verbose=false))
        s = clamp.(evaluate(svar), 0, 1) # for numeric stability
    elseif method==:asdp
        # todo
        error("ASDP not supported yet! sorry!")
    else
        error("modelX_gaussian: method can only be :equi, or :sdp, or :asdp")
    end
    X̃ = condition(X, μ, Σinv, Diagonal(s))
    return Knockoff(X, X̃, s, Σ, Σinv)
end

"""
    condition(x::AbstractVector, μ::AbstractVector, Σinv::AbstractMatrix, D::AbstractMatrix)

Samples a knockoff x̃ from Gaussian x using conditional distribution formulas:

If (x, x̃) ~ N((μ, μ), G) where G = [Σ  Σ - D; Σ - D  Σ], then we sample x̃ from
x̃|x = N(x - D*inv(Σ)(x - μ), 2D - D*inv(Σ)*D)

# todo: efficiency
"""
function condition(X::AbstractMatrix, μ::AbstractVector, Σinv::AbstractMatrix, D::AbstractMatrix)
    n, p = size(X)
    ΣinvD = Σinv * D
    new_V = Symmetric(2D - D * ΣinvD + 0.00001I) # small perturbation ensures positive eigvals
    L = cholesky(new_V).L
    return X - (X .- μ') * ΣinvD + randn(n, p) * L
end


text:
include("euler/euler.jl")
using .Primes: prime_sieve

function compute()::Int
    sieve = prime_sieve(1_000_000)
    primes = findall(sieve)[5:end]
    truncatable, result = 0, 0
    for prime in primes
        p, is_truncatable = 10, true
        while p < prime
            if sieve[prime % p] && sieve[prime ÷ p]
                p *= 10
            else
                is_truncatable = false
                break
            end
        end
        if is_truncatable
            truncatable += 1
            result += prime
            if truncatable == 11
                break
            end
        end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
include("gage_array_of_arrays.jl")
include("gage_array.jl")


text:
using Random, Statistics
N = 500
ts = cumsum(randn(N))
ts_nan = cumsum(randn(N))
ts_nan[1] = NaN
x = cos.(range(0, 20π, length = N)) .+ randn(N)*0.05
t = (0:N-1) + rand(N)

all_conceivable_methods = [
    PartialRandomization(0.3)
    PartialRandomization(0.8)
    WLS(rescale = false)
    WLS(AAFT(), rescale = true)
    WLS(rescale = false)
    WLS(CircShift(N), f = nothing)
    WLS(BlockShuffle(10), f = Statistics.cor)
    RandomCascade()
    PseudoPeriodic(3, 25, 0.05)
    BlockShuffle()
    BlockShuffle(4)
    RandomShuffle()
    AutoRegressive(1)
    AAFT()
    IAAFT()
    TFTS(0.05)
    TFTS(-0.05)
    TAAFT(0.05)
    TAAFT(-0.05)
    RandomFourier(true)
    RandomFourier(false)
    TFTD()
    TFTD(0.05)
    TFTDAAFT(0.03)
    TFTDIAAFT(0.03)
    CycleShuffle()
    IrregularLombScargle(t; tol = 10, n_total = 20000, n_acc = 5000)
]

methodnames = [string(nameof(typeof(x))) for x in all_conceivable_methods]

@testset "Reproducibility (rng)" begin
    @testset "$n" for (i, n) in enumerate(methodnames)
        method = all_conceivable_methods[i]
        rng = Random.MersenneTwister(1234)
        y = surrogate(x, method, rng)
        rng = Random.MersenneTwister(1234)
        z = surrogate(x, method, rng)
        @test y == z
    end
end


text:
function create_final_tp_df(pdir, wobs, incls; ofu=true)
  df = DataFrame()
  for wob in wobs
    for incl in incls
      filepath = joinpath(pdir, "tp", string(Int(wob)),
        string(Int(incl)), "tp.csv")
      tp_df = DataFrame(CSV.read(filepath))
      if size(df, 1) == 0
        df = DataFrame(tp_df[end, :])
      else
        append!(df, DataFrame(tp_df[end, :]))
      end
    end
  end
  df
end

function show_tp(pdir, wobs, incls; ofu=false)
  println()
  df = DataFrame()
  tp_df = DataFrame()
  for wob in wobs
    for incl in incls
      wobstring = string(Int(wob))
      inclstring = string(Int(incl))
      filepath = joinpath(pdir, "tp", wobstring, inclstring, "tp.csv")
      if ofu
        #tp_df = convert_tp_to_ofu_df(filepath)
      else
        tp_df = DataFrame(CSV.read(filepath))
      end
      if size(df, 1) == 0
        df = tp_df
      else
        append!(df, tp_df)
      end
    end
  end
  df
end

function show_solution(pdir, wobs, incls; tp=true, show=true)
  df = DataFrame()
  tp_df = DataFrame()
  for wob in wobs
    for incl in incls
      wobstring = string(Int(wob))
      inclstring = string(Int(incl))
      assembly = tp ? "tp" : "tp0"
      solpath = joinpath(pdir, assembly, wobstring, inclstring, "solution.csv")
      tppath = joinpath(pdir, assembly, wobstring, inclstring, assembly*".csv")
      tp && (tp_df = CSV.read(tppath))
      df = CSV.read(solpath)
      if show
        println("\n=====================================================")
        println("      $(assembly) (wob=$wobstring, incl=$inclstring)")
        println("=====================================================")
        println()
        size(tp_df, 1) > 0 && tp_df |> display
        size(tp_df, 1) > 0 && println()
        df |> display
        println()
      end
    end
  end
  (df, tp_df)
end

function show_solution(pdir, wob::Int, incl::Int; tp=true, show=true)
  df = DataFrame()
  tp_df = DataFrame()
  wobstring = string(wob)
  inclstring = string(incl)
  assembly = tp ? "tp" : "tp0"
  solpath = joinpath(pdir, assembly, wobstring, inclstring, "solution.csv")
  tppath = joinpath(pdir, assembly, wobstring, inclstring, assembly*".csv")
  #tp && (tp_df = convert_tp_to_ofu_df(tppath))
  #(df = convert_to_ofu_df(solpath, wobstring, inclstring))
  tp && (tp_df = CSV.read(tppath))
  df = CSV.read(solpath)
  if show
    println("\n=====================================================")
    println("      $(assembly) (wob=$wobstring, incl=$inclstring)")
    println("=====================================================")
    println()
    size(tp_df, 1) > 0 && tp_df |> display
    size(tp_df, 1) > 0 && println()
    df |> display
    println()
  end
  (df, tp_df)
end


text:
abstract type AbstractLabeledArray{T,N,AXES,TN} <: AbstractArray{T,N} end
abstract type AbstractDictArray{K,N,VS,SV} <: AbstractArray{LDict{K,SV},N} end


text:
##==============================================================================
## Common Accessors
##==============================================================================
# Common get and set methods

# NOTE this could be reduced with macros and function generation to even less code.

# Data levels
const DataLevel0 = Union{VariableDataLevel0, FactorDataLevel0}
const DataLevel1 = Union{VariableDataLevel1, FactorDataLevel1}
const DataLevel2 = Union{VariableDataLevel2, FactorDataLevel2}

##------------------------------------------------------------------------------
## label
##------------------------------------------------------------------------------

"""
$SIGNATURES

Return the label for a DFGNode.
"""
getLabel(v::DataLevel0) = v.label

##------------------------------------------------------------------------------
## tags
##------------------------------------------------------------------------------

"""
$SIGNATURES

Return the tags for a DFGNode.
"""
getTags(v::DataLevel0) = v.tags

"""
$SIGNATURES

Set the tags for a DFGNode.
"""
function setTags!(f::DataLevel0, tags::Union{Vector{Symbol},Set{Symbol}})
  empty!(f.tags)
  union!(f.tags, tags)
end

##------------------------------------------------------------------------------
## timestamp
##------------------------------------------------------------------------------

"""
$SIGNATURES

Get the timestamp of a DFGNode.
"""
getTimestamp(v::DataLevel1) = v.timestamp

"""
    $SIGNATURES

Set the timestamp of a Variable/Factor object in a factor graph.
Note:
Since `timestamp` is not mutable `setTimestamp!` calls `updateVariable!` internally.
See also [`setTimestamp`](@ref)
"""
function setTimestamp!(dfg::AbstractDFG, lbl::Symbol, ts::DateTime)
    if isVariable(dfg, lbl)
        return updateVariable!(dfg, setTimestamp(getVariable(dfg,lbl), ts))
    else
        return updateFactor!(dfg, setTimestamp(getFactor(dfg,lbl), ts))
    end
end

##------------------------------------------------------------------------------
## solvable
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Variables or factors may or may not be 'solvable', depending on a user definition.  Useful for ensuring atomic transactions.

Related:
- isSolveInProgress
"""
getSolvable(var::Union{DFGVariable, DFGFactor})::Int = var.solvable
#TODO DataLevel2

"""
    $SIGNATURES

Get 'solvable' parameter for either a variable or factor.
"""
function getSolvable(dfg::AbstractDFG, sym::Symbol)::Int
  if isVariable(dfg, sym)
    return getVariable(dfg, sym).solvable
  elseif isFactor(dfg, sym)
    return getFactor(dfg, sym).solvable
  end
end

#TODO data level2 for N
"""
    $SIGNATURES

Set the `solvable` parameter for either a variable or factor.
"""
function setSolvable!(node::N, solvable::Int)::Int where N <: DFGNode
  node.solvable = solvable
  return solvable
end

"""
    $SIGNATURES

Set the `solvable` parameter for either a variable or factor.
"""
function setSolvable!(dfg::AbstractDFG, sym::Symbol, solvable::Int)::Int
  if isVariable(dfg, sym)
    getVariable(dfg, sym).solvable = solvable
  elseif isFactor(dfg, sym)
    getFactor(dfg, sym).solvable = solvable
  end
  return solvable
end

"""
    $SIGNATURES

Variables or factors may or may not be 'solvable', depending on a user definition.
returns true if `getSolvable` > 0
Related:
- `getSolvable`(@ref)
"""
isSolvable(node::Union{DFGVariable, DFGFactor}) = getSolvable(node) > 0

##------------------------------------------------------------------------------
## solveInProgress
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Which variables or factors are currently being used by an active solver.  Useful for ensuring atomic transactions.

DevNotes:
- Will be renamed to `data.solveinprogress` which will be in VND, not DFGNode -- see DFG #201

Related

isSolvable
"""
function getSolveInProgress(var::Union{DFGVariable, DFGFactor}, solveKey::Symbol=:default)::Int
    # Variable
    var isa DFGVariable && return haskey(getSolverDataDict(var), solveKey) ? getSolverDataDict(var)[solveKey].solveInProgress : 0
    # Factor
    return getSolverData(var).solveInProgress
end

#TODO missing set solveInProgress and graph level accessor

isSolveInProgress(node::Union{DFGVariable, DFGFactor}, solvekey::Symbol=:default) = getSolveInProgress(node, solvekey) > 0

##==============================================================================
## Common Layer 2 CRUD and SET
##==============================================================================

##==============================================================================
## TAGS as a set, list, merge, remove, empty
##==============================================================================
"""
$SIGNATURES

Return the tags for a variable or factor.
"""
function listTags(dfg::AbstractDFG, sym::Symbol)
  getFnc = isVariable(dfg,sym) ? getVariable : getFactor
  getTags(getFnc(dfg, sym))
end
#alias for completeness
listTags(f::DataLevel0) = getTags(f)

"""
    $SIGNATURES

Merge add tags to a variable or factor (union)
"""
function mergeTags!(dfg::InMemoryDFGTypes, sym::Symbol, tags::Vector{Symbol})
  getFnc = isVariable(dfg,sym) ? getVariable : getFactor
  union!(getTags(getFnc(dfg, sym)), tags)
end
mergeTags!(f::DataLevel0, tags::Vector{Symbol}) = union!(f.tags, tags)

"""
$SIGNATURES

Remove the tags from the node (setdiff)
"""
function removeTags!(dfg::InMemoryDFGTypes, sym::Symbol, tags::Vector{Symbol})
  getFnc = isVariable(dfg,sym) ? getVariable : getFactor
  setdiff!(getTags(getFnc(dfg, sym)), tags)
end
removeTags!(f::DataLevel0, tags::Vector{Symbol}) = setdiff!(f.tags, tags)

"""
$SIGNATURES

Empty all tags from the node (empty)
"""
function emptyTags!(dfg::InMemoryDFGTypes, sym::Symbol)
  getFnc = isVariable(dfg,sym) ? getVariable : getFactor
  empty!(getTags(getFnc(dfg, sym)))
end
emptyTags!(f::DataLevel0) = empty!(f.tags)


text:
#################################################
# Two Step Reconstruction			#
#################################################
# Example for the paper
# "Simultaneous imaging of widely differing
#  particle concentrations in MPI:
#  problem statement and algorithmic proposal
#  for improvement"
# DOI: 10.1088/1361-6560/abf202
#################################################
# Phantom: Vessel and Kidneys (Dilution: 1 : 2^(-5))
# Julia: 1.6

##############################
# Activate local environment #
##############################
using Pkg
Pkg.activate(".")
Pkg.instantiate()

# load packages
using MPIReco
using Plots, Plots.PlotMeasures
using LazyArtifacts

pyplot() # use PyPlot backend for plotting

###################
# Parameter setup #
###################
## Regularization parameter
λ_high = 0.001
λ_low = 0.5
# SNR threshold
Θ_high = 2.0
Θ_low = 15
# iterations
ι_high = 3
ι_low = 1

# Threshold
Γ = 0.05

# Frequency selection
minFreq = 80e3
recChannels = [1,2,3]

# Plotting
vmin = 0.0
vmax = 0.0094
MIPplanes = 10

########
# Data #
########
@info "Load data"
# get calibration and experiment data
pathToMDFStore = artifact"MDFStore"
store = MDFDatasetStore(pathToMDFStore)
study = getStudies(store,"IncreasedDynamicRange")[1]
exp = getExperiments(study)

# load MPIFiles
bEmpty = MPIFile(exp[2]) # background scan
bMeas = MPIFile(exp[1]) # measurement
bSF = MPIFile(joinpath(calibdir(store),"1.mdf")) # system matrix

# Frequency selection
freq_high = filterFrequencies(bSF, minFreq=minFreq, SNRThresh=Θ_high, recChannels=recChannels)
freq_low = filterFrequencies(bSF, minFreq=minFreq, SNRThresh=Θ_low, recChannels=recChannels)

# load measurements
uEmpty_high = getMeasurementsFD(bEmpty,frequencies=freq_high, numAverages=acqNumFrames(bEmpty)) # background
uEmpty_low = getMeasurementsFD(bEmpty,frequencies=freq_low, numAverages=acqNumFrames(bEmpty)) # background
uMeas_high = getMeasurementsFD(bMeas,frequencies=freq_high,
				numAverages=acqNumFrames(bMeas), spectralLeakageCorrection=false) # measurement
uMeas_low = getMeasurementsFD(bMeas,frequencies=freq_low,
				numAverages=acqNumFrames(bMeas), spectralLeakageCorrection=false) # measurement
uMeas_high .-= uEmpty_high # subtract background
uMeas_low .-= uEmpty_low # subtract background

# load system matrix
S_high, grid = getSF(bSF, freq_high, nothing, "kaczmarz"; bgcorrection=true)
S_low, grid = getSF(bSF, freq_low, nothing, "kaczmarz"; bgcorrection=true)

###########################
# Two-Step Reconstruction #
###########################
@info "Start two-step reconstruction"

## Step 1: First reconstruction
cPre = reconstruction(S_high, uMeas_high, λ=λ_high, iterations=ι_high)

## Step 2: Thresholding
cThresh = copy(cPre)
cThresh[ abs.(cPre).< maximum(abs.(cPre))*Γ ] .= 0

## Step 3: Projection into raw data space
uProj = zero(uMeas_low)
uProj[:,1,1] = map(ComplexF32,S_low*vec(cThresh))

## Step 4: Subtraction
uCorr = uMeas_low - uProj

## Step 5: Second reconstruction
cPost = reconstruction(S_low, uCorr, λ=λ_low, iterations=ι_low)

## Step 6: Addition
c = cPost + cThresh

##############################################
# Reconstruction with P_low (for comparison) #
##############################################
@info "Start reconstruction with P_low"
cLow = reconstruction(S_low, uMeas_low, λ=λ_low, iterations=ι_low)

#################
# Visualization #
#################
# Reshaping
cPre = reshape(cPre[:,1], shape(grid)...)
cPost = reshape(cPost[:,1], shape(grid)...)
cLow = reshape(cLow[:,1], shape(grid)...)

# Create plots
p1 = heatmap(-21:2:21,-21:2:21,squeeze(reverse(maximum(cPre[:,:,MIPplanes],dims=3),dims=1)),
             clim = (vmin,maximum(cPre)), c = :viridis,
             colorbar_ticks = ([0,maximum(cPre)], [0,"κᵥ"]),
             title = "cPre", aspect_ratio = 1 )
p2 = heatmap(-21:2:21,-21:2:21,squeeze(reverse(maximum(cPost[:,:,MIPplanes],dims=3),dims=1)),
             clim=(vmin,vmax), c = :viridis,
             colorbar_ticks = ([0,vmax],[0,"κ₁"]),
             title = "cPost", aspect_ratio = 1 )
p3 = heatmap(-21:2:21,-21:2:21,squeeze(reverse(maximum(cLow[:,:,MIPplanes],dims=3),dims=1)),
             clim = (vmin,vmax), c = :viridis,
             colorbar_ticks = ([0,vmax],[0,"κ₁"]),
             title = "cLow", aspect_ratio = 1 )

# Figure
plot(p1, p2, p3, layout = (1, 3), size=(1000,300),
     xaxis = ("y / mm", (-21,21)),
     yaxis = ("x / mm", (-21,21)),
     colorbar_title = "concentration",
     left_margin = 5mm,
     tickfontsize = 10,
     colorbar_tickfontsize = 10)

max_stars_repo_name: IBIResearch/IncreasedDynamicRange

max_issues_repo_name: IBIResearch/IncreasedDynamicRange

max_forks_repo_name: IBIResearch/IncreasedDynamicRange


text:
# testing.jl

import Base: ==

==(lhs::Int, rhs::QueryResult) = lhs == rhs.info.value
==(lhs::AbstractString, rhs::QueryResult) = lhs == rhs.info.value
==(lhs::Float64, rhs::QueryResult) = lhs == rhs.info.value
==(lhs::Bool, rhs::QueryResult) = lhs == rhs.info.value
==(lhs::Vector{AbstractString}, rhs::QueryResult) = lhs == rhs.info.value
==(lhs::Dict, rhs::QueryResult) = lhs == rhs.info.value
==(lhs::Void, rhs::QueryResult) = lhs == rhs.info.value


text:
@testset "CRUD" begin
  println(testname("CRUD tests"))
  print("  [  ] Read a non-existing id ")
  @test_throws HTTP.ExceptionRequest.StatusError readdoc(db, "this-id-does-not-exist")
  println("\r  [OK] Read a non-existing id")

  print("  [  ] Create new doc ")
  data = createdoc(db, Dict("item" => "flange", "location" => "under-stairs cupboard"))
  @test haskey(data, "id")
  @test haskey(data, "rev")
  println("\r  [OK] Create new doc")

  print("  [  ] Create new doc with empty body should fail ")
  @test_throws ErrorException createdoc(db, Dict())
  println("\r  [OK] Create new doc with empty body should fail")

  print("  [  ] Read new doc by {id, rev} (note: bad idea usually) ")
  doc = readdoc(db, data["id"]; rev=data["rev"])
  @test haskey(doc, "item")
  @test doc["item"] == "flange"
  println("\r  [OK] Read new doc by {id, rev} (note: bad idea usually)")

  print("  [  ] Read doc with exotic params ")
  docs = readdoc(db, data["id"];
    open_revs=[data["rev"]],
    conflicts=true,
    attachments=true,
    atts_since=[data["rev"]],
    att_encoding_info=true,
    latest=true,
    meta=true,
    deleted_conflicts=true,
    revs=true,
    revs_info=true
  )

  doc = docs[1]["ok"]
  @test haskey(doc, "item")
  @test doc["item"] == "flange"
  println("\r  [OK] Read doc with exotic params")

  print("  [  ] Read doc with all open revs ")
  docs = readdoc(db, data["id"]; open_revs=["all"])
  doc = docs[1]["ok"]
  @test haskey(doc, "item")
  @test doc["item"] == "flange"
  println("\r  [OK] Read doc with all open revs")

  print("  [  ] Reading doc by id and bad rev should fail ")
  @test_throws HTTP.ExceptionRequest.StatusError readdoc(db, data["id"]; rev="3-63453748494907")
  println("\r  [OK] Reading doc by id and bad rev should fail")

  print("  [  ] Update existing doc ")
  doc = updatedoc(db; id=data["id"], rev=data["rev"], body=Dict("item" => "flange", "location" => "garage"))
  @test haskey(doc, "rev")
  @test occursin("2-", doc["rev"])
  println("\r  [OK] Update existing doc")

  print("  [  ] Delete doc ")
  newdoc = deletedoc(db; id=doc["id"], rev=doc["rev"])
  @test newdoc["ok"] == true
  println("\r  [OK] Delete doc")
end

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla


text:
module mynewcode
function viewLarIndices(W,EW)
	submodel = lar2hpc(W,hcat(EW...))
	VV = PyObject(Any[Any[k] for k in  0:size(W,2)-1])
	EZ = map(Array{Int32},EW-1)
	EV = PyObject(Any[PyObject(EZ[e])[:tolist]() for e=1:length(EZ) ])
	V = PyObject(Any[PyObject(W[:,v])[:tolist]() for v=1:size(W,2) ])
	size = maximum(p.SIZE(Any[1,2])(submodel))/6
	hpc = p.larModelNumbering(1,1,1)(V,PyObject([VV,EV]),submodel,size)
	p.VIEW(hpc)
end
end


text:
# NOTE: The `next_frame_id' and `frame_names' globals below are a hack, but they
# enable a significant reduction in allocations.
# Storing the names of all CartesianFrame3D objects in this frame_names vector instead
# of in the CartesianFrame3D and having CartesianFrame3D only contain an integer ID
# makes CartesianFrame3D an isbits (pointer free) type.
# This in turn makes it so that a lot of the geometry/dynamics types become isbits
# types, making them stack allocated and allowing all sorts of
# optimizations.
const next_frame_id = Ref(0)
const frame_names = Dict{Int64, String}()

"""
$(TYPEDEF)

A `CartesianFrame3D` identifies a three-dimensional Cartesian coordinate system.

`CartesianFrame3D`s are typically used to annotate the frame in which certain
quantities are expressed.
"""
struct CartesianFrame3D
    id::Int64

    @doc """
    $(SIGNATURES)

    Create a `CartesianFrame3D` with the given name.
    """ ->
    function CartesianFrame3D(name::String)
        ret = new(next_frame_id.x)
        next_frame_id.x = Base.Checked.checked_add(next_frame_id.x, 1)
        frame_names[ret.id] = name
        ret
    end

    @doc """
    $(SIGNATURES)

    Create an anonymous `CartesianFrame3D`.
    """ ->
    function CartesianFrame3D()
        ret = new(next_frame_id.x)
        next_frame_id.x = Base.Checked.checked_add(next_frame_id.x, 1)
        ret
    end
end

Base.print(io::IO, frame::CartesianFrame3D) = print(io, get(frame_names, frame.id, "anonymous"))
name_and_id(frame::CartesianFrame3D) = "\"$frame\" (id = $(frame.id))"
Base.show(io::IO, frame::CartesianFrame3D) = print(io, "CartesianFrame3D: $(name_and_id(frame))")

Base.in(x::CartesianFrame3D, y::CartesianFrame3D) = x == y

"""
$(SIGNATURES)

Check that `CartesianFrame3D` `f1` is one of `f2s`.

Note that if `f2s` is a `CartesianFrame3D`, then `f1` and `f2s` are simply checked for equality.

Throws an `ArgumentError` if `f1` is not among `f2s` when bounds checks
are enabled. `@framecheck` is a no-op when bounds checks are disabled.
"""
macro framecheck(f1, f2s)
    quote
        @boundscheck begin
            $(esc(f1)) ∉ $(esc(f2s)) && framecheck_fail($(QuoteNode(f1)), $(QuoteNode(f2s)), $(esc(f1)), $(esc(f2s)))
        end
    end
end

framecheck_string(expr, frame::CartesianFrame3D) = "$expr: $(name_and_id(frame))"

@noinline function framecheck_fail(expr1, expr2, f1::CartesianFrame3D, f2::CartesianFrame3D)
    throw(ArgumentError("$(framecheck_string(expr1, f1)) ≠ $(framecheck_string(expr2, f2))"))
end

@noinline function framecheck_fail(expr1, expr2, f1::CartesianFrame3D, f2s)
    buf = IOBuffer()
    print(buf, '(')
    first = true
    for f2 in f2s
        first || print(buf, ", ")
        first = false
        print(buf, name_and_id(f2))
    end
    print(buf, ')')
    throw(ArgumentError("$(framecheck_string(expr1, f1)) ∉ $(expr2): $(String(take!(buf)))"))
end


text:
function sdc(p::AbstractNFFTPlan{T,D,1}; iters=20) where {T,D}
    # Weights for sample density compensation.
    # Uses method of Pipe & Menon, 1999. Mag Reson Med, 186, 179.
    weights = ones(Complex{T}, p.M)
    weights_tmp = similar(weights)

    pre = p.params.precompute
    if pre == LINEAR
      scaling_factor = maximum(p.windowLinInterp)^2
    elseif pre == POLYNOMIAL
      scaling_factor = maximum(p.windowPolyInterp[1,:])^2
    elseif pre == TENSOR
      scaling_factor = maximum(maximum(p.windowTensor))^2
    elseif pre == FULL
      scaling_factor = maximum(p.B)^2
    end

    # Pre-weighting to correct non-uniform sample density
    for i in 1:iters
        convolve_transpose!(p, weights, p.tmpVec)
        p.tmpVec ./= scaling_factor
        convolve!(p, p.tmpVec, weights_tmp)
        weights_tmp ./= scaling_factor
        for j in 1:length(weights)
            weights[j] = weights[j] / (abs(weights_tmp[j]) + eps(T))
        end
    end
    # Post weights to correct image scaling
    # This finds c, where ||u - c*v||_2^2 = 0 and then uses
    # c to scale all weights by a scalar factor.
    u = ones(Complex{T}, p.N)
    f = p * u
    f = f .* weights # apply weights from above
    v = adjoint(p) * f
    c = vec(v) \ vec(u)  # least squares diff
    abs.(c * weights)
end


text:
struct MetricTerms
    ∂x
    ∂y
    ∂z
end

struct VolumeGeometry
    coordinates
    jacobian
    metricterms
    mass
    inversemass
end

struct LocalCommunication
    minus
    plus
end

struct MPICommunication
    interior
    exterior
end

struct SurfaceGeometry
    coordinates
    mass
    inversemass
end

struct ElementOperators
    mass
    inversemass
    DifferentiationMatrix
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
include("DeterministicInternal.jl")
include("Deterministic.jl")
include("ProbabilisticInternal.jl")
include("Probabilistic.jl")
include("ScenarioBasedInternal.jl")
include("ScenarioBased.jl")

export get_data
export get_horizon
export get_initial_time
export get_internal
export get_label
export get_percentiles
export get_resolution
export get_scenario_count
export get_time_series_uuid


text:
"""
    GibbsSampling(;ϵ::T=1e-5, nBurnin::Int=100, thinning::Int=1)

Draw samples from the true posterior via Gibbs Sampling.

## Keywords arguments
- `ϵ::T` : convergence criteria
- `nBurnin::Int` : Number of samples discarded before starting to save samples
- `thinning::Int` : Frequency at which samples are saved
"""
mutable struct GibbsSampling{T<:Real} <: SamplingInference{T}
    nBurnin::Int # Number of burnin samples
    thinning::Int # Frequency at which samples are saved
    ϵ::T # Convergence criteria
    n_iter::Integer # Number of samples computed
    HyperParametersUpdated::Bool # Flag for updating kernel matrices
    opt::SOptimizer
    sample_store::Vector{Vector{Vector{T}}}
    function GibbsSampling{T}(nBurnin::Int, thinning::Int, ϵ::Real) where {T}
        nBurnin >= 0 || error("nBurnin should be positive")
        thinning >= 0 || error("thinning should be positive")
        return new{T}(
            nBurnin,
            thinning,
            ϵ,
            0,
            false,
            SOptimizer(Descent()),
            Vector{Vector{Vector{T}}}[],
        )
    end
end

ρ(::GibbsSampling{T}) where {T} = one(T)

function GibbsSampling(; ϵ::T=1e-5, nBurnin::Int=100, thinning::Int=1) where {T<:Real}
    return GibbsSampling{T}(nBurnin, thinning, ϵ)
end

function Base.show(io::IO, ::GibbsSampling{T}) where {T<:Real}
    return print(io, "Gibbs Sampler")
end

function sample_local!(
    local_vars, l::AbstractLikelihood, y, f::Tuple{<:AbstractVector{T}}
) where {T}
    return sample_local!(local_vars, l, y, first(f))
end

function sample_global!(
    ∇E_μ::AbstractVector,
    ∇E_Σ::AbstractVector,
    X::AbstractVector,
    gp::SampledLatent{T},
    kernel_mat,
) where {T}
    gp.post.Σ .= inv(Symmetric(2.0 * Diagonal(∇E_Σ) + inv(kernel_mat.K)))
    rand!(MvNormal(cov(gp) * (∇E_μ + kernel_mat.K \ pr_mean(gp, X)), cov(gp)), gp.post.f)
    return copy(posterior(gp).f)
end


text:
# Functions for creating batch files to run on NIH Biowulf super computer

using DelimitedFiles

# transientgenelist = ["MYC", "DUSP5", "TRIB1", "PMAIP1", "SERPINE1", "SOX9", "ERRFI1", "NR4A2", "JUN", "GEM"]

"""
makeswarm(;G::Int=2,cell="HCT116",infolder="infolder",swarmfile::String="swarmfile",inlabel="inlabel",label="label",nsets=2,datafolder::String="datafolder",thresholdlow::Float64=0.,thresholdhigh::Float64=100000000.,conds::String="DMSO-AUXIN",result::String= "2021-03-11",batchsize=1000,maxtime = 3600. * 2,nchains::Int = 2,runcycle::Bool=true,transient::Bool=false,fittedparam::String="",fixedeffects::String="",juliafile::String="/home/carsonc/StochasticGene/runfiles/fitscript.jl",root="/home/carsonc/scrna/")

returns swarmfile executed on biowulf

"""
function makeswarm(;G::Int=2,cell="HCT116",infolder="infolder",swarmfile::String="swarmfile",inlabel="inlabel",label="label",nsets=2,datafolder::String="datafolder",thresholdlow::Float64=0.,thresholdhigh::Float64=100000000.,conds::String="DMSO-AUXIN",result::String= "2021-03-11",batchsize=1000,maxtime = 3600. * 2,nchains::Int = 2,runcycle::Bool=true,transient::Bool=false,fittedparam::String="",fixedeffects::String="",juliafile::String="/home/carsonc/StochasticGene/runfiles/fitscript.jl",root="/home/carsonc/scrna/")
    if occursin.("-",conds)
        cond = string.(split(conds,"-"))
    else
        cond = conds
    end
    genes = checkgenes(root,cond,datafolder,cell,thresholdlow,thresholdhigh)
    makeswarm(genes,G=G,cell=cell,infolder=infolder,swarmfile=swarmfile,inlabel=inlabel,label=label,nsets=nsets,datafolder=datafolder,conds=conds,result=result,batchsize=batchsize,maxtime=maxtime,nchains=nchains,runcycle=runcycle,transient=transient,fittedparam=fittedparam,fixedeffects=fixedeffects,juliafile=juliafile,root=root)
end
function makeswarm(genes::Vector;G::Int=2,cell="HCT116",infolder="infolder",swarmfile::String="swarmfile",inlabel="inlabel",label="label",nsets=2,datafolder::String,conds::String="DMSO-AUXIN",result::String="testout",batchsize=1000,maxtime=60.,nchains::Int=2,runcycle::Bool=true,transient::Bool=false,fittedparam::String="",fixedeffects="",juliafile::String="/home/carsonc/StochasticGene/runfiles/fitscript.jl",root="/home/carsonc/scrna/")
    resultfolder = joinpath("Results",result)
    infolder = joinpath("Results",infolder)
    ngenes = length(genes)
    println(ngenes)
    println(runcycle)
    if ngenes > batchsize
        batches = getbatches(genes,ngenes,batchsize)
        for batch in eachindex(batches)
            sfile = swarmfile * "_" * conds * "$batch" * ".swarm"
            writegenes(sfile,batches[batch],cell,nchains,juliafile,conds,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,fittedparam,fixedeffects)
        end
    else
        sfile = swarmfile * "_" * conds * ".swarm"
        f = open(sfile,"w")
        writegenes(sfile,genes,cell,nchains,juliafile,conds,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,fittedparam,fixedeffects)
    end
    make_fitfile(juliafile,fittedparam,fixedeffects)
end

"""
make_fitfile(fitfile,fittedparam,fixedeffects)

make the file the swarm file calls to execute julia code

"""
function make_fitfile(fitfile,fittedparam,fixedeffects)
        f = open(fitfile,"w")
        s =   '"'
        write(f,"@everywhere include($s/home/carsonc/StochasticGene/src/StochasticGene.jl$s)\n")
        write(f,"include($s/home/carsonc/StochasticGene/runfiles/scriptfunctions.jl$s)\n")
        if fittedparam == ""
            write(f,"@time fit_rna(parse(Int,ARGS[1]),ARGS[2],ARGS[3],ARGS[4],parse(Int,ARGS[5]),parse(Float64,ARGS[6]),ARGS[7],ARGS[8],ARGS[9],ARGS[10],ARGS[11],parse(Int,ARGS[12]),parse(Bool,ARGS[13]),parse(Bool,ARGS[14]))\n")
        elseif fixedeffects == ""
            write(f,"@time fit_rna(parse(Int,ARGS[1]),ARGS[2],ARGS[3],ARGS[15],ARGS[4],parse(Int,ARGS[5]),parse(Float64,ARGS[6]),ARGS[7],ARGS[8],ARGS[9],ARGS[10],ARGS[11],parse(Int,ARGS[12]),parse(Bool,ARGS[13]),parse(Bool,ARGS[14]))\n")
        else
            write(f,"@time fit_rna(parse(Int,ARGS[1]),ARGS[2],ARGS[3],ARGS[15],ARGS[16],ARGS[4],parse(Int,ARGS[5]),parse(Float64,ARGS[6]),ARGS[7],ARGS[8],ARGS[9],ARGS[10],ARGS[11],parse(Int,ARGS[12]),parse(Bool,ARGS[13]),parse(Bool,ARGS[14]))\n")
        end

        close(f)
end

function getbatches(genes,ngenes,batchsize)
    nbatches = div(ngenes,batchsize)
    batches = Vector{Vector{String}}(undef,nbatches+1)
    println(batchsize," ",nbatches)
    for i in 1:nbatches
        batches[i] = genes[batchsize*(i-1)+1:batchsize*(i)]
    end
    batches[end] = genes[batchsize*nbatches+1:end]
    return batches
end

function writegenes(sfile,genes,cell,nchains,juliafile,cond,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,fittedparam,fixedeffects)
    f = open(sfile,"w")
    for gene in genes
        gene = check_genename(gene)
        writedlm(f,["julia -p" nchains juliafile nchains gene cell cond G maxtime infolder resultfolder datafolder inlabel label nsets runcycle transient fittedparam fixedeffects])
    end
    close(f)
end

function check_genename(gene)
    if occursin("(",gene)
        gene = replace(gene,"(" => "[")
        gene = replace(gene,")" => "]")
    end
    return gene
end

function checkgenes(root,conds::Vector,datafolder,celltype::String,thresholdlow::Float64,thresholdhigh::Float64)
    genes = Vector{Vector}(undef,2)
    if occursin.("-",datafolder)
        datafolder = string.(split(datafolder,"-"))
        for i in 1:2
            genes[i] = checkgenes(root,conds[i],datafolder[i],celltype,thresholdlow,thresholdhigh)
        end
    else
        for i in 1:2
            genes[i] = checkgenes(root,conds[i],datafolder,celltype,thresholdlow,thresholdhigh)
        end
    end
    intersect(genes[1],genes[2])
end

function checkgenes(root,cond::String,datafolder,cell::String,thresholdlow::Float64,thresholdhigh::Float64)
    genes = intersect(get_halflives(root,cell,thresholdlow,thresholdhigh), get_genes(root,cond,datafolder))
    alleles = get_alleles(root,cell)
    if alleles != nothing
        return intersect(genes,alleles)
    else
        return genes
    end
end

function get_genes(root,cond,datafolder)
    genes = Vector{String}(undef,0)
    files = readdir(joinpath(root,datafolder))
    for file in files
        if occursin(cond,file)
            push!(genes,split(file,"_")[1])
        end
    end
    return genes
end

function get_halflives(root,cell,thresholdlow::Float64,thresholdhigh::Float64)
    file = get_file(root,"data/halflives",cell)
    get_halflives(file,thresholdlow,thresholdhigh)
end

function get_halflives(file,thresholdlow::Float64,thresholdhigh::Float64)
    genes = Vector{String}(undef,0)
    halflives = readdlm(file,',')
    for row in eachrow(halflives)
        if typeof(row[2]) <: Number
            if thresholdlow <= row[2] < thresholdhigh
                push!(genes,string(row[1]))
            end
        end
    end
    return genes
end

function get_alleles(root,cell)
    file = get_file(root,"data/alleles",cell)
    if file !=nothing
        return readdlm(file)[2:end,1]
    else
        return nothing
    end
end

function get_file(root,folder,type)
    folder = joinpath(root,folder)
    files = readdir(folder)
    for file in files
        if occursin(type,file)
            path = joinpath(folder,file)
            return path
        end
    end
    nothing
end

fix(folder) = writeruns(fixruns(findjobs(folder)))

function findjobs(folder)
    files = readdir(folder)
    files = files[occursin.("swarm_",files)]
    for (i,file) in enumerate(files)
        files[i] = split(file,"_")[2]
    end
    unique(files)
end

function fixruns(jobs,message="FAILED")
    runlist = Vector{String}(undef,0)
    for job in jobs
        if occursin(message,read(`jobhist $job`,String))
            swarmfile = findswarm(job)
            list = readdlm(swarmfile,',')
            runs =  chomp(read(pipeline(`jobhist $job`, `grep $message`),String))
            runs = split(runs,'\n')
            println(job)
            for run in runs
                linenumber = parse(Int,split(split(run," ")[1],"_")[2]) + 1
                while linenumber < length(list)
                    a = String(list[linenumber])
                    linenumber += 1000
                    println(a)
                    push!(runlist,a)
                end
            end
        end
    end
    return runlist
end

function writeruns(runs,outfile="fitfix.swarm")

    f = open(outfile,"w")
    for run in runs
        writedlm(f,[run],quotes=false)
    end
    close(f)

end

function findswarm(job)
    sc = "Swarm Command"
    line = read(pipeline(`jobhist $job`, `grep $sc`),String)
    list = split(line," ")
    list[occursin.(".swarm",list)][1]
end

function collate(folders,type="rates_scRNA_T120")
    genes = Vector{Vector}(undef,0)
    for folder in folders
        files = readdir(folder)
        for file in files
            if occursin(type,file)
                contents=readdlm(file,',')
                for row in eachrow(contents)
                    row[1] .== genes
                end
            end
        end
    end
end

large_deviance(measurefile,threshold) = filter_gene(measurefile,"Deviance",threshold)

function filter_gene(measurefile,measure,threshold)
    genes = Vector{String}(undef,0)
    measures,header = readdlm(measurefile,',',header=true)
    println(length(measures[:,1]))
    col = findfirst(header[1,:] .== measure)
    for d in eachrow(measures)
        if d[col] > threshold || isnan(d[col])
            push!(genes,d[1])
        end
    end
    println(length(genes))
    return genes
end

function filter_gene_nan(measurefile,measure)
    genes = Vector{String}(undef,0)
    measures,header = readdlm(measurefile,',',header=true)
    println(length(measures[:,1]))
    col = findfirst(header[1,:] .== measure)
    for d in eachrow(measures)
        if isnan(d[col])
            push!(genes,d[1])
        end
    end
    println(length(genes))
    return genes
end

function get_missing_genes(folder,type,label,cond,model)
    genes = checkgenes("DMSO",folder,0.,100000000.,"../")
    genes1=StochasticGene.getgenes(folder,type,label,cond,model)
    union(setdiff(genes1,genes),setdiff(genes,genes1))
end

function get_missing_genes(genes,folder,type,label,cond,model)
    genes1=StochasticGene.getgenes(folder,type,label,cond,model)
    union(setdiff(genes1,genes),setdiff(genes,genes1))
end

function scan_swarmfiles(jobid,folder=".")
    if ~(typeof(jobid) <: String)
        jobid = string(jobid)
    end
    genes = Array{String,1}(undef,0)
    files = readdir(folder)
    files = files[occursin.(jobid,files)]
    for file in files
        genes = vcat(genes,scan_swarmfile(file))
    end
    return genes
end

function scan_swarmfile(file)
    genes = Array{String,1}(undef,0)
    contents = readdlm(file,'\t')
    lines = contents[occursin.("[\"",string.(contents))]
    for line in lines
        push!(genes,split.(string(line)," ")[1])
    end
    return genes
end

function scan_fitfile(file,folder=".")
    genes = Array{String,1}(undef,0)
    joinpath(folder,file)
    file = readdlm(file,'\t')
    for line in eachrow(file)
        push!(genes,line[4])
    end
    return genes
end

function replace_yield(G,folder1,folder2,cond1,cond2,outfolder)
    if typeof(G) <: Number
        G = string(G)
    end
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    files1 = getratefile(folder1,G,cond1)
    files2 = getratefile(folder2,G,cond2)
    for file1 in files1
        gene = StochasticGene.getgene(file1)
        file2 = getratefile(files2,gene)
        outfile = joinpath(outfolder,file2)
        r1 = StochasticGene.readrates(joinpath(folder1,file1))
        r2 = StochasticGene.readrates(joinpath(folder2,file2))
        r2[end] = r1[end]
        f = open(outfile,"w")
        writedlm(f,[r2],',')
        close(f)
    end

end

"""
assemble_r(G,folder1,folder2,cond1,cond2,outfolder)

Combine rates from two separate fits into a single rate vector

"""

function assemble_r(G,folder1,folder2,cond1,cond2,outfolder)
    if typeof(G) <: Number
        G = string(G)
    end
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    files1 = getratefile(folder1,G,cond1)
    files2 = getratefile(folder2,G,cond2)
    for file1 in files1
        gene = StochasticGene.getgene(file1)
        file2 = getratefile(files2,gene)
        if file2 != 0
            file2=joinpath(folder2,file2)
        else
            file2=joinpath(folder1,file1)
        end
        name = replace(file1, cond1 => cond1 * "-" * cond2)
        outfile = joinpath(outfolder,name)
        assemble_r(joinpath(folder1,file1),file2,outfile)
    end
end

function  assemble_r(ratefile1,ratefile2,outfile)
    r1 = StochasticGene.readrates(ratefile1,2)
    r2 = StochasticGene.readrates(ratefile2,2)
    r1[end] = clamp(r1[end],eps(Float64),1-eps(Float64))
    r = vcat(r1[1:end-1],r2[1:end-1],r1[end])
    f = open(outfile,"w")
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    close(f)
end

function assemble_r(gene,G,folder1,folder2,cond1,cond2,outfolder)
    file1 = getratefile(gene,G,folder1,cond1)[1]
    file2 = getratefile(gene,G,folder2,cond2)[1]
    name = replace(file1, cond1 => cond1 * "-" * cond2)
    println(name)
    outfile = joinpath(outfolder,name)
    println(outfile)
    assemble_r(joinpath(folder1,file1),joinpath(folder2,file2),outfile)
end

function getratefile(files,gene)
    files = files[occursin.("_"*gene*"_",files)]
    if length(files) > 0
        return files[1]
    else
        # println(gene)
        return 0
    end
end

function getratefile(folder,G,cond)
    files = readdir(folder)
    files = files[occursin.("rates_",files)]
    files = files[occursin.("_"*cond*"_",files)]
    files[occursin.("_"*G*"_",files)]
end

getratefile(gene,G,folder,cond) = getfile("rate",gene,G,folder,cond)

function getfile(type,gene::String,G::String,folder,cond)
    files = readdir(folder)
    files = files[occursin.(type,files)]
    files = files[occursin.("_"*gene*"_",files)]
    files = files[occursin.("_"*G*"_",files)]
    files[occursin.("_"*cond*"_",files)]
end

function change_name(folder,oldname,newname)
    files = readdir(folder)
    files = files[occursin.(oldname,files)]
    for file in files
        newfile = replace(file, oldname => newname)
        mv(joinpath(folder,file),joinpath(folder,newfile),force=true)
    end
end

function make_halflife(infile,outfile,col=4)
    f = open(outfile,"w")
    writedlm(f,["Gene" "Halflife"],',')
    contents,rows = readdlm(infile,',',header=true)
    for row = eachrow(contents)
        gene = string(row[1])
        gene = strip(gene,'*')
        h1 = float(row[col])
        h2 = float(row[col+1])
        if h1 > 0 || h2 > 0
            h = (h1 + h2)/(float(h1>0) + float(h2>0))
            writedlm(f,[gene h],',')
        end
    end
    nothing
end

function make_datafiles(infolder,outfolder,label)
    histograms = readdir(infolder)
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    for file in histograms
        newfile = replace(file,label => "")
        cp(joinpath(infolder,file),joinpath(outfolder,newfile))
    end
end

function make_dataframe(folder::String,winners::String = "")
    rs = Array{Any,2}(undef,0,8)
    files =readdir(folder)
    n = 0
    for file in files
        if occursin("rate",file)
            t = parse(Float64,split(split(file,"T")[2],"_")[1])
            r,head = readdlm(joinpath(folder,file),',',header=true)
            r = [vcat(r[:,[1,2,3,4,5,10]], r[:,[1,6,7,8,9,10]])  [zeros(size(r)[1]); ones(size(r)[1])]  t*ones(2*size(r)[1])/60.]
            rs = vcat(rs,r)
            n += 1
        end
    end
    if winners != ""
        w = get_winners(winners,2*n)
        return DataFrame(Gene = rs[:,1],on = float.(rs[:,2]),off=float.(rs[:,3]),eject=float.(rs[:,4]),decay=float.(rs[:,5]),yield=float.(rs[:,6]),cond = Int.(rs[:,7]),time = float.(rs[:,8]),winner = w)
    else
        return DataFrame(Gene = rs[:,1],on = float.(rs[:,2]),off=float.(rs[:,3]),eject=float.(rs[:,4]),decay=float.(rs[:,5]),yield=float.(rs[:,6]),cond = Int.(rs[:,7]),time = float.(rs[:,8]))
    end
end

function get_winners(winners::String,n::Int)
    m,h = readdlm(winners,',',header=true)
    winner = repeat(m[:,end],n)
end

max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
geps = 1e-5

"""
Function to define inf. velocity of a horseshoe vortex

* `a`: first kin (relative to collocation point: `colpt .- kin`)
* `b`: second kin (relative to collocation point: `colpt .- kin`)
* `x`: wind direction (unit vector for adimensionalized results, equal to vectorial
freestream velocity if dimensional results are desired)

* return: the influence velocity over the collocation point
"""
function hshoe_vinf(a::AbstractVector, b::AbstractVector, x::AbstractVector)
    na=norm(a)
    nb=norm(b)
    ab=dot(a, b)

    d1=(na*nb+ab)
    if d1>-geps && d1<=0.0
        d1=-geps
    elseif d1<geps && d1>0.0
        d1=geps
    end

    d2=(na-dot(a, x))
    if d2>-geps && d2<=0.0
        d2=-geps
    elseif d2<geps && d2>0.0
        d2=geps
    end

    d3=(nb-dot(b, x))
    if d3>-geps && d3<=0.0
        d3=-geps
    elseif d3<geps && d3>0.0
        d3=geps
    end

    vbound = (cross(a, b)/d1)*(1.0/na+1.0/nb)
    vinf = cross(a, x)/(na*d2).-cross(b, x)/(nb*d3)

    vinf./=(4*pi)
    vbound./=(4*pi)

    return vinf, vbound
end

"""
Function to define inf. velocity of a finite vortex line segment

* `a`: first kin (relative to collocation point: `colpt .- kin`)
* `b`: second kin (relative to collocation point: `colpt .- kin`)

* return: the influence velocity over the collocation point
"""
function vline_segment_vinf(a::AbstractVector, b::AbstractVector)
    na=norm(a)
    nb=norm(b)
    ab=dot(a, b)

    d1=(na*nb+ab)
    if d1>-geps && d1<=0.0
        d1=-geps
    elseif d1<geps && d1>0.0
        d1=geps
    end

    vinf=(cross(a, b)/d1)*(1.0/na+1.0/nb)

    vinf/=(4*pi)

    return vinf
end

"""
Build linear aerodynamic problem for an aircraft
"""
function linear_problem(acft, kins, cols, norms, us, xhat, PGβ)

    b = [
        begin
            nu = norm(u)

            (
                - nu * strip.afl.CL0 - u ⋅ n * strip.afl.CLα
            )
        end * PGβ for (i, (u, n, strip)) in enumerate(zip(us, norms, acft.wing_strips))
    ]

    Fg = eltype(b)

    As = [
        Matrix{Fg}(undef, length(b), length(b)) for i = 1:3
    ]

    A = [
        begin
            vinfl, vbound = hshoe_vinf(cols[i] .- kins[j][1], cols[i] .- kins[j][2], xhat)

            for (id, A) in enumerate(As)
                A[i, j] = vinfl[id] + vbound[id]
            end

            C = acft.wing_strips[i].afl.CLα * PGβ * (
                (vinfl .+ vbound) ⋅
                norms[i]
            )

            if i == j
                C -= 2.0 / acft.wing_strips[i].c
            end

            C
        end for i = 1:length(acft.wing_strips), j = 1:length(acft.wing_strips)
    ]

    A, b, As

end

"""
Add downwash to velocities at collocation points
"""
function add_downwash!(us, Γs, As)

    dVs = [
        A * Γs for A in As
    ]

    for (id, dV) in enumerate(dVs)
        for (i, u) in enumerate(us)
            u[id] += dV[i]
        end
    end

end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
@doc raw"""
θ-pinch equilibrium in (x,y,z) coordinates with covariant components of
the vector potential given by
```math
A (x,y,z) = \frac{B_0}{2} \big( - y , \, x , \, 0 \big)^T
```
resulting in the magnetic field with covariant components
```math
B (x,y,z) = \big( 0 , \, 0 , \, B_0 \big)^T .
```

Parameters:
    B₀: B-field at magnetic axis
"""
struct ThetaPinch{T <: Number} <: CartesianEquilibrium
    name::String
    B₀::T

    function ThetaPinch{T}(B₀::T) where T <: Number
        new("ThetaPinchEquilibrium", B₀)
    end
end

ThetaPinch(B₀::T=1.0) where T <: Number = ThetaPinch{T}(B₀)

function Base.show(io::IO, equ::ThetaPinch)
    print(io, "θ-Pinch Equilibrium in (x,y,z) Coordinates with\n")
    print(io, "  B₀ = ", equ.B₀)
end

function R(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    r(x,equ)
end

function r(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    sqrt(r²(x,equ))
end

function r²(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    X(x,equ)^2 + Y(x,equ)^2
end

function θ(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    atan2(Y(x,equ), X(x,equ))
end

function ϕ(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    θ(x,equ)
end

function A₁(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    - equ.B₀ * Y(x,equ) / 2
end

function A₂(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    + equ.B₀ * X(x,equ) / 2
end

function A₃(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number}
    zero(eltype(x))
end

macro zpinch_equilibrium(R₀, B₀)
    generate_equilibrium_code(ThetaPinch(R₀, B₀); output=false)
end


text:
#   This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith
#   Copyright (c) 2021 Jingyue Wang
#
#   TODO: type-check. for example, if we define `b : matrix([1,2],[3,4])`, then
#   `cauchy_matrix([a,b,c])` gives very wierd result that can not be rendered
#   correctly as LaTeX on DomTerm.
export
    MExpr,
    @mx_str,
    @jx_str,
    parse,
    convert,
    get_mstr,
    clean_matrix_env,
    mexpr_to_tex,
    jexpr_to_mstr

import
    Base: parse, convert, *

using MathMLRepr

const infix_ops = [:+, :-, :*, :/, :^, ://, :.*, :∘]

isinfix(args) = args[1] in infix_ops && length(args) > 2

const jsym_to_m = Dict(
    :α           => "alpha",
    :β           => "beta",
    :δ           => "delta",
    :ϵ           => "epsilon",
    :ζ           => "zeta",
    :η           => "eta",
    :θ           => "theta",
    :ι           => "iota",
    :κ           => "kappa",
    :λ           => "lambda",
    :μ           => "mu",
    :ν           => "nu",
    :ξ           => "xi",
    :ρ           => "rho",
    :σ           => "sigma",
    :τ           => "tau",
    :υ           => "upsilon",
    :χ           => "chi",
    :ψ           => "psi",
    :ω           => "omega",
    :ℯ           => "%e",
    :pi          => "%pi",
    :π           => "%pi",
    :γ           => "%gamma",
    :eulergamma  => "%gamma",
    :golden      => "%phi",
    :φ           => "%phi",
    :im          => "%i",
    :Inf         => "inf"
)

jexpr_to_mstr(x::Symbol) = x in keys(jsym_to_m) ? jsym_to_m[x] : string(x)
jexpr_to_mstr(x::Rational) = string(x.num)*"/"*string(x.den)
jexpr_to_mstr(x::Union{AbstractFloat, Integer}) = string(x)
jexpr_to_mstr(x::String) = x
jexpr_to_mstr(x::LineNumberNode) = ""
jexpr_to_mstr(x::QuoteNode) = error("Can not process QuoteNode")
# return a string containing a valid Maxima statement
# We assume expr is NOT a block expression, so basically just one statement.
# The function can handle block expressions inside a function object.
function jexpr_to_mstr(expr::Expr; defeqn::Bool=false)
    o = ""
    if expr.head == :call
        operator = isinfix(expr.args) ? " $(expr.args[1]) " : ", "
        # operators that need special treatments to make a Maxima expression
        if operator == " // "; operator = "/"; end
        if operator == " ∘ "; operator = "."; end
        if operator == " .* "; operator = "*"; end
        if !isinfix(expr.args)
            o *= jexpr_to_mstr(expr.args[1])
        end
        o *= "("* (length(expr.args) > 1 ? jexpr_to_mstr(expr.args[2]) : "")
        for arg in expr.args[3:end] # if s > end args[3:end] is empty
            o *= operator*jexpr_to_mstr(arg)
        end
        o *= ")" #JW: we can replace `()` by `{}` for a component behind `^` in TeX string
    elseif expr.head == :(=)
        #TODO: save variable names and their type info
        if defeqn === true
            o *= jexpr_to_mstr(expr.args[1])*" = "*jexpr_to_mstr(expr.args[2])
        elseif expr.args[1] isa Symbol
            o *= jexpr_to_mstr(expr.args[1])*": "*jexpr_to_mstr(expr.args[2])
        else # should be a :call that implies a functin definition, e.g., f(x)=x^2
            fun_body = join(jexpr_to_mstrv(expr.args[2]), ",")
            o *= jexpr_to_mstr(expr.args[1])*":= "*fun_body
        end
    elseif expr.head == :block
        # we deal a limited case of :block, namely only one expression inside
        length(expr.args) > 2 && error("The RHS of `=` is too compiated to process")
        # args[1] is assumed to be a LineNumberNode
        !(expr.args[1] isa LineNumberNode) && error("Unknown expression in RHS")
        o *= jexpr_to_mstr(expr.args[2])
    elseif expr.head == :function
        o *= jexpr_to_mstr(expr.args[1])*" := block([]"
        # args[2] of a function always returns a Julia block expression, even for an
        # empty function (`function f(x) end`), it returns an empty block
        args = jexpr_to_mstrv(expr.args[2])
        for a in args; o *= ", "*a; end
        o *= ")"
    elseif expr.head == :row || expr.head == :vect
        o = "["*jexpr_to_mstr(expr.args[1])
        for a in expr.args[2:end]
            o *= ", "*jexpr_to_mstr(a)
        end
        o *= "]"
    elseif expr.head == :vcat
        o = "matrix("
        for a in expr.args[1:end-1]
            o *= jexpr_to_mstr(a)*", "
        end
        o *= jexpr_to_mstr(expr.args[end])*")"
    elseif expr.head == Symbol("'")
        o = "transpose("*jexpr_to_mstr(expr.args[1])*")"
    else
        # can not deal with QuoteNode either
        error("Nested :$(expr.head) block structure not supported by Maxima.jl")
    end

    return o
end

# for a single expression, return an array that only contains one `mstring`
# for a block expression, return an array that contains all expressions in the block
function jexpr_to_mstrv(expr::Expr)
    if expr.head in (:block, :toplevel)
        strvec = String[]
        for a in expr.args
            push!(strvec, jexpr_to_mstr(a))
        end
        # a LineNumeberNode yields an empty string.
        return clean_empty_ele(strvec)
    else
        return String[jexpr_to_mstr(expr)]
    end
end

function jexpr_to_mstr(x::Vector{T}) where T
    o = "["*jexpr_to_mstr(x[1])
    for a in x[2:end]
        o *= ", "*jexpr_to_mstr(a)
    end
    o *= "]"
end

function jexpr_to_mstr(Mat::Matrix{T}) where T
    m, n = size(Mat)
    o = "matrix("
    # add the first row
    o *= jt2ms_addrow(Mat[1, :])
    # add each of the rest of rows
    for i = 2:m
        o *=  ", "*jt2ms_addrow(Mat[i, :])
    end
    o *= ")"
    return o
end

# add a Julia vector as a row to a Maxima table
function jt2ms_addrow(v::Vector{T}) where T
    o = "["*jexpr_to_mstr(v[1])
    for j = 2:length(v)
        o *= ", "*jexpr_to_mstr(v[j])
    end
    o *= "]"
    return o
end

"""

A Maxima expression

## Summary:

type MExpr <: Any

## Fields:

strvec :: Vector{String}
"""
struct MExpr
    strvec::Vector{String}
end
MExpr(x::Rational) = MExpr("$(x.num)"*"/"*"$(x.den)")
MExpr(x::Number) = MExpr(string(x))
MExpr(x::String) = MExpr(push!(String[], x))
MExpr(x::Symbol) = MExpr(string(x))
"""
    MExpr(expr::Expr)

Convert a Julia expression `expr` to an equivalent Maxima expression

## Examples
```julia
julia> MExpr(:(sin(x*im) + cos(y*φ)))

                           cos(%phi y) + %i sinh(x)

```
"""
MExpr(expr::Expr) = MExpr(jexpr_to_mstrv(expr))

macro mx_str(x) MExpr(x) end
macro jx_str(x) MExpr(Meta.parse(x)) end

*(x::MExpr,  y::MExpr)  = MExpr(append!(String[], x.strvec, y.strvec))
*(x::MExpr,  y::String) = MExpr(append!(String[], x.strvec, [y]))
*(x::String, y::MExpr)  = MExpr(append!(String[x], y.strvec))

get_mstr(m::MExpr, i) = m.strvec[i]

clean_empty_ele(m::Vector{String}) = filter(!isempty, m)
allspace(x::AbstractString) = all(isspace, x)
clean_space_ele(m::Vector{String}) = filter(!allspace, m)
# for each mstring in a MExpr, if it is a command sequence, split it into
# a vector of single-commands. Maxima commands are separated by `$` or `;` in
# a sequence.
function normalize(m::MExpr)
    r = String[]
    for s in m.strvec
        p = split(s, ['$', ';'])
        # NOTE: `split()` returns `p`, a Vector{SubString{String}}, NOT a
        #Vector{String}. `append!()` makes `r` a Vector{String}
        append!(r, p)
    end
    return MExpr(r |> clean_empty_ele |> clean_space_ele)
end

const m_to_jstr = Dict(
    "alpha"   =>  "α",
    "beta"    =>  "β",
    "delta"   =>  "δ",
    "epsilon" =>  "ϵ",
    "zeta"    =>  "ζ",
    "eta"     =>  "η",
    "theta"   =>  "θ",
    "iota"    =>  "ι",
    "kappa"   =>  "κ",
    "lambda"  =>  "λ",
    "mu"      =>  "μ",
    "nu"      =>  "ν",
    "xi"      =>  "ξ",
    "rho"     =>  "ρ",
    "sigma"   =>  "σ",
    "tau"     =>  "τ",
    "upsilon" =>  "υ",
    "chi"     =>  "χ",
    "psi"     =>  "ψ",
    "omega"   =>  "ω",
    "%e"      =>  "ℯ",
    "%pi"     =>  "π",
    "%i"      =>  "im",
    "%gamma"  =>  "eulergamma",
    "%phi"    =>  "φ",
    "inf"     =>  "Inf",
    "minf"    =>  "-Inf"
)
"""
    parse(mexpr::MExpr)

Convert a Maxima expression into an equivalent Julia expression

## Examples
```julia
julia> parse(mx\"sin(%i*x)\")
:(im * sinh(x))

```
"""
#FIXME: I should do real parsing instead of naive text substitutions.
#NOTE: `parse()` lambdify a symbolic function, convert a symbolic number into
#      a concrete type (Int64, Float32,...).
function parse(m::MExpr)
    exprs = Any[]
    for ex in normalize(m).strvec
        for key in keys(m_to_jstr)
            #FIXME: we should parse the MExpr, not simply doing a text substitution.
            ex = replace(ex, key => m_to_jstr[key])
        end
        if occursin( ":=", ex)
            args = split(ex, ":=")
            push!(
                exprs,
                Expr(
                    :function,
                    parse(args[1]),  # e.g. parse("f(x)")
                    args[2] |> String |> MExpr |> parse))
        elseif occursin( ":", ex)
            args = split(ex, ":")
            push!(
                exprs,
                Expr(
                    :(=),
                    parse(args[1]),
                    args[2] |> String |> MExpr |> parse))
        elseif occursin( "block([], ", ex)
            x = replace(ex, "block([],", "") |> strip |> chop  # chop the last `)`
            xs = split(x, ",")
            l = Vector{Any}(length(x))
            for s in xs
                !isempty(s) && push!(l, s |> String |> MExpr |> parse)
            end
            push!(exprs, Expr(:block, l...))
        else
            push!(exprs, Meta.parse(ex))
        end
    end
    return length(exprs) == 1 ? exprs[1] : Expr(:block, exprs...)
end

convert(::Type{MExpr}, m::MExpr) = m
convert(::Type{Array{String,1}}, m::MExpr) = m.strvec
convert(::Type{String}, m::MExpr) = join(m.strvec, "; ")
convert(::Type{T}, m::MExpr) where T<:Number = eval(parse(m))
#IMPORTANT! we strip the leading and traling `$`s generated by Maxima
function convert(::Type{LATEX}, m::MExpr)
    check = "'(" * replace(convert(String, m), r";" => ")\$\n'(") * ")"
    write(ms.input, "$check\$\n print(ascii(4))\$")
    out = (readuntil(ms.output, EOT) |> String
                                     |> str -> rstrip(str, EOT))

    s = ""
    if occursin(synerr, out) || occursin(maxerr, out)
        @warn "Invalid Maxima expression"
        s = "Invalid Maxima Expression: "*out
    else
        write(ms, "tex1('("*replace(convert(String, m), r";" => "))\$\ntex1('(") * "))")
        texstr = read(ms)
        s = replace(texstr, r"\nfalse\n" => "")
        # Warning! quick hacks inside this function
        s = clean_matrix_env(s)
        s = String(strip(s, '$'))
    end
    return LATEX(s)
end
#JW: In stdlib 1.6.1 `display(x::MExpr) dereferences x into Any
convert(::Type{Any}, m::MExpr) = m
mexpr_to_tex(m::MExpr) = convert(LATEX, m)

# All tricks to post-process the TeX code generated by maxima are wrapped
# in this function
function clean_matrix_env(x::String)
    x = strip(x)
    x = chop(x, head=1, tail=1) # x retuend by `tex1()` is encluded by double-quote `"`
    if occursin("pmatrix", x)
        x = replace(x, r"\\\\ifx.*?\\\\fi"=>"") # del all the \ifx
        x = replace(x, r"\\\\cr +$"=>"") # del the last empty \cr
        x = replace(x, r"\\\\cr"=>"\\\\") # \cr -> \\
        x = replace(x, r"\\\\([a-z])"=>s"\\\1") # any double-slash command \\cmd => \cmd
        x = replace(x, r"\\\\([,!])"=>s"\\\1") # any double-slash command \\opr => \opr
        x = "\\begin{pmatrix}"*x*"\\end{pmatrix}"
    else
        x = replace(x, raw"\\\\"=>"\\")
    end
    return x
end