text:
using .JSort
using PyCall
import PyPlot; const plt = PyPlot;

function localalign(detector::ChannelDetector,
                    region; nbins=5000, emin=-150, emax=15000,
                    kwargs...)
    full, edges = histogram(detector, nbins=nbins)
    chns = goodchannels(detector)

    signal = pyimport("scipy.signal")
    smoother(x) = signal.savgol_filter(x, 51, 4)

    eclip(x) = clip(x, emin, emax)

    kw = get(kwargs, :a1, Dict())
    if get(kw, :plot, false)
        plothist(full, edges)
    end
    delta = get(kw, :delta, [3500, 3500])
    searchregion = eclip(region[1]-delta[1]):eclip(region[2]+delta[2])

    println("Aligning")

    coeff, _ = alignspectra(full, edges, lowregion=searchregion, lowsmoother=smoother,
                            lownumregions=get(kw, :nregions, 2),
                            lowregionwidth=get(kw, :swidth, 200),
                            lowsearchwidth=get(kw, :lwidth, 60),
                            plot=get(kw, :plot, false), order=1)
    coefficients = repeat([0.0 1.0], outer=(33, 1))
    coefficients[chns, :] .= hcat(coeff...)'
    calib = GammaCalibrator(coefficients)
    aligned = calib(detector, inplace=false)

    if get(kw, :plot, false)
        plot(aligned)
    end

    alignedhist = histogram(aligned, edges)

    kw = get(kwargs, :a2, Dict())
    delta = get(kw, :delta, [3500, 3500])
    searchregion = eclip(region[1]-delta[1]):eclip(region[2]+delta[2])

    coeff, _ = alignspectra(alignedhist, edges, lowregion=searchregion, lowsmoother=smoother,
                            lownumregions=get(kw, :nregions, 3),
                            lowregionwidth=get(kw, :swidth, 200),
                            lowsearchwidth=get(kw, :lwidth, -10),
                            plot=get(kw, :plot, false), order=1)

    coefficients = repeat([0.0 1.0], outer=(33, 1))
    coefficients[chns, :] .= hcat(coeff...)'
    combine!(calib, GammaCalibrator(coefficients))
    aligned = calib(detector, inplace=false)

    if get(kw, :plot, false)
      plot(aligned)
    end

    #S = sum(aligned)
    #hist = histogram(S, edges)
    #plot(S)
    # Energy - index is wrong
    #peak = findpeaks(hist[region], numpeaks=1)[1]
    #peak += first(region)
    #@show peak
    #@show peaks = map(x -> x(peak), invcalib.poly)
    #ax.scatter(peaks[chns], 1.5:length(chns)+0.5, marker="x", color="r", s=2)

    return calib
end

function regionalign(d::ChannelDetector, regions; nbins=1000,
                     emin=100, emax=2500, kwargs...)
    strips, edges = JSort.histogram(d, nbins=nbins)
    shift, gain = first(edges), step(edges)
    indextoe(x) = shift + gain*x
    etoindex(x) = round(Int, -shift/gain + 1/gain*x)
    #regions = [map(etoindex, region) for region in regions]
    eclip(x) = clip(x, emin, emax)
    delta = 400

    peaks = []
    for region in regions
        low, high = region[1], region[2]
        # Align all strips to the reference strip
        lowregion = eclip(low-delta):eclip(high+delta)
        @show lowregion
        println("Aligning")
        al, calib, _ = align(d, lowregion=lowregion, lownumregions=2,
                             lowregionwidth=200,
                             lowsearchwidth=80, order=1, kwargs...);
        println("Done aliging")
        println("Region peaking")
        # Find the peak(s) in the aligned region
        peak = regionpeaks(al, [region])
        println("Done region peaking")
        # Backshift the peaks to the original coordinate system
        peak = backshift(peak, goodchannels(d), invert(calib))
        push!(peaks, peak)
    end
    peaks
end

function regionpeaks(d::ChannelDetector, regions; nbins=5000)
    strips, edges = JSort.histogram(d, nbins=nbins)
    shift, gain = first(edges), step(edges)
    indextoe(x) = shift + gain*x
    etoindex(x) = round(Int, -shift/gain + 1/gain*x)
    regions = [map(etoindex, region) for region in regions]
    println("Regioning for $regions")
    peaks = regionpeaks(strips, regions)
    [map(indextoe, peak) for peak in peaks]
end

function regionpeaks(strips, regions)
    peaks = Vector{Int}[]
    for strip in strips
        peaks_ = regionpeaks(strip, regions)
        push!(peaks, peaks_)
    end
    peaks
end

function regionpeaks(strip::AbstractVector{<:Number}, regions)::Vector{Int}
    peaks = Int[]
    for region in regions
        peak = regionpeak(strip, region[1], region[2])
        push!(peaks, peak)
    end
    return peaks
end

function regionpeak(strip::AbstractVector{<:Number}, lower::Int, upper::Int)::Int
    region = strip[lower:upper]
    peak = findpeaks(region, numpeaks=1)[1]
    return peak + lower - 2 #BUG -2 is indicative of a offset bug
end

function backshift(points, channels, calib)
    shiftedpoints = similar(points)
    for i in 1:length(points)
        shiftedpoints[i] = calib[channels[i]](points[i])
    end
    shiftedpoints
end

function backshift(points::AbstractVector, coefficients::AbstractVector)
    shiftedpoints = similar(points)
    for i in 1:length(points)
        shift, gain = coefficients[i]
        poly(x) = x/gain - shift/gain
        shiftedpoints[i] = poly(points[i])
    end
    shiftedpoints
end

function backshift(point::Number, coefficients::AbstractVector)
    shiftedpoints = []
    for i in 1:length(points)
        shift, gain = coefficients[i]
        poly(x) = x/gain - shift/gain
        push!(shiftedpoints, poly(point))
    end
    shiftedpoints
end

function clip(x, lower, upper)
    if x < lower
        lower
    elseif x > upper
        upper
    else
        x
    end
end

function plothist(strips, edges; points=nothing)
    mat = cat(strips'..., dims=1)
    fig, ax = plt.subplots()
    ax.pcolormesh(edges, 1:length(strips)+1, mat .|> log10)
    if !isnothing(points)
        for points_ in points
            for (i, point) in enumerate(points_)
                ax.scatter(point, i+0.5, marker="x", color="r", s=3)
            end
        end
    end
    ax
end

function plothist(det::ChannelDetector; nbins=1000, points=nothing)
    det_strips, det_edges = JSort.histogram(det, nbins=nbins)
    plothist(det_strips, det_edges, points=points)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
# Pkg.add.(["TensorFlow", "Distributions", "ProgressMeter", "MLLabelUtils", "MLDataUtils"])
# # Pkg.add("MLDatasets")  # not registered yet
# Pkg.clone("https://github.com/JuliaML/MLDatasets.jl.git")
# Pkg.checkout("MLDataUtils", "dev")
# using Conda
# # Conda.add("tensorflow")
# # using PyCall
# # @pyimport Conda
# # @show Conda.PYTHONDIR

using Base.Test
using Distributions
using MLLabelUtils
using MLDataUtils
using MLDatasets
using ProgressMeter
using TensorFlow

#Training Hyper Parameter
const learning_rate = 0.001
const training_iters = 2 #Just two, becuase I don't have anything to stop overfitting and I don't got all day
const batch_size = 256
const display_step = 100 #How often to display the

# Network Parameters
const n_input = 28 # MNIST data input (img shape: 28*28)
const n_steps = 28 # timesteps
const n_hidden = 128 # hidden layer num of features
const n_classes = 10; # MNIST total classes (0-9 digits)

const traindata_raw, trainlabels_raw = MNIST.traindata();
@show size(traindata_raw)
@show size(trainlabels_raw)

imshow(x) = join(mapslices(join, (x->x ? 'X': ' ').(x'.>0), 2), "\n") |> print
@show trainlabels_raw[8]
imshow(traindata_raw[:,:,8])

"""Makes 1 hot, row encoded labels."""
encode(labels_raw) = convertlabel(LabelEnc.OneOfK, labels_raw, LabelEnc.NativeLabels(collect(0:9)),  LearnBase.ObsDim.First())

"""Prepares the data by encoding the labels and batching"""
prepared_batchs(data_raw, labels_raw) = eachbatch((data_raw, encode(labels_raw)), #Will zip these
                               batch_size,
                               (MLDataUtils.ObsDim.Last(), MLDataUtils.ObsDim.First())) #Slicing dimentions for each

@testset "data prep" begin

    @test encode([4,1,2,3,0]) == [0 0 0 0 1 0 0 0 0 0
                                  0 1 0 0 0 0 0 0 0 0
                                  0 0 1 0 0 0 0 0 0 0
                                  0 0 0 1 0 0 0 0 0 0
                                  1 0 0 0 0 0 0 0 0 0]

    data_b1, labels_b1 = first(prepared_batchs(traindata_raw, trainlabels_raw))
    @test size(data_b1) == (n_steps, n_input, batch_size)
    @test size(labels_b1) == (batch_size, n_classes)
end;

sess = Session(Graph())
X = placeholder(Float32, shape=[n_steps, n_input, batch_size])
Y_obs = placeholder(Float32, shape=[batch_size, n_classes])

variable_scope("model", initializer=Normal(0, 0.5)) do
    global W = get_variable("weights", [n_hidden, n_classes], Float32)
    global B = get_variable("bias", [n_classes], Float32)
end;

# Prepare data shape to match `rnn` function requirements
# Current data input shape: (n_steps, n_input, batch_size) from the way we declared X (and the way the data actually comes)
# Required shape: 'n_steps' tensors list of shape (batch_size, n_input)

x = transpose(X, Int32.([1, 3, 2].-1)) # Permuting batch_size and n_steps. (the -1 is to use 0 based indexing)
x = reshape(x, [n_steps*batch_size, n_input]) # Reshaping to (n_steps*batch_size, n_input)
x = split(1, n_steps, x) # Split to get a list of 'n_steps' tensors of shape (batch_size, n_input)
@show get_shape.(x);

Hs, states = nn.rnn(nn.rnn_cell.LSTMCell(n_hidden), x; dtype=Float32);
Y_pred = nn.softmax(Hs[end]*W + B)

@show get_shape(Y_obs)
@show get_shape(Y_pred);

cost = reduce_mean(-reduce_sum(Y_obs.*log(Y_pred), reduction_indices=[1])) #cross entropy
@show get_shape(Y_obs.*log(Y_pred))
@show get_shape(cost) #Should be [] as it should be a scalar

optimizer = train.minimize(train.AdamOptimizer(learning_rate), cost)

correct_prediction = indmax(Y_obs, 2) .== indmax(Y_pred, 2)
@show get_shape(correct_prediction)
accuracy = reduce_mean(cast(correct_prediction, Float32));

run(sess, initialize_all_variables())

kk=0
for jj in 1:training_iters
    for (xs, ys) in prepared_batchs(traindata_raw, trainlabels_raw)
        run(sess, optimizer,  Dict(X=>xs, Y_obs=>ys))
        kk+=1
        if kk % display_step == 1
            train_accuracy, train_cost = run(sess, [accuracy, cost], Dict(X=>xs, Y_obs=>ys))
            info("step $(kk*batch_size), loss = $(train_cost),  accuracy $(train_accuracy)")
        end
    end
end

ext: jl


text:
"""
    PyIO(o; own=false, text=missing, buflen=4096)

Wrap the Python byte-based IO stream `o` as a Julia IO stream.

When this goes out of scope and is finalized, it is automatically flushed. If `own=true` then it is also closed.

If `text=false` then `o` must be a binary stream and arbitrary binary I/O is possible. If `text=true` then `o` must be a text stream and only UTF-8 must be written (i.e. use `print` not `write`). If `text` is not specified then it is chosen automatically. If `o` is a text stream and you really need a binary stream, then often `PyIO(o.buffer)` will work.

For efficiency, reads and writes are buffered before being sent to `o`. The size of the buffer is `buflen`.
"""
mutable struct PyIO <: IO
    ref::PyRef
    # true to close the file automatically
    own::Bool
    # true if `o` is text, false if binary
    text::Bool
    # true if we are definitely at the end of the file; false if we are not or don't know
    eof::Bool
    # input buffer
    ibuflen::Int
    ibuf::Vector{UInt8}
    # output buffer
    obuflen::Int
    obuf::Vector{UInt8}

    function PyIO(
        o;
        own::Bool = false,
        text::Union{Missing,Bool} = missing,
        buflen::Integer = 4096,
        ibuflen = buflen,
        obuflen = buflen,
    )
        io = new(
            PyRef(o),
            own,
            text === missing ? pyisinstance(o, pyiomodule().TextIOBase) : text,
            false,
            ibuflen,
            UInt8[],
            obuflen,
            UInt8[],
        )
        finalizer(io) do io
            if CONFIG.isinitialized
                io.own ? close(io) : flush(io)
            end
        end
        io
    end
end
export PyIO

ispyreftype(::Type{PyIO}) = true
pyptr(io::PyIO) = pyptr(io.ref)
Base.unsafe_convert(::Type{CPyPtr}, io::PyIO) = checknull(pyptr(io))
C.PyObject_TryConvert__initial(o, ::Type{PyIO}) = C.putresult(PyIO(pyborrowedref(o)))

"""
    PyIO(f, o; ...)

Equivalent to `f(PyIO(o; ...))` except the stream is automatically flushed or closed according to `own`.

For use in a `do` block, like `PyIO(o, ...) do io; ...; end`.
"""
function PyIO(f::Function, o; opts...)
    io = PyIO(o; opts...)
    try
        return f(io)
    finally
        io.own ? close(io) : flush(io)
    end
end

# If obuf is non-empty, write it to the underlying stream.
function putobuf(io::PyIO)
    if !isempty(io.obuf)
        @py `$io.write($(io.text ? pystr(io.obuf) : pybytes(io.obuf)))`
        empty!(io.obuf)
    end
    nothing
end

# If ibuf is empty, read some more from the underlying stream.
# After this call, if ibuf is empty then we are at EOF.
# TODO: in binary mode, `io.readinto()` to avoid copying data
function getibuf(io::PyIO)
    if isempty(io.ibuf)
        append!(io.ibuf, @pyv `$io.read($(io.ibuflen))`::Vector{UInt8})
    end
    nothing
end

function Base.flush(io::PyIO)
    putobuf(io)
    @py `$io.flush()`
    nothing
end

Base.close(io::PyIO) = (flush(io); @py `$io.close()`; nothing)

function Base.eof(io::PyIO)
    if io.eof
        true
    elseif !isempty(io.ibuf)
        false
    else
        getibuf(io)
        io.eof = isempty(io.ibuf)
    end
end

Base.fd(io::PyIO) = @pyv `$io.fileno()`::Int

Base.isreadable(io::PyIO) = @pyv `$io.readable()`::Bool

Base.iswritable(io::PyIO) = @pyv `$io.writable()`::Bool

Base.isopen(io::PyIO) = @pyv `not $io.closed`::Bool

function Base.unsafe_write(io::PyIO, ptr::Ptr{UInt8}, n::UInt)
    ntodo = n
    while true
        nroom = max(0, io.obuflen - length(io.obuf))
        if ntodo < nroom
            append!(io.obuf, unsafe_wrap(Array, ptr, ntodo))
            return n
        else
            append!(io.obuf, unsafe_wrap(Array, ptr, nroom))
            putobuf(io)
            ptr += nroom
            ntodo -= nroom
        end
    end
end

function Base.write(io::PyIO, c::UInt8)
    push!(io.obuf, c)
    if length(io.obuf) ≥ io.obuflen
        putobuf(io)
    end
end

function Base.unsafe_read(io::PyIO, ptr::Ptr{UInt8}, n::UInt)
    ntodo = n
    while true
        navail = length(io.ibuf)
        if ntodo ≤ navail
            unsafe_wrap(Array, ptr, ntodo) .= splice!(io.ibuf, 1:ntodo)
            return
        else
            unsafe_wrap(Array, ptr, navail) .= io.ibuf
            ptr += navail
            ntodo -= navail
            empty!(io.ibuf)
            eof(io) && throw(EOFError())
        end
    end
end

function Base.read(io::PyIO, ::Type{UInt8})
    eof(io) && throw(EOFError())
    popfirst!(io.ibuf)
end

function Base.seek(io::PyIO, pos::Integer)
    putobuf(io)
    empty!(io.ibuf)
    io.eof = false
    @py `$io.seek($pos, 0)`
    io
end

function Base.truncate(io::PyIO, pos::Integer)
    seek(io, position(io))
    @py `$io.truncate($pos)`
    io
end

function Base.seekstart(io::PyIO)
    putobuf(io)
    empty!(io.ibuf)
    io.eof = false
    @py `$io.seek(0, 0)`
    io
end

function Base.seekend(io::PyIO)
    putobuf(io)
    empty!(io.ibuf)
    io.eof = false
    @py `$io.seek(0, 2)`
    io
end

function Base.skip(io::PyIO, n::Integer)
    putobuf(io)
    if io.text
        if n ≥ 0
            read(io, n)
        else
            error("`skip(io, n)` for text PyIO streams only implemented for positive `n`")
        end
    else
        if 0 ≤ n ≤ io.ibuflen
            read(io, n)
        else
            @py `$io.seek($(n - length(io.ibuf)), 1)`
            empty!(io.ibuf)
            io.eof = false
        end
    end
end

function Base.position(io::PyIO)
    putobuf(io)
    if io.text
        if isempty(io.ibuf)
            @pyv `$io.tell()`::Int
        else
            error(
                "`position(io)` text PyIO streams only implemented for empty input buffer (e.g. do `read(io, length(io.ibuf))` first)",
            )
        end
    else
        (@pyv `$io.tell()`::Int) - length(io.ibuf)
    end
end


text:
#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2021-Jul-08: add leaf level SIF simulation
#     2021-Jul-08: use mat_b and mat_f for SIF at backward and forward directions
#     2021-Aug-05: add option to sumulate SIF in photon to photon mode
#     2021-Oct-22: refactor the function to leaf_SIF to return the SIFs directly
#     2022-Jan-13: use LeafBiophysics directly in the function rather than Leaf
#     2022-Feb-02: fix documentation
#     2022-Feb-02: unpack CONSTANTS only
#
#######################################################################################################################################################################################################
"""

    leaf_SIF(bio::LeafBiophysics{FT}, wls::WaveLengthSet{FT}, rad::HyperspectralRadiation{FT}, ϕ::FT = FT(0.01); ϕ_photon::Bool = true) where {FT<:AbstractFloat}

Return the leaf level SIF at backward and forward directions, given
- `bio` `ClimaCache.LeafBiophysics` type struct that contains leaf biophysical parameters
- `wls` `ClimaCache.WaveLengthSet` type struct that contains wave length bins
- `rad` `ClimaCache.HyperspectralRadiation` type struct that contains incoming radiation information
- `ϕ` Fluorescence quantum yield
- `ϕ_photon` If true (default), convert photon to photon when computing SIF; otherwise, convert energy to energy

---
# Examples
```julia
wls = WaveLengthSet{Float64}();
bio = LeafBiophysics{Float64}(wls);
rad = HyperspectralRadiation{Float64}();
sif_b,sif_f = leaf_SIF(bio, wls, rad, 0.01);
sif_b,sif_f = leaf_SIF(bio, wls, rad, 0.01; ϕ_photon=false);
```
"""
function leaf_SIF(bio::LeafBiophysics{FT}, wls::WaveLengthSet{FT}, rad::HyperspectralRadiation{FT}, ϕ::FT = FT(0.01); ϕ_photon::Bool = true) where {FT<:AbstractFloat}
    @unpack IΛ_SIFE, ΔΛ_SIFE, Λ_SIF, Λ_SIFE = wls;

    # calculate the excitation energy and photons
    _e_excitation = (view(rad.e_direct , IΛ_SIFE) .+ view(rad.e_diffuse, IΛ_SIFE)) .* ΔΛ_SIFE;

    # convert energy to energy using the matrices
    if !ϕ_photon
        _sif_b = bio.mat_b * _e_excitation * ϕ / FT(pi);
        _sif_f = bio.mat_f * _e_excitation * ϕ / FT(pi);

        return _sif_b, _sif_f
    end;

    # convert energy to photon
    _phot_excitation = photon.(Λ_SIFE, _e_excitation);

    # convert photon to photon using the matrices
    _phot_b = bio.mat_b * _phot_excitation * ϕ / FT(pi);
    _phot_f = bio.mat_f * _phot_excitation * ϕ / FT(pi);

    # convert photon to back to energy
    _sif_b = energy.(Λ_SIF, _phot_b);
    _sif_f = energy.(Λ_SIF, _phot_f);

    return _sif_b, _sif_f
end


text:
# unifbase.jl
# common functions for unify on dvc or naive

#include("misc.jl")
#include("primitives.jl")
#include("types.jl")

## unification
"""
putasubst puts a pair of v,t to a substituion
"""
function putasubst(vars::Vlist, vt::Symbol, tm::Term, subst::Tlist)
 rs = map(x->if x==vt;tm else x end, subst)
end

function loopcheck(vt::Symbol, et::Number; istop=false)::Bool
  return false
end

function loopcheck(vt::Symbol, et::Symbol; istop=false)::Bool
  return false
end

function loopcheck(vt::Symbol, et::Expr; istop=true)::Bool
  for arg in et.args
   if typeof(arg) == Symbol
    (vt == arg) && throw(Loop(vt,et,:loopcheckse))
   else
    loopcheck(vt, arg, istop=false)&& return throw(Loop(vt,et,:loopcheckse))
   end
  end
  return false
end

function maketlist(vars::Vlist, t1::Symbol, t2::Symbol)::Tlist
  map(x->if x==t1;t2 else t1 end, vars)
end

"""
replace the ix'th element of tlist with t
without typechecking
"""
function putarray(tlist::Array, ix::Number, t::Any)
  map(i->if i==ix; t else tlist[i] end, 1:length(tlist))
end

function vindex(vars::Vlist, v::Symbol)
 ix = findfirst(x->x==v,vars)
 if ix == nothing; return 0 end
 return ix
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
#module Vector2d

importall Base
export Vector2D, convert, dot, getindex, norm, normsq, normalize

@doc """A simple type for a 2-component vector.
Previous tests seemed to show that this version was faster than
e.g. ImmutableArrays. But these performance tests should be redone.

This should presumably be a subtype of an abstract array type to avoid
so much rewriting of preexisting methods. """ ->
immutable Vector2D{T}
	x::T
	y::T
end

Vector2D{T}(v::Vector{T}) = Vector2D{T}(v[1], v[2])
#Vector2D{T}(v::Vector) = Vector2D{T}(v[1], v[2])
#Vector2D(v::Vector) = Vector2D(v[1], v[2])

+(v::Vector2D, w::Vector2D) = Vector2D(v.x+w.x, v.y+w.y)
+{T}(v::Vector2D{T}, w::Vector2D{T}) = Vector2D{T}(v.x+w.x, v.y+w.y)
#+{S,T}(v::Vector2D{S}, w::Vector2D{T}) = Vector2D{promote_rule(T,S)}(v.x+w.x, v.y+w.y)
-{T}(v::Vector2D{T}, w::Vector2D{T}) = Vector2D{T}(v.x-w.x, v.y-w.y)
*{T}(v::Vector2D{T}, lamb::Number) = Vector2D{T}(v.x*lamb, v.y*lamb)
*{T}(lamb::Number, v::Vector2D{T}) = v*lamb
/{T}(v::Vector2D{T}, lamb::Number) = Vector2D{T}(v.x/lamb, v.y/lamb)

=={T}(v::Vector2D{T}, w::Vector2D{T}) = v.x==w.x && v.y==w.y

size{T}(v::Vector2D{T}) = 2
dot{T}(v::Vector2D{T}, w::Vector2D{T}) = v.x*w.x + v.y*w.y
norm{T}(v::Vector2D{T}) = √(v ⋅ v)
normsq{T}(v::Vector2D{T}) = v ⋅ v
normalize{T}(v::Vector2D{T}) = v / norm(v)

getindex{T}(v::Vector2D{T}, i) = (i==1) ? v.x : v.y

atan2(v::Vector2D) = atan2(v[2], v[1])  # no need to remember the order in atan2

#There are problems with this `convert`
convert{T}(::Type{Vector2D{T}}, v::Vector{T}) = Vector2D{T}(v[1], v[2])
convert(::Type{Vector2D}, v::Vector) = Vector2D(v[1], v[2])

#end

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
using CSV, DataFrames, Dates

function load_board(filename)
    board = CSV.read(filename, DataFrame, header = 0) |> Array{Int8, 2}
    return board
end

function print_board(board::Array{Int8, 2})
    for i in 1:9
        if (i-1) % 3 == 0 && i != 0
            println("- - - - - - - - - - - -")
        end

        for j in 1:9
            if (j-1) % 3 == 0 && j != 1
                print(" | ")
            end
            if j == 9
                println(board[i, j])
            else
                print("$(board[i, j]) ")
            end
        end
    end
end

function solve(board::Array{Int8, 2})
    find = find_empty(board)
    if find == (10, 10)
        return true
    else
        row, col = find
        for i in 1:9
            if valid(board, i, (row, col))
                board[row, col] = i
                if solve(board)
                    return true
                else
                    board[row, col] = 0
                end
            end
        end
    end
    return false
end

function valid(board::Array{Int8, 2}, num::Int8, pos::Tuple{Int8, 1})
    # Check column
    num_in_col = findall( x -> x == num, board[:pos[2]] )
    if length(num_in_col) > 0
        return false
    end

    # Check row
    num_in_row = findall( x -> x == num, board[pos[1]:end] )
    if length(num_in_row) > 0
        return false
    end

    # Check box
    coordinates = Dict( 1 => 1, 2 => 1, 3 => 1,
                        4 => 4, 5 => 4, 6 => 4,
                        7 => 7, 8 => 7, 9 => 7)
    box = board[coordinates[pos[1]]:(coordinates[pos[1]]+2),
          coordinates[pos[2]]:(coordinates[pos[2]]+2)]
    num_in_box = findall( x -> x == num, box )

    if length(num_in_box) > 0
        return false
    end

    return true
end

function find_empty(board::Array{Int8, 2})
    for i in 1:9
        for j in 1:9
            if board[i, j] == 0
                return (i, j)  # row, col
            end
        end
    end
    return (10,10)
end

sudoku = load_board(raw"C:\Users\kaszt\github\Julia_Sudoku\9x9.csv")
println("\nOriginal:\n")
print_board(sudoku)
start_time = Time(Dates.now())
solve(sudoku)
end_time = Time(Dates.now())
println("\nSolved:\n")
print_board(sudoku)
time_taken = end_time - start_time
println(Dates.Millisecond(time_taken))


text:
export indep_poly

"""
`indep_poly(G)` returns the independence polynomial of the
`SimpleGraph` `G`.
"""
function indep_poly(G::SimpleGraph, cache_flag::Bool=true)
    if cache_flag && cache_check(G,:indep_poly)
        return cache_recall(G,:indep_poly)
    end

    if NV(G)==0
        p =  Poly([1])
        if cache_flag
            cache_save(G,:indep_poly,p)
        end
        return p
    end
    if NE(G)==0
        p = Poly([1,1])^NV(G)
        if cache_flag
            cache_save(G,:indep_poly,p)
        end
        return p
    end

    if is_connected(G)
        v = first(G.V)  # get any edge
        Nv = G[v]

        G1 = deepcopy(G)
        SimpleGraphs.delete!(G1,v)
        p1 = indep_poly(G1,false)
        for w in Nv
            SimpleGraphs.delete!(G1,w)
        end
        p2 = indep_poly(G1,false)

        p = p1 + Poly([0,1])*p2
        if cache_flag
            cache_save(G,:indep_poly,p)
        end
        return p
    end

    comps = parts(components(G))
    p = Poly([1])
    for S in comps
        H = induce(G,S)
        pH = indep_poly(H,false)
        p *= pH
    end
    if cache_flag
        cache_save(G,:indep_poly,p)
    end
    return p
end


text:
type Objective
    func
    args
    dim
    inherit
    constraint
    history

    function Objective(dim; func= Nullable(), args=Nullable(),
        inherit=obj_default_inherit, constraint=Nullable())
        return new(func, args, dim, obj_default_inherit, constraint, Float64[])
    end
end

function obj_construct_solution(objective::Objective, x; parent=Nullable())
    sol = Solution()
    sol.x = x
    sol.attach = objective.inherit(parent=parent)
    return sol
end

# evaluate the objective function of a solution
function obj_eval(objective, solution)
    solution.value = objective.func(solution, objective.args)
    push!(objective.history, solution.value)
    return solution.value
end

function obj_eval_constraint(objective, solution)
  #Todo
end

function obj_default_inherit(; parent=Nullable())
    return Nullable()
end

function get_history_bestsofar(objective)
    history_bestsofar = Float64[]
    bestsofar = Inf
    for i in 1:length(objective.history)
        if objective.history[i] < bestsofar
            bestsofar = objective.history[i]
        end
        push!(history_bestsofar, bestsofar)
    end
    return history_bestsofar
end

function obj_clean_history(objective::Objective)
    objective.history = []
end


text:
using SciMLTutorials
SciMLTutorials.weave_file(".","test.jmd")


text:
using NLopt
"""
`iscausebf(ω, c, e, iset; kwargs...)`

But-for causality.

Is `c` (provably) the actual cause of `e` in the causal world `ω`, allowing
interventions on variables defined in  `isset`
"""
function iscausebf(ω::Ω, c::RandVar, e::RandVar, iset;
                   sizes = [length(i(ω)) for i in iset],
                   proj = identity,
                   optargs...)
  @pre Bool(c(ω)) && Bool(e(ω)) "Both cause and effect must be true in ω"

  # Construct a loss function that is 1 when the constraint is satisfied
  loss = let rngs = splitvec(sizes), butfor = Omega.ciid(ω -> !c(ω) & !e(ω)), proj = proj
    function loss(vec, grad)
      # Since optimization requries vectors, we must do some book keeping
      replmap = Dict(iset[i] => proj(vec[rng]) for (i, rng) in enumerate(rngs))
      # In loss function we construct intervened mode where intervenable variables
      # in iset have been set to values under consideration
      butforint = replace(butfor, replmap)

      # Get a score of the degree tow hich the effect `e` stil holds (in log scale)
      Omega.logerr(butforint(ω))
    end
  end
  # optimize this loss function
  nlopt(loss, ω, iset, sizes; optargs...)
end

"Sequence of ranges of lengths in lengths"
function splitvec(lengths)
  i = 1
  ll = UnitRange{Int}[]
  for l in lengths
    b = i + l
    push!(ll, i:b - 1)
    i = b
  end
  ll
end

function nlopt(loss, ω, iset, sizes;
               n = prod(sizes),
               eq = isapprox,
               NlOptAlg = :LN_COBYLA,
               init = () -> rand(n))
  # n = prod(sizes)
  opt = NLopt.Opt(NlOptAlg, n)
  NLopt.max_objective!(opt, loss)
  NLopt.stopval!(opt, -0.0)
  NLopt.initial_step!(opt, 0.5)
  initreplaces = init()
  # loss(initreplaces, rand(3))
  (minf, ωvecoptim, ret) = NLopt.optimize(opt, initreplaces)
  if ret == :FORCED_STOP
    error("Forced Stop")
  end
  eq(minf, 0.0)
end

# Specify map from vector of values to intervention values
# Separate optimization criteria
# Return true false
# Have some tolerance in equality?
# Is the fact that the temperature is warm
# Implement Pearl things.  leep other things constant?
#


text:
module CellAnalyzerOld
using ImageView
using Images
using ColorTypes
using FixedPointNumbers
using HDF5

export myImgConvert, myDoG, myNormalize, myInvert, myShow, myNaiveMaxDetect, myMaxProject, myMeanImageFromFiles, mySegmentation, mySizeFilter, processFrame, extractMeanSignal, exportSignal, extractSignal, myStack2HDF5

function myImgConvert(img::Images.Image{ColorTypes.Gray{FixedPointNumbers.UfixedBase{UInt16,16}},3,Array{ColorTypes.Gray{FixedPointNumbers.UfixedBase{UInt16,16}},3}})

    nx, ny, nz = size(img)
    dat = data(img)
    mydat = zeros(Float64, (nx,ny,nz))
    for x in 1:nx
        for y in 1:ny
            for z in 1:nz
                mydat[x,y,z] = Float64(dat[x,y,z])
            end
        end
    end
    mydat
end

function myLabelConvert(img::Images.Image{ColorTypes.Gray{FixedPointNumbers.UfixedBase{UInt16,16}},3,Array{ColorTypes.Gray{FixedPointNumbers.UfixedBase{UInt16,16}},3}})

    nx, ny, nz = size(img)
    dat = raw(img)
    mydat = zeros(Float64, (nx,ny,nz))
    for x in 1:nx
        for y in 1:ny
            for z in 1:nz
                mydat[x,y,z] = Float64(dat[x,y,z])
            end
        end
    end
    mydat
end

function myDoG(img::Array{Float64,3};dx=0.46, dy=0.46, dz = 1., rmin = 5,rmax = 10)
    println("compute scale-normalized DoG")
    scaling = dx / dz
    large = imfilter_gaussian(img, [rmax, rmax, rmax * scaling])
    small = imfilter_gaussian(img, [rmin, rmin, rmin * scaling])

    (rmin/(rmax-rmin))*(large - small)
end

function myNormalize(img::Array{Float64,3})
    minval = minimum(img)
    maxval = maximum(img)
    valrange = maxval - minval
    map(x-> (x-minval)/valrange, img)
end

function myInvert(img::Array{Float64,3})
    map(x->1-x, img)
end

function myShow(img::Array{Float64,3})
    view(grayim(img))
end

function myNaiveMaxDetect(img::Array{Float64,3})
    nx,ny,nz = size(img)
    maximg = Array(Float64, size(img))

    for x in 2:nx-1
        for y in 2:ny-1
            for z in 2:nz-1
                if img[x,y,z] >= maximum(
                                        [img[x+1,y,z],img[x-1,y,z],
                                         img[x,y+1,z],img[x,y-1,z],
                                         img[x+1,y+1,z],img[x+1,y-1,z],img[x-1,y-1,z],img[x-1,y+1,z],
                                         img[x,y,z+1],img[x+1,y,z+1],img[x-1,y,z+1],
                                         img[x,y+1,z+1],img[x,y-1,z+1],
                                         img[x+1,y+1,z+1],img[x+1,y-1,z+1],img[x-1,y-1,z+1],img[x-1,y+1,z+1],
                                         img[x,y,z-1],img[x+1,y,z-1],img[x-1,y,z-1],
                                         img[x,y+1,z-1],img[x,y-1,z-1],
                                         img[x+1,y+1,z-1],img[x+1,y-1,z-1],img[x-1,y-1,z-1],img[x-1,y+1,z-1],
                                         ])
                    maximg[x,y,z] = 1.
                else
                    maximg[x,y,z] = 0.
                end
            end
        end
    end
    maximg
end

function myBinarize(img::Array{Float64,3}, t::Float64)
    map(x-> x < t ? 0. : 1., img)
end

function myBlobDetect(img::Array{Float64,3},mask::Array{Float64,3})

    nx, ny, nz = size(img)
    maximg = Array(Float64, size(img))

    pixlist = Array(Float64, (Int64(sum(mask[2:nx-1,2:ny-1,2:nz-1])),4))

    i = 1

    for x in 2:nx-1
        for y in 2:ny-1
            for z in 2:nz-1
                if mask[x,y,z] == 1.
                    pixlist[i,:] = [x,y,z,img[x,y,z]]
                    i = i + 1
                end

            end
        end
    end

    pixlist = sortrows(pixlist, by = x->(x[4]), rev = true)
    maxlabel = 0

    for k in 1:size(pixlist,1)
        x,y,z = Int64(pixlist[k,1]),Int64(pixlist[k,2]),Int64(pixlist[k,3])
        println((x,y,z))
        val = pixlist[k,4]
        localmaxlabel = -1

        for xx in -1:1
            for yy in -1:1
                for zz in -1:1
                    if xx == yy == zz == 0
                        continue
                    else
                        if img[x+xx,y+yy,z+zz] < val
                            continue
                        elseif img[x+xx,y+yy,z+zz] > val && localmaxlabel == -1
                            localmaxlabel = maximg[x+xx,y+yy,z+zz]
                        elseif img[x+xx, y+yy, z+zz] > val && maximg[x+xx, y+yy, z+zz] == 0
                            localmaxlabel = 0
                            break
                        elseif img[x+xx,y+yy,z+zz] > val && localmaxlabel >0 && localmaxlabel == maximg[x+xx,y+yy,z+zz]
                            continue
                        elseif img[x+xx,y+yy,z+zz] > val && localmaxlabel >0 && localmaxlabel != maximg[x+xx,y+yy,z+zz]
                            localmaxlabel = 0
                            break
                        elseif img[x+xx,y+yy,z+zz] == val
                            localmaxlabel = -2
                        end
                    end
                end
            end
        end

        if localmaxlabel == -1
            maximg[x,y,z] = maxlabel + 1
            maxlabel = maxlabel + 1
        elseif localmaxlabel == -2
            maximg[x,y,z] = 0
        else
            maximg[x,y,z] = localmaxlabel
        end
    end
    maximg
end

function myStackExport(img::Array{Float64,3},dir::AbstractString)

    nx, ny, nz = size(img)

    for z in 1:nz
        ImageView.imwrite(ImageView.grayim(convert(Array{UInt16,2},img[:,:,z])),string(dir,lpad(string(z),4,"0"),".tif"))
    end
end

function myStack2HDF5(img::Array{Float64,3},fname::AbstractString)
    fid = HDF5.h5open(fname, "w")
    write(fid, "img", convert(Array{UInt16,3},img))
    close(fid)
end

function myMaxProject(img::Array{Float64,3})
    nx,ny,nz = size(img)
    mmax = zeros(Float64,(nx,ny))
    for x in 1:nx
        for y in 1:ny
            mmax[x,y] = maximum(img[x,y,:])
        end
    end
    mmax
end

function myMeanImageFromFiles(fileNames::Array{ASCIIString,1})
    img0 = ImageView.imread(fileNames[1]);
    img0dat = myImgConvert(img0);
    nx,ny,nz = size(img0dat);
    meanimg = img0dat * 1/length(fileNames)

    for f in fileNames[2,:]
            img = ImageView.imread(f);
            imgdat = myImgConvert(img);
            meanimg = meanimg + imgdat * 1/length(fileNames)
    end
    meanimg
 end

function myImagePreprocess(dir)
    #fnames = map(x->string(dir, x), readdir(dir)[2:end])
    fnames = map(x->string(dir,x), filter(x->contains(x,".tif"), readdir(dir)));
    println(fnames)
    mimg = myMeanImageFromFiles(fnames);
    nimg = myNormalize(mimg);
    bg = Images.imfilter_gaussian(nimg, [50,50,20])
    fg = Images.imfilter_gaussian(nimg, [2, 2, 1])
    img = fg - bg;
    myNormalize(img)
end

function mySegmentation(dir::AbstractString; t = 0.2)
    img = myImagePreprocess(dir)
    dog = myNormalize(myInvert(myDoG(img)))
    bin = myBinarize(dog, t)
    res = myBlobDetect(dog, bin)
end

function mySizeFilter(limg::Array{Float64,3}; minSize = 100, maxSize = 50000)
    #labels = unqiue(limg)
    sizes  = Dict{Float64,Int64}()
    #nx, ny, nz = size(limg)
    for lab in limg
        sizes[lab] = get(sizes, lab, 0) + 1
    end
    sizes

    map(x -> sizes[x] > minSize && sizes[x] < maxSize ? x : 0., limg)
end

function extractMeanSignal(img::Array{Float64,3}, label::Array{Float64,3})

    signal = Dict{Float64,Float64}()
    nlab = Dict{Float64,Int64}()

    nx,ny,nz = size(img)

    for x in 1:nx
        for y in 1:ny
            for z in 1:nz
                if label[x,y,z] == 0.
                    continue
                else
                    signal[label[x,y,z]] = get(signal, label[x,y,z], 0) + img[x,y,z]
                    nlab[label[x,y,z]] = get(nlab, label[x,y,z],0) + 1
                end
            end
        end
    end

    for k in keys(signal)
        signal[k] = signal[k] / nlab[k]
    end

    signal
end

function processFrame(fname::AbstractString, label::Array{Float64,3})
    sig = Images.imread(fname);
    sig = myImgConvert(sig);
    res = extractMeanSignal(sig, label);
end

function extractSignals(dir::AbstractString, label::Array{Float64,3})

    fnames = filter(x->contains(x,".tif"), readdir(dir));

    signalDict = [id => zeros(length(fnames),1) for id in unique(label)]

    for (t, f) in enumerate(fnames)
        print(string(t), " ")
        signalsAtT = processFrame(string(dir,f), label)
        for lab in keys(signalsAtT)
            signalDict[lab][t] = signalsAtT[lab]
        end

    end

    signalDict
end

function exportSignal(fname::AbstractString, signalDict::Dict)
    fid = h5open(fname, "w")
    for id in keys(signalDict)
        if id == 0.0
            continue
        else
            write(fid, string("l_",Int64(id)), signalDict[id][:])
        end
    end

    close(fid)
end

end

max_stars_repo_name: nscherf/optical-cell-mapping

max_issues_repo_name: nscherf/optical-cell-mapping

max_forks_repo_name: nscherf/optical-cell-mapping


text:
export CR3BP, HandcodedCR3BP
export jacobi

#---------------------------#
# CR3BP EQUATIONS OF MOTION #
#---------------------------#
struct _CR3BP_ODEFunctions{S,F,F2} <: Abstract_ModelODEFunctions
    ode_system :: S
    ode_f      :: F
    ode_stm_f  :: F2
end

function ModelingToolkit.ODESystem(::Type{_CR3BP_ODEFunctions})
    # Build from the ER3BP equations (with eccentricity = 0 for circular)
    eqs_er3bp = ODESystem(_ER3BP_ODEFunctions)
    (μ, e) = parameters(eqs_er3bp)
    eqs = [eq.lhs ~ simplify(substitute(eq.rhs, e => 0)) for eq in equations(eqs_er3bp)]
    return ODESystem(eqs, independent_variable(eqs_er3bp), states(eqs_er3bp), [μ])
end

# Build the equations at pre-compile time
const CR3BP_ODEFunctions = _CR3BP_ODEFunctions()

#-------------#
# CR3BP MODEL #
#-------------#
struct CR3BP{handcoded,O<:_CR3BP_ODEFunctions,P<:R3BPSystemProperties} <: Abstract_R3BPModel
    ode   :: O
    props :: P
end

# TODO: Move hand-coded functions into :analytical of CR3BP
CR3BP(args...; kwargs...) = CR3BP{Val{false}}(R3BPSystemProperties(args...; kwargs...))
function CR3BP{nothandcoded}(props::R3BPSystemProperties; kwargs...) where {nothandcoded<:Val{false}}
    if length(kwargs) == 0
        ode = CR3BP_ODEFunctions
    else
        ode = _CR3BP_ODEFunctions(; kwargs...)
    end
    CR3BP{false, typeof(ode), typeof(props)}(ode, props)
end

# XXX: The parameters overload below gives us a performance boost.
ModelingToolkit.parameters(model::CR3BP{false}) = SVector(model.props.μ)

# Hand-coded CR3BP
HandcodedCR3BP(args...; kwargs...) = CR3BP{Val{true}}(R3BPSystemProperties(args...; kwargs...))
function CR3BP{handcoded}(props::R3BPSystemProperties) where {handcoded<:Val{true}}
    ode = _CR3BP_ODEFunctions(nothing, handcodedCR3BP, handcodedCR3BP_withSTM)
    CR3BP{true, typeof(ode), typeof(props)}(ode, props)
end
ModelingToolkit.parameters(model::CR3BP{true}) = SVector(model.props.μ)

#---------#
# METHODS #
#---------#
@doc "Jacobi integral = 2U - V^2 = -2E [Oshima 2019, Eq. 3; Koon, Eq. 2.3.13; Wakker, Eq. 3.54]"
jacobi(μ, vel, pos) = 2*centrifugal_potential(μ, pos) - norm(vel)^2

function handcodedCR3BP(du, u, p, t)
    @inbounds begin
        μ = p[1]
        x, y, z, x̂, ŷ, ẑ = u
        r1, r2 = (1 - μ, μ) ./ (distance_vectors(μ, (x, y, z)).^3)
        du[1] = x̂
        du[2] = ŷ
        du[3] = ẑ
        du[4] =  2ŷ  + x - r1*(μ + x) + r2*(1 - μ - x)
        du[5] = -2x̂  + y - r1*y       - r2*y
        du[6] =          - r1*z       - r2*z
    end
end

function handcodedCR3BP_withSTM(du, u, p, t)
    @inbounds begin
        μ = p[1]
        x, y, z, x̂, ŷ, ẑ = u
        (d, r) = distance_vectors(μ, (x, y, z))
        d3, d5 = d^3, d^5
        r3, r5 = r^3, r^5
        r1, r2 = (1 - μ, μ) ./ ((d, r).^3)
        A, B = 3*(1-μ)/d5, 3*μ/r5

        du[1] = x̂
        du[2] = ŷ
        du[3] = ẑ
        du[4] =  2ŷ  + x - r1*(μ + x) + r2*(1 - μ - x)
        du[5] = -2x̂  + y - r1*y       - r2*y
        du[6] =          - r1*z       - r2*z

        XX = 1 - (1-μ)/d3 - μ/r3 + A*(x+μ)^2 + B*(x - 1 + μ)^2
        XY = A*(x+μ)*y + B*(x-1+μ)*y
        XZ = A*(x+μ)*z + B*(x-1+μ)*z
        YY = 1 - (1-μ)/d3 - μ/r3 + A*y^2 + B*y^2
        YZ = A*y*z + B*y*z
        ZZ = -(1-μ)/d3 - μ/r3 + A*z^2 + B*z^2

        du[7] = u[10]
        du[8] = u[11]
        du[9] = u[12]
        du[10] = XX*u[7] + XY*u[8] + XZ*u[9] + 2*u[11]
        du[11] = XY*u[7] + YY*u[8] + YZ*u[9] - 2*u[10]
        du[12] = XZ*u[7] + YZ*u[8] + ZZ*u[9]

        du[13] = u[16]
        du[14] = u[17]
        du[15] = u[18]
        du[16] = XX*u[13] + XY*u[14] + XZ*u[15] + 2*u[17]
        du[17] = XY*u[13] + YY*u[14] + YZ*u[15] - 2*u[16]
        du[18] = XZ*u[13] + YZ*u[14] + ZZ*u[15]

        du[19] = u[22]
        du[20] = u[23]
        du[21] = u[24]
        du[22] = XX*u[19] + XY*u[20] + XZ*u[21] + 2*u[23]
        du[23] = XY*u[19] + YY*u[20] + YZ*u[21] - 2*u[22]
        du[24] = XZ*u[19] + YZ*u[20] + ZZ*u[21]

        du[25] = u[28]
        du[26] = u[29]
        du[27] = u[30]
        du[28] = XX*u[25] + XY*u[26] + XZ*u[27] + 2*u[29]
        du[29] = XY*u[25] + YY*u[26] + YZ*u[27] - 2*u[28]
        du[30] = XZ*u[25] + YZ*u[26] + ZZ*u[27]

        du[31] = u[34]
        du[32] = u[35]
        du[33] = u[36]
        du[34] = XX*u[31] + XY*u[32] + XZ*u[33] + 2*u[35]
        du[35] = XY*u[31] + YY*u[32] + YZ*u[33] - 2*u[34]
        du[36] = XZ*u[31] + YZ*u[32] + ZZ*u[33]

        du[37] = u[40]
        du[38] = u[41]
        du[39] = u[42]
        du[40] = XX*u[37] + XY*u[38] + XZ*u[39] + 2*u[41]
        du[41] = XY*u[37] + YY*u[38] + YZ*u[39] - 2*u[40]
        du[42] = XZ*u[37] + YZ*u[38] + ZZ*u[39]
    end
end


text:
using Test, Documenter, Quaternionic

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


text:
using Test

n = 10;
A = rand(n, n); A = A + A';

include("ldlt_symm.jl")

# L = Array(Float64, 2)
# D = Array(Float64, 2)
# P = Array(Float64, 2)
# pp = Array(Int, 1)
global L     = Matrix{Float64}(undef, n, n)
global D     = Matrix{Float64}(undef, n, n)
global P     = Matrix{Float64}(undef, n, n)
global pp    = Vector{Float64}(undef, n)
global ρ     = Float64
global ncomp = Int64

try (L, D, pp, ρ, ncomp) = ldlt_symm(A)
catch except
    @test rethrow(except)
end

#@test_approx_eq P*A*P'  L*D*L'
# @test_approx_eq A[pp,pp]  L*D*L'
@test A[pp, pp] ≈ L * D * L'

v = rand(n,1);
#@test_approx_eq P*v  v[pp]
#@test_approx_eq P'*v  v[invperm(pp)]


text:
@testset "Katz" begin
    g5 = DiGraph(4)
    add_edge!(g5, 1, 2); add_edge!(g5, 2, 3); add_edge!(g5, 1, 3); add_edge!(g5, 3, 4)
    for g in testdigraphs(g5)
        z = @inferred(katz_centrality(g, 0.4))
        @test round.(z, 2) == [0.32, 0.44, 0.62, 0.56]
    end
end


text:
function Initi_model_struct(H::T,Hc::T,N::T,Wf::T,B::T,αf::T,β₁::T,β₂::T,
    Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    cons = Constants()
    env = EnvironmentStruct()
    kinetic = PhotoKineticRates()
    photo = PhotoPar(kinetic,env.Tₐ)
    hydPar = HydraulicsPar(;Kₓₗ₀=Kₓₗ₀,i=i)
    treepar = TreePar(;αf=αf,β₁=β₁,β₂=β₂,rₘ=rₘ,Nₛ=Nₛ,a_Jmax=a_Jmax,b_Jmax=b_Jmax)

    Hs = H-Hc
    As = Wf/treepar.αf
    Ww = treepar.ρw*As*(treepar.β₁*H+treepar.β₂*Hs)
    treesize = TreeSize(Wf,Ww,H,Hs,As,B,N)

    model = CCPHStruct(cons,env,treepar,treesize,photo,hydPar)

    return model,kinetic
end

function Initi_model_struct(j::Integer,data::RoData,
    αf::T,β₁::T,β₂::T,Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    data_ind = Find_data_ind(j)
    H = data.H[data_ind]
    N = data.N[data_ind]
    Wf = data.Wf[data_ind]
    B = data.B[data_ind]
    Hc = data.Hc
    model,kinetic = Initi_model_struct(H,Hc,N,Wf,B,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

    return model,kinetic
end

function Get_Data_RO_C_F_CCPH(par::Array{Float64,1},weatherts_F::WeatherTS,
    weatherts_C::WeatherTS,data_F::RoData,data_C::RoData;sim_steps::Integer=83)

    αf,β₁,β₂,b_Jmax = 460.0,1.27,-0.27,0.0
    Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i = par[1:6]

    #Fertilized stand
    GPP_model_F = zeros(sim_steps+1)
    EC_model_F = zeros(sim_steps+1)
    gₛ_model_F = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_F,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_F[j],EC_model_F[j], modeloutput, gₛ_model_F[j], Nₘ_f_opt = CalcModelOutput!(j,model,weatherts_F,kinetic)
    end

    #Control stand
    GPP_model_C = zeros(sim_steps+1)
    EC_model_C = zeros(sim_steps+1)
    gₛ_model_C = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_C,αf,β₁,β₂,Nₛ_C,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_C[j],EC_model_C[j], modeloutput, gₛ_model_C[j], Nₘ_f_opt = CalcModelOutput!(j,model,weatherts_C,kinetic)
    end

    return GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,gₛ_model_F,gₛ_model_C
end

function run_validation_RO_2019()
    file_name = "adaptive_rwm_RO_C_GPP_Ec_20220228"
    #Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i,X0,τ,Smax,τ_C,a_GPP,b_GPP,a_Ec,b_Ec
    par =  [0.0028160038624083067, 27.33276184313149,
     0.028643648082103606, 0.09336572808750652,
      0.0019727630599060625, 2.2083338895077254,
       -6.817991596443347, 13.184687442117964,
        20.387232607829905, 2.6564297126274803,
         0.434540406001315, 0.06910534451564838,
          0.13628204854679268, 0.1383012394158192]

    weather_file="Weather_RO_2"

    RO_data = Load_RO_data(weather_file=weather_file)

    X0,τ,Smax,τ_C = par[7:10]
    weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax,weather_file= weather_file)
    weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax,weather_file= weather_file)

    data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

    n_data = length(weatherts_F.date)

    GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,gₛ_model_F,gₛ_model_C = Get_Data_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C;sim_steps=n_data-1)

    ind = 85:n_data

    @show GPP_R²_F =  calcR²(data_F.GPP[ind],GPP_model_F[ind])
    @show GPP_cor_F = cor(data_F.GPP[ind],GPP_model_F[ind])
    @show gₛ_R²_F =  calcR²(data_F.E_C[ind],EC_model_F[ind])
    @show gₛ_cor_F = cor(data_F.E_C[ind],EC_model_F[ind])

    @show GPP_R²_C =  calcR²(data_C.GPP[ind],GPP_model_C[ind])
    @show GPP_cor_C = cor(data_C.GPP[ind],GPP_model_C[ind])
    @show gₛ_R²_C =  calcR²(data_C.E_C[ind],EC_model_C[ind])
    @show gₛ_cor_C = cor(data_C.E_C[ind],EC_model_C[ind])

    plot(weatherts_F.date,data_F.GPP,label="Data",ylabel="GPP")
    pl1 = plot!(weatherts_F.date,GPP_model_F,label="Model")
    plot(weatherts_C.date,data_C.GPP,label="Data",ylabel="GPP")
    pl2 = plot!(weatherts_C.date,GPP_model_C,label="Model")

    plot(weatherts_F.date,data_F.E_C,label="Data",ylabel="E_C")
    pl3 = plot!(weatherts_F.date,EC_model_F,label="Model")
    plot(weatherts_C.date,data_C.E_C,label="Data",ylabel="E_C")
    pl4 = plot!(weatherts_C.date,EC_model_C,label="Model")

    plot(pl1,pl3,pl2,pl4,layout=(2,2),legends=false)

    a_GPP,b_GPP,a_Ec,b_Ec = par[11:14]

    σ_GPP_F = a_GPP.+GPP_model_F*b_GPP
    σ_GPP_C = a_GPP.+GPP_model_C*b_GPP

    σ_EC_F = a_Ec.+b_Ec*EC_model_F
    σ_EC_C = a_Ec.+b_Ec*EC_model_C

    c = -log(0.025*2) #Use for calculating 95% credible intervals

    plot(weatherts_F.date[ind],data_F.GPP[ind],label="Data",ylabel="GPP")
    plot!(weatherts_F.date[ind],GPP_model_F[ind]+c*σ_GPP_F[ind])
    plot!(weatherts_F.date[ind],GPP_model_F[ind]-c*σ_GPP_F[ind])
    pl1 = plot!(weatherts_F.date[ind],GPP_model_F[ind],label="Model")
    plot(weatherts_C.date[ind],data_C.GPP[ind],label="Data",ylabel="GPP")
    plot!(weatherts_C.date[ind],GPP_model_C[ind]+c*σ_GPP_C[ind])
    plot!(weatherts_C.date[ind],GPP_model_C[ind]-c*σ_GPP_C[ind])
    pl2 = plot!(weatherts_C.date[ind],GPP_model_C[ind],label="Model")

    plot(weatherts_F.date[ind],data_F.E_C[ind],label="Data",ylabel="E_C")
    plot!(weatherts_F.date[ind],EC_model_F[ind]+c*σ_EC_F[ind])
    plot!(weatherts_F.date[ind],EC_model_F[ind]-c*σ_EC_F[ind])
    pl3 = plot!(weatherts_F.date[ind],EC_model_F[ind],label="Model")
    plot(weatherts_C.date[ind],data_C.E_C[ind],label="Data",ylabel="E_C")
    plot!(weatherts_C.date[ind],EC_model_C[ind]+c*σ_EC_C[ind])
    plot!(weatherts_C.date[ind],EC_model_C[ind]-c*σ_EC_C[ind])
    pl4 = plot!(weatherts_C.date[ind],EC_model_C[ind],label="Model")

    plot(pl1,pl2,pl3,pl4,layout=(2,2),legends=false)
    savefig("./plots/"*file_name*"_validation.svg")
end

run_validation_RO_2019()


text:
##
using Loom
using HDF5

##
file = h5open("test.h5", "w")
typeof(file)
file["matrix"] = randn(10, 10)
file["row_attrs"] = randn(10, 2)
file["col_attrs"] = randn(10, 3)
close(file)

##
lmfile = Loom.open("test.h5"; mode="r")
##
lmfile.matrix

##
Loom.matrixsize(lmfile)

##
show(lmfile.conn.ra)

typeof(lmfile.matrix)

##
close(lmfile)

##
lmfile[1:2, 1:2]


text:
using LibCURL

# init a curl handle
curl = curl_easy_init();

# set the URL and request to follow redirects
curl_easy_setopt(curl, CURLOPT_URL, "http://www.amisllp.uk")
curl_easy_setopt(curl, CURLOPT_FOLLOWLOCATION, 1);

# setup the callback function to recv data
function curl_write_cb(curlbuf::Ptr{Nothing}, s::Csize_t, n::Csize_t, p_ctxt::Ptr{Nothing})
    sz = s * n
    data = Array{UInt8,1}(undef,sz)
    ccall(:memcpy, Ptr{Nothing}, (Ptr{Nothing}, Ptr{Nothing}, UInt64), data, curlbuf, sz)
    println("recd: ", String(data))
    sz::Csize_t
end

c_curl_write_cb = @cfunction curl_write_cb  Csize_t (Ptr{Nothing}, Csize_t, Csize_t, Ptr{Nothing});
curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, c_curl_write_cb);

# execute the query
res = curl_easy_perform(curl);
println("curl url exec response : ", res);

# retrieve HTTP code
http_code = Array{Clong,1}(undef,1);
curl_easy_getinfo(curl, CURLINFO_RESPONSE_CODE, http_code);
println("httpcode : ", http_code);

# release handle
curl_easy_cleanup(curl);


text:
module ExtendableSparse
using SparseArrays
using LinearAlgebra
using SuiteSparse
using Requires

using DocStringExtensions

import SparseArrays: rowvals, getcolptr, nonzeros
import Base: copy

include("sparsematrixcsc.jl")
include("sparsematrixlnk.jl")
include("extendable.jl")

export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, updateindex!, rawupdateindex!, colptrs

include("factorizations.jl")
export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner, ParallelILU0Preconditioner, reorderlinsys
export AbstractFactorization,LUFactorization, CholeskyFactorization
export issolver
export factorize!, update!
export ILUTPreconditioner, AMGPreconditioner
export PardisoLU, MKLPardisoLU

include("simple_iteration.jl")
export simple,simple!

include("sprand.jl")
export sprand!,sprand_sdd!, fdrand,fdrand!,fdrand_coo

function __init__()
    @require Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" include("pardiso_lu.jl")
    @require IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" include("ilut.jl")
    @require AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" include("amg.jl")
end

end # module


text:
export TD, FD
import Base: unsafe_convert
using Base: size_to_strides
using Knet.KnetArrays: DevArray
using AutoGrad: Value

using CUDA.CUDNN: handle,
    cudnnTensorDescriptor_t,
    cudnnCreateTensorDescriptor,
    cudnnSetTensor4dDescriptor,
    cudnnSetTensor4dDescriptorEx,
    cudnnGetTensor4dDescriptor,
    cudnnSetTensorNdDescriptor,
    cudnnSetTensorNdDescriptorEx,
    cudnnGetTensorNdDescriptor,
    cudnnGetTensorSizeInBytes,
    cudnnDestroyTensorDescriptor,

    cudnnFilterDescriptor_t,
    cudnnCreateFilterDescriptor,
    cudnnSetFilter4dDescriptor,
    cudnnGetFilter4dDescriptor,
    cudnnSetFilterNdDescriptor,
    cudnnGetFilterNdDescriptor,
    cudnnGetFilterSizeInBytes,
    cudnnDestroyFilterDescriptor,

    cudnnDataType_t,
        CUDNN_DATA_FLOAT,
        CUDNN_DATA_DOUBLE,
        CUDNN_DATA_HALF,
        CUDNN_DATA_INT8,
        CUDNN_DATA_INT32,
        CUDNN_DATA_INT8x4,
        CUDNN_DATA_UINT8,
        CUDNN_DATA_UINT8x4,
        CUDNN_DATA_INT8x32,

    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,
        CUDNN_TENSOR_NHWC,
        CUDNN_TENSOR_NCHW_VECT_C

# cudnn tensor/filter descriptors: need to use mutable structs in order to have finalizers.

mutable struct TensorDescriptor; ptr::cudnnTensorDescriptor_t; end
unsafe_convert(::Type{<:Ptr}, td::TensorDescriptor)=td.ptr
TD(a) = TD(eltype(a),size(a))
function TD(T::Type, dims::Dims{N}) where {N}
    ptr = cudnnTensorDescriptor_t[C_NULL]
    cudnnCreateTensorDescriptor(ptr)
    sz = Cint[reverse(dims)...]
    st = Cint[reverse(size_to_strides(1,dims...))...]
    cudnnSetTensorNdDescriptor(ptr[1], DT(T), N, sz, st)
    td = TensorDescriptor(ptr[1])
    finalizer(x->cudnnDestroyTensorDescriptor(x.ptr), td)
    return td
end

mutable struct FilterDescriptor; ptr::cudnnFilterDescriptor_t; end
unsafe_convert(::Type{<:Ptr}, fd::FilterDescriptor)=fd.ptr
FD(a) = FD(eltype(a),size(a),CUDNN_TENSOR_NCHW)
function FD(T::Type, dims::Dims{N},format::cudnnTensorFormat_t) where N
    ptr = cudnnFilterDescriptor_t[C_NULL]
    cudnnCreateFilterDescriptor(ptr)
    sz = Cint[reverse(dims)...]
    cudnnSetFilterNdDescriptor(ptr[1], DT(T), format, N, sz)
    fd = FilterDescriptor(ptr[1])
    finalizer(x->cudnnDestroyFilterDescriptor(x.ptr), fd)
    return fd
end

DT(::Type{T}) where T = get(DataTypes, T) do; error("CUDNN does not support $T"); end
DataTypes = Dict{Type,cudnnDataType_t}(
    Float32 => CUDNN_DATA_FLOAT,
    Float64 => CUDNN_DATA_DOUBLE,
    Float16 => CUDNN_DATA_HALF,
    Int8 => CUDNN_DATA_INT8,
    UInt8 => CUDNN_DATA_UINT8,
    Int32 => CUDNN_DATA_INT32,
    # The following are 32-bit elements each composed of 4 8-bit integers, only supported with CUDNN_TENSOR_NCHW_VECT_C (TODO)
    # CUDNN_DATA_INT8x4,
    # CUDNN_DATA_UINT8x4,
    # CUDNN_DATA_INT8x32,
)

nothing


text:
module StaticSites

using Dates
using Logging
using ..Remark
export StaticSite, generate_page, prepare

struct StaticSite
    source_root     :: String
    target_root     :: String
    static_root     :: String
    default_modules :: Vector{Module}
    default_page    :: String

    function StaticSite(;
                        source_root="remark/",
                        target_root="target/",
                        static_root="static/",
                        default_modules=[],
                        default_page="main.rem")
        new(source_root,
            target_root,
            static_root,
            default_modules,
            default_page)
    end
end

"""
Prepare the site for page generation by:

 1. Building a directory structure skeleton
 2. Copying static files
"""
function prepare(site::StaticSite)
    mkpath(site.target_root)
    for file in readdir(site.static_root)
        @info "Copying static file" file
        cp(joinpath(site.static_root, file),
           joinpath(site.target_root, file);
           force=true)
    end
end

function generate_page(site::StaticSite,
                       destination, page=site.default_page;
                       data::AbstractDict{Symbol}=Dict{Symbol}(),
                       modules=[])
    start = Dates.now()
    description = isempty(destination) ? "Index page" : destination
    @debug "Started generating remark page" start page=description

    location = joinpath(site.target_root, destination)
    mkpath(location)

    source = joinpath(site.source_root, page)

    data[:currentpage] = destination
    open(joinpath(location, "index.html"), "w") do f
        Remark.tohtml(f, source, data;
                                 modules=vcat(site.default_modules, modules))
        println(f)
    end

    finish = Dates.now()
    @info "Generated remark page" elapsed=(finish - start) page=description
end

end


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Hinton diagram generation using vector graphics.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-12-24

# TODO: Consider a version with JavaScript injection.

using Colors
using Compose

immutable HintonVecStyle
    pos   ::Color
    neg   ::Color
    bg    ::Color
    margin::Float64
end
pos(o::HintonVecStyle)    = o.pos
neg(o::HintonVecStyle)    = o.neg
bg(o::HintonVecStyle)     = o.bg
margin(o::HintonVecStyle) = o.margin
const VECSTYLE = HintonVecStyle(
    colorant"white",
    RGB(map(x -> x/255, (51 , 51 , 50 ))...), # Dark grey, "almost" black.
    RGB(map(x -> x/255, (132, 132, 130))...), # Battleship grey.
    0.002,
    )

function hintonvec{T<:Real}(w::AbstractMatrix{T}, sty::HintonVecStyle=VECSTYLE)
    mabs   = maximum(abs.(w))
    (m, n) = size(w)
    mgn    = margin(sty)
    step   = (1 - (max(m, n) + 1)*mgn)/max(m, n)
    cont   = compose(context(),
        rectangle(0, 0, (n + 1)*mgn + n*step, (m + 1)*mgn + m*step),
        fill(bg(sty)),
        )
    for i in 1:n, j in 1:m
        v      = w[j, i]/mabs
        colour = v < 0.0 ? neg(sty) : pos(sty)
        cont   = compose(cont, compose(
            context(i*mgn + (i - 1)step, j*mgn + (j - 1)step, step, step),
            rectangle((1 - v)/2, (1 - v)/2, v, v),
            fill(colour),
            ))
    end
    cont
end


text:
# Implements the legendre interface for the unit normalization case

abstract type AbstractLegendreSphereNorm <: AbstractLegendreNorm end

"""
    struct LegendreFourPiNorm <: AbstractLegendreNorm end

Trait type denoting the ``4\\pi`` (geodesy) normalization of the associated Legendre
functions ``F_\\ell^m(x)``.
The orthonormal normalization is defined such that
```math
    \\int_{-1}^{1} \\left[ F_\\ell^m(x) \\right]^2 \\,dx = 2
```
The normalization factor with respect to the standard (unnormalized) Legendre
functions ``P_\\ell^m(x)`` ([`LegendreUnitNorm`](@ref)) is given by
```math
    F_\\ell^m(x) \\equiv \\sqrt{2\\pi(2\\ell+1) \\frac{(\\ell-m)!}{(\\ell+m)!}} P_\\ell^m(x)
```
"""
struct LegendreFourPiNorm <: AbstractLegendreSphereNorm end

"""
    struct LegendreOrthoNorm <: AbstractLegendreNorm end

Trait type denoting the orthonormal (full) normalization of the associated Legendre
functions ``O_\\ell^m(x)``.
The orthonormal normalization is defined such that
```math
    \\int_{-1}^{1} \\left[ O_\\ell^m(x) \\right]^2 \\,dx = 1
```
The normalization factor with respect to the standard (unnormalized) Legendre
functions ``P_\\ell^m(x)`` ([`LegendreUnitNorm`](@ref)) is given by
```math
    O_\\ell^m(x) \\equiv \\sqrt{\\frac{2\\ell+1}{2} \\frac{(\\ell-m)!}{(\\ell+m)!}} P_\\ell^m(x)
```
"""
struct LegendreOrthoNorm <: AbstractLegendreSphereNorm end

"""
    struct LegendreSphereNorm <: AbstractLegendreNorm end

Trait type denoting the spherical-harmonic normalization of the associated Legendre
functions ``\\lambda_\\ell^m(x)``.
The spherical-harmonic normalization is defined such that
```math
    \\int_{-1}^{1} \\left[ \\lambda_\\ell^m(x) \\right]^2 \\,dx = \\frac{1}{2\\pi}
```
The normalization factor with respect to the standard (unnormalized) Legendre
functions ``P_\\ell^m(x)`` ([`LegendreUnitNorm`](@ref)) is given by
```math
    \\lambda_\\ell^m(x) \\equiv \\sqrt{\\frac{2\\ell+1}{4\\pi} \\frac{(\\ell-m)!}{(\\ell+m)!}} P_\\ell^m(x)
```
"""
struct LegendreSphereNorm <: AbstractLegendreSphereNorm end

@inline function initcond(::LegendreFourPiNorm, ::Type{T}) where T
    return one(T)
end
@inline function initcond(::LegendreOrthoNorm, ::Type{T}) where T
    return sqrt(inv(T(2)))
end
@inline function initcond(::LegendreSphereNorm, ::Type{T}) where T
    # comparing this against
    #   convert(T, inv(sqrt(4*convert(BigFloat, π))))
    # shows that this is exact within Float64 precision
    return inv(sqrt(4 * convert(T, π)))
end

# Version of sqrt() which skips the domain (x < 0) check for the IEEE floating point types.
# For nonstandard number types, just fall back to a regular sqrt() since eliminating the
# domain check is probably no longer the dominant contributor to not vectorizing.
unchecked_sqrt(x::T) where {T <: Base.IEEEFloat} = Base.sqrt_llvm(x)
unchecked_sqrt(x::T) where {T <: Integer} = unchecked_sqrt(float(x))
unchecked_sqrt(x) = Base.sqrt(x)

@inline function
coeff_μ(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer) where T
    # The direct derivation of the normalization constant gives
    #     return sqrt(one(T) + inv(convert(T, 2l)))
    # but when comparing results for T ∈ (Float64,BigFloat), the Float64 results differ by
    # ±1 ulp at various points.
    #
    # Instead, consider the following which limits to 0 (i.e. as (2l)^-1) rather than 1.
    # ```math
    #     μ'_ℓ = μ_ℓ - 1 = \sqrt{1 + \frac{1}{2ℓ}} - 1
    # ```
    # Now complete the square by multiplying through by a factor of unity to instead write
    # ```math
    #   \left( \sqrt{1 + \frac{1}{2ℓ}} - 1 \right) \times
    #       \frac{\sqrt{1 + \frac{1}{2ℓ}} + 1}{\sqrt{1 + \frac{1}{2ℓ}} + 1}
    #   = \frac{1}{\sqrt{2ℓ(2ℓ + 1) + 2ℓ}}
    # ```
    # Then rewrite ``μ_ℓ`` in terms of ``μ'_ℓ`` to arrive at
    # ```math
    #   μ_ℓ = 1 + \left[ 2ℓ + \sqrt{2ℓ(2ℓ + 1)} \right]^{-1}
    # ```
    # The square root is calculated on a growing quantity rather than one limiting to 1,
    # so the significance of the low-order digits are better preserved. This expression
    # allows Float64 calculations to be correctly rounded values when compared to BigFloat
    # to very high ℓ.
    xT = convert(T, 2l)
    return one(T) + inv(xT + unchecked_sqrt(muladd(xT, xT, xT)))
end

@inline function
coeff_ν(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer) where T
    return unchecked_sqrt(convert(T, 2l + 1))
end

@inline function
coeff_α(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer, m::Integer) where T
    lT = convert(T, l)
    mT = convert(T, m)
    # Write factors in two pieces to make compiler's job easier. In the case where
    # both coeff_α and coeff_β are called and inlined, the next line from both functions
    # should be merged and shared.
    fac1 = (2lT + 1) / ((2lT - 3) * (lT^2 - mT^2))
    fac2 = 4*(lT - 1)^2 - 1
    return unchecked_sqrt(fac1 * fac2)
end

@inline function
coeff_β(::AbstractLegendreSphereNorm, ::Type{T}, l::Integer, m::Integer) where T
    lT = convert(T, l)
    mT = convert(T, m)
    # Write factors in two pieces to make compiler's job easier. In the case where
    # both coeff_α and coeff_β are called and inlined, the next line from both functions
    # should be merged and shared.
    fac1 = (2lT + 1) / ((2lT - 3) * (lT^2 - mT^2))
    fac2 = (lT - 1)^2 - mT^2
    return unchecked_sqrt(fac1 * fac2)
end

# Extra functions

"""
    N = Nlm([T=Float64], l, m)

Computes the normalization constant
```math
    N_\\ell^m \\equiv \\sqrt{\\frac{2\\ell+1}{4\\pi} \\frac{(\\ell-m)!}{(\\ell+m)!}}
```
which defines the Spherical Harmonic normalized functions ``\\lambda_\\ell^m(x)`` in
terms of the standard unit normalized ``P_\\ell^m(x)``
```math
    \\lambda_\\ell^m(x) \\equiv N_\\ell^m P_\\ell^m(x)
```
using numbers of type `T`.

See also [`Plm`](@ref) and [`λlm`](@ref).
"""
function Nlm(::Type{T}, l::Integer, m::Integer) where T
    fac1 = one(T)
    for ii in (l-m+1):(l+m)
        fac1 *= convert(T, ii)
    end
    num = 2*convert(T, l) + 1
    den = 4*convert(T, π)
    return sqrt( num * inv(den) * inv(fac1) )
end
Nlm(l::Integer, m::Integer) = Nlm(Float64, l, m)


text:
"""
`LoaderError` should be thrown when loader library code fails, and other libraries should
be given the chance to recover from the error.  Reports the library name and an error message:
LoaderError("ImageMagick", "Foo not available")
"""
struct LoaderError <: Exception
    library::String
    msg::String
    ex
end
Base.showerror(io::IO, e::LoaderError) = println(io, e.library, " load error: ",
                                                 e.msg, "\n  due to ", e.ex, "\n  Will try next loader.")

"""
`WriterError` should be thrown when writer library code fails, and other libraries should
be given the chance to recover from the error.  Reports the library name and an error message:
WriterError("ImageMagick", "Foo not available")
"""
struct WriterError <: Exception
    library::String
    msg::String
    ex
end
Base.showerror(io::IO, e::WriterError) = println(
    io, e.library, " writer error: ",
    e.msg, "\n  due to ", e.ex, "\n  Will try next loader."
)

struct SpecError <: Exception
    mod::Module
    call::Symbol
end
Base.showerror(io::IO, e::SpecError) = print(io, e.mod, " is missing $(e.call) and fileio_$(e.call)")

"""
Handles a list of thrown errors after no IO library was found working
"""
function handle_exceptions(exceptions::Vector, action)
    # first show all errors when there are more then one
    multiple = length(exceptions) > 1
    println(stderr, "Error$(multiple ? "s" : "") encountered while $action.")
    if multiple
        println("All errors:")
        println("===========================================")
        for (err, file, bt) in exceptions
            showerror(stdout, err)
            println("\n===========================================")
        end
    end
    # then handle all errors.
    # this way first fatal exception throws and user can still see all errors
    # TODO, don't throw, if it contains a NotInstalledError?!
    println(stderr, "\nFatal error:")
    for exception in exceptions
        continue_ = handle_error(exception...)
        continue_ || break
    end
end

handle_error(e, q, bt) = throw(CapturedException(e, stacktrace(bt)))


text:
module HydrostaticBoussinesq

export HydrostaticBoussinesqModel, Forcing

using StaticArrays
using LinearAlgebra: dot, Diagonal
using CLIMAParameters.Planet: grav

using ..Ocean
using ...VariableTemplates
using ...MPIStateArrays
using ...Mesh.Filters: apply!
using ...Mesh.Grids: VerticalDirection
using ...Mesh.Geometry
using ...DGMethods
using ...DGMethods: init_state_auxiliary!
using ...DGMethods.NumericalFluxes
using ...DGMethods.NumericalFluxes: RusanovNumericalFlux
using ...BalanceLaws

import ..Ocean: coriolis_parameter
import ...DGMethods.NumericalFluxes: update_penalty!

import ...BalanceLaws:
    vars_state,
    init_state_prognostic!,
    init_state_auxiliary!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_state!,
    boundary_conditions,
    update_auxiliary_state!,
    update_auxiliary_state_gradient!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

import ..Ocean:
    ocean_init_state!,
    ocean_init_aux!,
    ocean_boundary_state!,
    _ocean_boundary_state!

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)
⊗(a::SVector, b::SVector) = a * b'

include("hydrostatic_boussinesq_model.jl")
include("bc_velocity.jl")
include("bc_temperature.jl")
include("LinearHBModel.jl")
include("Courant.jl")

end


text:
struct Rectangle
    x::Float64
    y::Float64
    w::Float64
    h::Float64
end

Rectangle() = Rectangle(0.0, 0.0, 0.0, 0.0)

function contains(r1::Rectangle, r2::Rectangle)
    r1.x <= r2.x && r1.x+r1.w >= r2.x+r2.w && r1.y <= r2.y && r1.y+r1.h >= r2.y+r2.h
end

function overlaps(r1::Rectangle, r2::Rectangle)
    r1.x < r2.x+r2.w && r2.x < r1.x+r1.w && r1.y < r2.y+r2.h && r2.y < r1.y+r1.h
end

function merge(rects::Vector{Rectangle})
    minx, miny = rects[1].x, rects[1].y
    maxx, maxy = rects[1].x+rects[1].w, rects[1].y+rects[1].h
    for r in rects
        r.x < minx && (minx = r.x)
        r.y < miny && (miny = r.y)
        r.x+r.w > maxx && (maxx = r.x+r.w)
        r.y+r.h > maxy && (maxy = r.y+r.h)
    end
    Rectangle(minx, miny, maxx-minx, maxy-miny)
end
merge(rects::Rectangle...) = merge([rects...])

function Base.parse(::Type{Rectangle}, s::String)
    items = Vector{String}(split(s," "))
    @assert length(items) == 4
    xywh = map(x -> parse(Float64,x), items)
    Rectangle(xywh...)
end

Base.string(r::Rectangle) = join([r.x,r.y,r.w,r.h], " ")


text:
# ---
# title: 74. Search a 2D Matrix
# id: problem74
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/search-a-2d-matrix/description/>
# hidden: true
# ---
#
# Write an efficient algorithm that searches for a value in an `m x n` matrix.
# This matrix has the following properties:
#
#   * Integers in each row are sorted from left to right.
#   * The first integer of each row is greater than the last integer of the previous row.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/05/mat.jpg)
#
#
#
#     Input: matrix = [[1,3,5,7],[10,11,16,20],[23,30,34,50]], target = 3
#     Output: true
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/05/mat2.jpg)
#
#
#
#     Input: matrix = [[1,3,5,7],[10,11,16,20],[23,30,34,50]], target = 13
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: matrix = [], target = 0
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[i].length`
#   * `0 <= m, n <= 100`
#   * `-104 <= matrix[i][j], target <= 104`
#
#
## @lc code=start
using LeetCode

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


text:
include("euler/euler.jl")
using .Calculus: champernowne_digits
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    result = 1
    for i ∈ 0:n
        j, k = 10 ^ i, 0
        while champernowne_digits(k + 1) < j
            k += 1
        end
        k += 1
        remainder = j - champernowne_digits(k - 1)
        number = 10 ^ (k - 1) + remainder ÷ k - 1
        if remainder % k == 0
            result *= parse(Int, string(number)[end])
        else
            result *= parse(Int, string(number + 1)[remainder % k])
        end
    end
    return result
end

compute(6)

@benchmark compute(6)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using MetabolicFluxAnalysis

"""
  WX_Transform
represents a parametrization in terms of net / xch flux for reversible fluxes.

Usage to resolve indeces of the parametrization:
  get_wxflux_net(net,wxt,fid) : returns the position of the net flux for the
                                WX-pair which contains flux fid
  get_wxflux_xch(net,wxt,fid) : returns the position of the xch flux for the
                                WX-pair which contains flux fid

#NOTE: The indeces of the WX parametrization work as follows:
first, the wx-pair is resolved for a given flux, then the position of the FIRST
flux of that pair is always used to parametrize the net flux, the SECOND one
always to parametrize the XCH flux.
"""
struct WX_Transform
    wparams::Array{Int64,2} # wparams[i ,1/2] is index of ith net / sx flux (parametrization)
                             # by convention, the net flux is at the pos of the first pos. original flux,
                             # and the sx flux is at the second pos. original flux

    f::Any
end

function eval_wx_transform(wt::WX_Transform,xp::Array{Float64,2})
    xf = copy(xp)
    for zi in 1:size(wt.wparams,1)
        wti = wt.wparams[zi,:]
        for zx in 1:size(xp,2)
            (wfa::Float64,wfb::Float64) = wt.f( xp[ wti[1] , zx ] , xp[ wti[2] , zx ] )
            xf[ wti , zx ] = [wfa;wfb]
        end
    end
    return xf
end

"""
  get_wxflux_net(net::Network , wt::WX_Transform , fid::String )
returns the index of the net-flux for the fid flux.

This function searches for the flux fid in the original network, uses this to
resolve the WX-Pair, and returns the first index of the pair.

# NOTE
If the flux is NOT part of a WX-pair, then it returns the flux.
"""
function get_wxflux_net(net::Network , wt::WX_Transform , fid::String )
    idx_flux = find( x -> x==fid , net.fluxes )

    findpair = (wt.wparams .== idx_flux)
    swp      = sum(findpair[:])
    if(swp==1)
        pair = find( any(findpair,2) )[1]
        return wt.wparams[pair,1]
    else
        if(swp==0) # flux not in wx pair, return idx itself
            return idx_flux
        else
            error("Something wrong.. swp = $(swp)")
        end
    end
end

"""
  get_wxflux_xch(net::Network , wt::WX_Transform , fid::String )
returns the index of the net-flux for the fid flux.

This function searches for the flux fid in the original network, uses this to
resolve the WX-Pair, and returns the second index of the pair.

# NOTE
If the flux is NOT part of a WX-pair, then it returns -1
"""
function get_wxflux_xch(net::Network , wt::WX_Transform , fid::String )
    idx_flux = find( x -> x==fid , net.fluxes )

    findpair = (wt.wparams .== idx_flux)
    swp      = sum(findpair[:])
    if(swp==1)
        pair = find( any(findpair,2) )[1]
        return wt.wparams[pair,2]
    else
        if(swp==0) # flux not in wx pair, return idx itself
            return -1
        else
            error("Something wrong.. swp = $(swp)")
        end
    end
end

"""
  s_wx_flux(wxt::WX_Transform, flux::String)
returns whether given flux has wx parametrization or not
"""
function is_wx_flux(net::Network , wt::WX_Transform , fid::String )
    idx_flux = find( x -> x==fid , net.fluxes )
    if( any( wt.wparams[:] .== idx_flux ) )
        return true
    else
        return false
    end
end

"""
  find_wx_flux(wxt::WX_Transform, flux::String)
returns positions (flux[n],flux[x]) .
Note that either of the two original fluxes can be supplied here.
"""
function get_wx_flux_pair(net::Network , wt::WX_Transform , fid::String )
    idx_flux = find( x -> x==fid , net.fluxes )
    if( any( wt.wparams[:] .== idx_flux ) )
        find_pair = find( any( (wt.wparams .== idx_flux) , 2) )
        pair      = wt.wparams[find_pair,:]
        # sanity check
        if(pair[2]<pair[1]); error("Ordering of wx fluxes wrong.."); end
        return (pair[1],pair[2])
    else
        return (-1,-1)
    end
end

"""
  find_wx_flux(wxt::WX_Transform, flux::String)
returns the index of the given wx paramtrization parameter.
Notation is: flux_id[n] , flux_id[x] for net / xchange flux
"""
function find_wx_flux( net::Network, wxt::WX_Transform, flux::String)
    if(flux[end]==']' && flux[end-2]=='[')
        if( flux[end-1]=='n' || flux[end-1]=='N' )
            return get_wxflux_net(net::Network , wxt::WX_Transform , flux[1:end-3] )
        elseif( flux[end-1]=='x' || flux[end-1]=='X' )
            return get_wxflux_xch(net::Network , wxt::WX_Transform , flux[1:end-3] )
        end
    end
    # if no [x] / [n] is specified, return the flux index..
    return find(fi -> fi==flux ,net.fluxes)
end

"""
  create_all_wx_fluxnames(net::Network,wxt::WX_Transform)
creates all wx flux names
Notation is: flux_id[n] , flux_id[x] for net / xchange flux
"""
function create_all_wx_fluxnames(net::Network,wxt::WX_Transform)
    flux_indeces = collect(1:length(net.fluxes))

    test_wx_flux = (x::Int64) -> any( wxt.wparams[:] .== x )

    flux_names   = Vector{String}(length(net.fluxes))

    for zi in 1:length(net.fluxes)
        if(test_wx_flux(zi))
            wxp = get_wx_flux_pair(net,wxt,net.fluxes[zi])
            flux_names[ wxp[1] ] = string( net.fluxes[wxp[1]] , "[n]")
            flux_names[ wxp[2] ] = string( net.fluxes[wxp[1]] , "[x]")
        else
            flux_names[ zi ]     = net.fluxes[zi]
        end
    end
    return flux_names
end

"""
lc2polyhedron(net::Network , wx::WX_Transform , lcs::Vector{LinConstr})
returns Polyhedron
creates polyhedron from the given constraints

# NOTE: Format of  constraints: flux_id[x] is xchange flux, flux_id[n] is net flux
"""
function lc2polyhedron(net::Network , wxt::WX_Transform , lcs::Vector{LinConstr} )
    n = length(net.fluxes)

    # Init the flux resolution dictionary:
    # NOTE: for fluxes which are not wx-parametrized, and which therefore basically do not have
    # the [x]/[n] postfix, we map the [n] version to the non-postfix flux
    d_s::Dict{String,Int64}   = Dict{String,Int64}()
    #for zi=1:length(net.fluxes); d_s[net.fluxes[zi]] = zi; end
    wx_fluxes = create_all_wx_fluxnames(net,wxt)
    for zi in 1:length(wx_fluxes); print(string(zi," --> " , wx_fluxes[zi] ,"\n")) ;d_s[ wx_fluxes[zi] ] = zi; end
    # add the [n] aliases for non wx-parametrized fluxes
    for fi in wx_fluxes
        if(fi[end]!=']'); d_s[string(fi,"[n]")] = d_s[fi]; end
    end

    d_s["1"] = length(wx_fluxes)+1 # to collect the "1" data

    G::Array{ Float64 , 2} = Array{ Float64 }(0,n)
    h::Array{ Float64 , 2} = Array{ Float64 }(0,1)
    A::Array{ Float64 , 2} = Array{ Float64 }(0,n)
    b::Array{ Float64 , 2} = Array{ Float64 }(0,1)

    v_G =  Array{ RowVector{Float64,Array{Float64,1}},1 }()
    v_A =  Array{ RowVector{Float64,Array{Float64,1}},1 }()
    v_h =  Array{ Float64 , 1 }()
    v_b =  Array{ Float64 , 1 }()

    G_s = spzeros(0,n)
    A_s = spzeros(0,n)

    for lci in lcs
        row_lhs = resolve_lc( n+1 , d_s , lci.lhs )
        row_lhs_one = row_lhs[n+1]
        row_lhs = row_lhs[1:n]
        row_rhs = resolve_lc( n+1 , d_s , lci.rhs )
        row_rhs_one = row_rhs[n+1]
        row_rhs = row_rhs[1:n]

        # copy all to rhs, except for "one" part:
        p_row_lhs = row_lhs - 1.0 * row_rhs
        p_row_rhs = row_rhs_one - row_lhs_one

        if(lci.constr=='=')
            #A = [A;p_row_lhs']; b = [b;p_row_rhs]; continue;
            #push!( v_A , p_row_lhs'); push!( v_b, p_row_rhs); continue;
            A_s = vcat( A_s , p_row_lhs' ); push!( v_b, p_row_rhs); continue;
        end
        if(lci.constr=='<')
            #G = [G;p_row_lhs']; h = [h;p_row_rhs]; continue;
            #push!( v_G , p_row_lhs'); push!(v_h,p_row_rhs); continue;
            G_s = vcat( G_s , p_row_lhs' ); push!(v_h,p_row_rhs); continue;
        end
        if(lci.constr=='>')
            #G = [G;-p_row_lhs']; h = [h;-p_row_rhs]; continue;
            #vG = push!(-p_row_lhs'); push!(v_h,-p_row_rhs); continue;
            G_s = vcat( G_s , -p_row_lhs' ); push!(v_h,-p_row_rhs); continue;
        end
    end
    #print("\nDONE!\n")
    #println("\n\n $(G) \n\n $(h) \n\n $(A) \n\n $(b) \n\n ")
    #G = reduce( (x,y) -> [x;y] , v_G )
    #A = reduce( (x,y) -> [x;y] , v_A )
    h = [ h ; v_h]
    b = [ b ; v_b]
    G = full(G_s)
    A = full(A_s)
    return Polyhedron(G,h,A,b)
end

"""
  get_N(net::Network,wxt::WX_Transform)
returns stoichiometric matrix for network with wx transform
"""
function get_N(net::Network,wxt::WX_Transform)
    # create N, then sort out [xch] fluxes..
    N = get_N(net)

    # remove all constraints involing "second" fluxes of wxparams
    for iwx in wxt.wparams[:,2]
        #icr = find( N[:,iwx] .!= 0 )
        #N = N[setdiff(collect(1:size(N,1)),icr),:]
        N[:,iwx] = 0 # this should be ok, because in every constraint of N, the forward flux should always be non-zero, if this flux is non-zero..
    end
    N_all_zero = find( all( N .==0 ,2) )
    N = N[ setdiff(collect(1:size(N,1)),N_all_zero),:]
    return N
end

"""
  get_flux_polyhedron(net::Network , wxt::WX_Transform , constraints::Vector{LinConstr} ; max_val=Inf )
creates the flux polyhedron for network with wx_transform and wx constraints
CONVENTION:
Here, all (net) fluxes are reversible by default
We limit the xch fluxes to the interval [0,1] by default.
"""
function get_flux_polyhedron(net::Network , wxt::WX_Transform , constraints::Vector{LinConstr} ; max_val=Inf )

     m_n = get_N(net,wxt)
     n   = size(m_n,2)
     p_0 = lc2polyhedron(net , wxt , constraints)

     d_s::Dict{String,Int64}   = Dict{String,Int64}()
     for zi=1:length(net.fluxes); d_s[net.fluxes[zi]] = zi; end

     # create [0,1] bounding box for xch fluxes
     ub =  Inf*ones(n)
     lb = -Inf*ones(n)
     lb[wxt.wparams[:,2]] = 0.0
     ub[wxt.wparams[:,2]] = 1.0

     # if max_val, create bouding box. [-max_val,max_val] around everything..
     if(max_val<Inf)
         ub[ub.>max_val] = max_val
         lb[ get_idx_sx_fluxes(net) ] = -max_val
     end

     p_1 =  Polyhedron( p_0 , bounding_box(lb,ub) )
     p_2 =  Polyhedron( p_1 , Polyhedron( zeros(0,n),zeros(0,1),m_n,zeros(size(m_n,1),1) ) )
     return p_2
end

struct optwx_parametrization
    rfd_p::rfd_parametrization
    wx_t::WX_Transform
end

function eval_opt2raw( tf::optwx_parametrization , x::Array{Float64,2})
    x_pwx   = tf.rfd_p.at_opt2raw(x)
    x_flux  = MetabolicFluxAnalysis.eval_wx_transform(tf.wx_t,x_pwx)
end

function get_optimal_wx_parametrization(net::Network,lcs::Vector{LinConstr},max_flux::Float64)
    (p_pwx,wtx)  = MetabolicFluxAnalysis.parametrization_wiechert_exchange( net, lcs, max_flux )

    # round wx_parameterspace polytope:
    rfd_parametrization_pwx = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( p_pwx ; remove_red_constr=1)

    # sample dat polytope unif:
    parametrization = optwx_parametrization(rfd_parametrization_pwx,wtx)
    return parametrization
end

# find flux  pairs with "exchange flux" ( <-> ), and reparametrize them by
# "net-flux" + "exchange flux".
# n: stoichiometric matrix
# sx_param optiosn:
# "log":  is "log4" , means, sxp(0) = 0.0001*max, sxp(1) = max
function parametrization_wiechert_exchange( net::MetabolicFluxAnalysis.Network , lcs::Vector{MetabolicFluxAnalysis.LinConstr} , max_flux::Float64 ; sx_param="log")#::Polyhedron )

    if(max_flux>1e6); print("Set max flux to 1e6.. need finite bound\n"); end
    max_flux = min(max_flux,1e6)

    p_constr = MetabolicFluxAnalysis.lc2polyhedron(net , lcs::Vector{LinConstr} )

    print("Wiechert SX Parametrization:\n")

    # find pairs in stoich. matrix:
    N = MetabolicFluxAnalysis.get_N(net)
    w_pairs = zeros(Int64,0,2)

    # store parametrized fluxes, to avoid problems with multiple reversible fluxes with same stoichiometry
    # e.g. to handle --> fum_a,fum_a_bw,fum_b,fum_b_bw
    parametrized_fluxes = Vector{Int64}()

    for zi in 1:size(N,2)
        for zj in zi+1:size(N,2)
            if( in(zi,parametrized_fluxes) || in(zj,parametrized_fluxes) ) # skip!
                continue;
            end

            if( rank(N[:,[zi,zj]])==1 )
                # and check if is in opposite direction:
                if( sign.(N[:,zi] ) == -sign.(N[:,zj]) )
                    print("pair: ($(zi),$(zj))  -   $(net.fluxes[zi]) , $(net.fluxes[zj]) ")

                    #check if there are other constraints on it:
                    cfree = !any( p_constr.A[:,[zi,zj]][:] .!=0 ) && !any( p_constr.G[:,[zi,zj]][:] .!=0 )
                    if(!cfree)
                        print(" --> skip parametrization, found in constraints..")
                    else
                        w_pairs = [w_pairs;[zi zj]]
                        push!( parametrized_fluxes , zi); push!( parametrized_fluxes , zj)
                        print(" --> parametrized!")
                    end
                    print("\n");
                end
                #w_pairs = [ w_pairs ; [zi zi]]
            end
        end
    end
    if( sx_param=="log" || sx_param=="log4")
        fxch = (xpxch::Float64) ->  exp( xpxch * (log(10000.))  ) * (max_flux/10000.)
    end
    if( sx_param=="log5")
        fxch = (xpxch::Float64) ->  exp( xpxch * (log(100000.))  ) * (max_flux/100000.)
    end
    if( sx_param=="log3")
        fxch = (xpxch::Float64) ->  exp( xpxch * (log(1000.))  ) * (max_flux/1000.)
    end
    if(sx_param=="abs")
        fxch = (xpxch::Float64) -> fxch * max_flux
    end

    ft = (xnet::Float64,xxch::Float64) -> (  (xnet>0) ? xnet+xxch : xxch , (xnet>0) ? xxch : -xnet + xxch )
    ftp = (xnet::Float64,pxxch::Float64) -> ft( xnet , fxch(pxxch) )

    wtx =WX_Transform(w_pairs,ftp)

    # create parameter space polytope:
    # add [0,1] bounds for all wx-parametrized sx fluxes
    # add [-max , max] for wx-parametrized net flux
    # and do the default thing: bound all (non-wx-parametrized) non-sx fluxes by [0,max_flux]
    # and all sx fluxes by [-max_flux,max_flux]
    # !! AND !!
    # Remove the sx parametrized flux from the stoichiometric matrix..
    p_v = get_flux_polyhedron_vanilla(net::Network , lcs::Vector{LinConstr} )

    # remove fluxes w_pairs[:,2] from the stoichiometric matrix:
    p_v.A[:,w_pairs[:,2]] = 0 #zeros(size(wx,1))

    # identify non-system-exchange fluxes
    fnsx = get_idx_non_sx_fluxes(net)
    #print(fnsx)

    lb =  Inf * -ones(length(net.fluxes))
    ub =  Inf *  ones(length(net.fluxes))
    lb[fnsx] = 0
    # if max_val, create bouding box. [0,max_val] for nsx fluxes, AND [-max_val,max_val] for sx fluxes
    if(max_flux<Inf)
        ub[ub.>max_flux] = max_flux
        lb[ get_idx_sx_fluxes(net) ] = -max_flux
    end
    lb[w_pairs[:,2]] = 0
    ub[w_pairs[:,2]] = 1
    lb[w_pairs[:,1]] = -max_flux
    ub[w_pairs[:,1]] = max_flux

    pp = Polyhedron( p_v , bounding_box(lb,ub) )

    return (pp,wtx)
end


text:
using IntervalArithmetic
using Test

@testset "Tests with rational intervals" begin

    a = Interval(1//2, 3//4)
    b = Interval(3//7, 9//12)

    @test a + b == Interval(13//14, 3//2)  # exact

    @test sqrt(a + b) == Interval(0.9636241116594314, 1.2247448713915892)

    X = Interval(1//3)
    @test sqrt(X) == Interval(0.5773502691896257, 0.5773502691896258)

end

setprecision(Interval, 64)

@testset "Rounding rational intervals" begin

    X = Interval(big(1)//3)
    Y = Interval(big"5.77350269189625764452e-01", big"5.77350269189625764561e-01")
    @test sqrt(X) == Y
end

@testset "Tests with float intervals" begin

    c = @floatinterval(0.1, 0.2)

    @test isa(@floatinterval(0.1), Interval) == true
    @test c == Interval(0.09999999999999999, 0.2)

    @test widen(c) == Interval(0.09999999999999998, 0.20000000000000004)

    @test @floatinterval(pi) == Interval(3.141592653589793, 3.1415926535897936)
end

@testset "Testing functions of intervals" begin
    f(x) = x + 0.1

    c = @floatinterval(0.1, 0.2)
    @test f(c) == Interval(0.19999999999999998, 0.30000000000000004)

    d = @interval(0.1, 0.2)
    @test f(d) == @biginterval(0.2, 0.3)
end

@testset "Testing conversions" begin
    f = @interval(0.1, 0.2)
    @test @floatinterval(f) == Interval(0.09999999999999999, 0.2)

    g = @floatinterval(0.1, 0.2)
    @test @interval(g) == Interval(9.9999999999999992e-02, 2.0000000000000001e-01)
end


text:
module Gen2DAgentMotion

include("scene.jl")
include("planner.jl")
include("motion.jl")
include("distributions.jl")

end # module


text:
function step!(unext, uprev, parameters)
    wL, wC, wR = parameters
    numpts = length(unext)
    for i = 2:numpts-1
        unext[i] = wL * uprev[i-1] + wC * uprev[i] + wR * uprev[i]
    end
end

function single_step_error(parameters, u0, uexact)
    T = eltype(parameters)
    unext = zeros(T,length(u0))
    step!(unext, u0, parameters)
    err = sum((unext - uexact) .^ 2)
    return err
end

numpts = 10
u0 = ones(numpts)
uexact = range(0.0,stop=1.0,length=numpts)
parameters = [0.,0.5,0.8]
err = single_step_error(parameters,u0,uexact)

using ForwardDiff
df = ForwardDiff.gradient(x->single_step_error(x,u0,uexact),parameters)

max_stars_repo_name: ArjunNarayanan/ExploreJacobi

max_issues_repo_name: ArjunNarayanan/ExploreJacobi

max_forks_repo_name: ArjunNarayanan/ExploreJacobi


text:
using TermWin
using Dates

TermWin.initsession()
#s = "a test"
v = newTwCalendar( rootTwScreen, today(), posy=:random, posx=:random )
v.title = "Input: "
activateTwObj( rootTwScreen )
ret = v.value
TermWin.endsession()
println( "You entered ", string( ret ) )

v = nothing
TermWin.initsession()
v = newTwCalendar( rootTwScreen, today(), posy=:random, posx=:random, ncalStyle=false )
activateTwObj( rootTwScreen )
ret = v.value
TermWin.endsession()
println( "You entered ", string( ret ) )