text:
using Ferrite, SparseArrays, LinearAlgebra, FerriteViz

struct J2Plasticity{T, S <: SymmetricTensor{4, 3, T}}
    G::T  # Shear modulus
    K::T  # Bulk modulus
    σ₀::T # Initial yield limit
    H::T  # Hardening modulus
    Dᵉ::S # Elastic stiffness tensor
end;

function J2Plasticity(E, ν, σ₀, H)
    δ(i,j) = i == j ? 1.0 : 0.0 # helper function
    G = E / 2(1 + ν)
    K = E / 3(1 - 2ν)

    Isymdev(i,j,k,l) = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
    temp(i,j,k,l) = 2.0G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l))
    Dᵉ = SymmetricTensor{4, 3}(temp)
    return J2Plasticity(G, K, σ₀, H, Dᵉ)
end;

mutable struct MaterialState{T, S <: SecondOrderTensor{3, T}}
    # Store "converged" values
    ϵᵖ::S # plastic strain
    σ::S # stress
    k::T # hardening variable

    # Store temporary values used during equilibrium iterations
    temp_ϵᵖ::S
    temp_σ::S
    temp_k::T
end

function MaterialState()
    return MaterialState(
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                0.0,
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                0.0)
end

function update_state!(state::MaterialState)
    state.ϵᵖ = state.temp_ϵᵖ
    state.σ = state.temp_σ
    state.k = state.temp_k
end;

function vonMises(σ)
    s = dev(σ)
    return sqrt(3.0/2.0 * s ⊡ s)
end;

function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticity, state::MaterialState)
    # unpack some material parameters
    G = material.G
    K = material.K
    H = material.H

    # We use (•)ᵗ to denote *trial*-values
    σᵗ = material.Dᵉ ⊡ (ϵ - state.ϵᵖ) # trial-stress
    sᵗ = dev(σᵗ)         # deviatoric part of trial-stress
    J₂ = 0.5 * sᵗ ⊡ sᵗ   # second invariant of sᵗ
    σᵗₑ = sqrt(3.0*J₂)   # effetive trial-stress (von Mises stress)
    σʸ = material.σ₀ + H * state.k # Previous yield limit

    φᵗ  = σᵗₑ - σʸ # Trial-value of the yield surface

    if φᵗ < 0.0 # elastic loading
        state.temp_σ = σᵗ
        return state.temp_σ, material.Dᵉ
    else # plastic loading
        h = H + 3G
        μ =  φᵗ / h   # plastic multiplier

        c1 = 1 - 3G * μ / σᵗₑ
        s = c1 * sᵗ           # updated deviatoric stress
        σ = s + vol(σᵗ)       # updated stress

        # Compute algorithmic tangent stiffness ``D = \frac{\Delta \sigma }{\Delta \epsilon}``
        κ = H * (state.k + μ) # drag stress
        σₑ = material.σ₀ + κ  # updated yield surface

        δ(i,j) = i == j ? 1.0 : 0.0
        Isymdev(i,j,k,l)  = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
        Q(i,j,k,l) = Isymdev(i,j,k,l) - 3.0 / (2.0*σₑ^2) * s[i,j]*s[k,l]
        b = (3G*μ/σₑ) / (1.0 + 3G*μ/σₑ)

        Dtemp(i,j,k,l) = -2G*b * Q(i,j,k,l) - 9G^2 / (h*σₑ^2) * s[i,j]*s[k,l]
        D = material.Dᵉ + SymmetricTensor{4, 3}(Dtemp)

        # Store outputs in the material state
        Δϵᵖ = 3/2 *μ / σₑ*s            # plastic strain
        state.temp_ϵᵖ = state.ϵᵖ + Δϵᵖ  # plastic strain
        state.temp_k = state.k + μ     # hardening variable
        state.temp_σ = σ               # updated stress
        return state.temp_σ, D
    end
end

function create_values(interpolation)
    # setup quadrature rules
    qr      = QuadratureRule{3,RefTetrahedron}(2)
    face_qr = QuadratureRule{2,RefTetrahedron}(3)

    # create geometric interpolation (use the same as for u)
    interpolation_geom = Lagrange{3,RefTetrahedron,1}()

    # cell and facevalues for u
    cellvalues_u = CellVectorValues(qr, interpolation, interpolation_geom)
    facevalues_u = FaceVectorValues(face_qr, interpolation, interpolation_geom)

    return cellvalues_u, facevalues_u
end;

function create_dofhandler(grid, interpolation)
    dh = DofHandler(grid)
    dim = 3
    push!(dh, :u, dim, interpolation) # add a displacement field with 3 components
    close!(dh)
    return dh
end

function create_bc(dh, grid)
    dbcs = ConstraintHandler(dh)
    # Clamped on the left side
    dofs = [1, 2, 3]
    dbc = Dirichlet(:u, getfaceset(grid, "left"), (x,t) -> [0.0, 0.0, 0.0], dofs)
    add!(dbcs, dbc)
    close!(dbcs)
    return dbcs
end;

function doassemble(cellvalues::CellVectorValues{dim},
                    facevalues::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
                    dh::DofHandler, material::J2Plasticity, u, states, t) where {dim}
    r = zeros(ndofs(dh))
    assembler = start_assemble(K, r)
    nu = getnbasefunctions(cellvalues)
    re = zeros(nu)     # element residual vector
    ke = zeros(nu, nu) # element tangent matrix

    for (cell, state) in zip(CellIterator(dh), states)
        fill!(ke, 0)
        fill!(re, 0)
        eldofs = celldofs(cell)
        ue = u[eldofs]
        assemble_cell!(ke, re, cell, cellvalues, facevalues, grid, material,
                       ue, state, t)
        assemble!(assembler, eldofs, re, ke)
    end
    return K, r
end

function assemble_cell!(Ke, re, cell, cellvalues, facevalues, grid, material,
                        ue, state, t)
    n_basefuncs = getnbasefunctions(cellvalues)
    reinit!(cellvalues, cell)

    for q_point in 1:getnquadpoints(cellvalues)
        # For each integration point, compute stress and material stiffness
        ∇u = function_gradient(cellvalues, q_point, ue)
        ϵ = symmetric(∇u) # Total strain
        σ, D = compute_stress_tangent(ϵ, material, state[q_point])

        dΩ = getdetJdV(cellvalues, q_point)
        for i in 1:n_basefuncs
            δϵ = symmetric(shape_gradient(cellvalues, q_point, i))

            re[i] += (δϵ ⊡ σ) * dΩ # add internal force to residual
            for j in 1:i
                Δϵ = symmetric(shape_gradient(cellvalues, q_point, j))
                Ke[i, j] += δϵ ⊡ D ⊡ Δϵ * dΩ
            end
        end
    end
    symmetrize_lower!(Ke)

    # Add traction as a negative contribution to the element residual `re`:
    for face in 1:nfaces(cell)
        if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "right")
            reinit!(facevalues, cell, face)
            for q_point in 1:getnquadpoints(facevalues)
                dΓ = getdetJdV(facevalues, q_point)
                for i in 1:n_basefuncs
                    δu = shape_value(facevalues, q_point, i)
                    re[i] -= (δu ⋅ t) * dΓ
                end
            end
        end
    end

end

function symmetrize_lower!(K)
    for i in 1:size(K,1)
        for j in i+1:size(K,1)
            K[i,j] = K[j,i]
        end
    end
end;

function solve(liveplotting=false)
    # Define material parameters
    E = 200.0e9 # [Pa]
    H = E/20   # [Pa]
    ν = 0.3     # [-]
    σ₀ = 200e6  # [Pa]
    material = J2Plasticity(E, ν, σ₀, H)

    L = 10.0 # beam length [m]
    w = 1.0  # beam width [m]
    h = 1.0  # beam height[m]
    n_timesteps = 100
    u_max = zeros(n_timesteps)
    traction_magnitude = 1.e7 * range(0.5, 1.0, length=n_timesteps)

    # Create geometry, dofs and boundary conditions
    n = 2
    nels = (10n, n, 2n) # number of elements in each spatial direction
    P1 = Vec((0.0, 0.0, 0.0))  # start point for geometry
    P2 = Vec((L, w, h))        # end point for geometry
    grid = generate_grid(Tetrahedron, nels, P1, P2)
    interpolation = Lagrange{3, RefTetrahedron, 1}() # Linear tet with 3 unknowns/node

    dh = create_dofhandler(grid, interpolation) # JuaFEM helper function
    dbcs = create_bc(dh, grid) # create Dirichlet boundary-conditions

    cellvalues, facevalues = create_values(interpolation)

    # Pre-allocate solution vectors, etc.
    n_dofs = ndofs(dh)  # total number of dofs
    u  = zeros(n_dofs)  # solution vector
    u_history = Vector{Vector{Float64}}()
    if liveplotting
        plotter = MakiePlotter(dh,u)
        fig = ferriteviewer(plotter)
        display(fig)
    end
    Δu = zeros(n_dofs)  # displacement correction
    r = zeros(n_dofs)   # residual
    K = create_sparsity_pattern(dh); # tangent stiffness matrix

    # Create material states. One array for each cell, where each element is an array of material-
    # states - one for each integration point
    nqp = getnquadpoints(cellvalues)
    states = [[MaterialState() for _ in 1:nqp] for _ in 1:getncells(grid)]

    # states = [MaterialState() for _ in 1:nqp for _ in 1:getncells(grid)]
    # temp_states = [MaterialState() for _ in 1:nqp for _ in 1:getncells(grid)]

    # Newton-Raphson loop
    NEWTON_TOL = 1 # 1 N

    for timestep in 1:n_timesteps
        t = timestep # actual time (used for evaluating d-bndc)
        traction = Vec((0.0, 0.0, traction_magnitude[timestep]))
        newton_itr = -1
        update!(dbcs, t) # evaluates the D-bndc at time t
        apply!(u, dbcs)  # set the prescribed values in the solution vector

        while true; newton_itr += 1

            if newton_itr > 8
                error("Reached maximum Newton iterations, aborting")
                break
            end
            K, r = doassemble(cellvalues, facevalues, K, grid, dh, material, u,
                             states, traction);
            norm_r = norm(r[Ferrite.free_dofs(dbcs)])

            if norm_r < NEWTON_TOL
                break
            end

            apply_zero!(K, r, dbcs)
            Δu = Symmetric(K) \ r
            u -= Δu
        end
        if liveplotting
            FerriteViz.update!(plotter,u)
            sleep(0.1)
        end
        push!(u_history,u)

        # Update all the material states after we have reached equilibrium
        for cell_states in states
            foreach(update_state!, cell_states)
        end
        u_max[timestep] = max(abs.(u)...) # maximum displacement in current timestep
    end

    # ## Postprocessing
    # Only a vtu-file corrsponding to the last time-step is exported.
    #
    # The following is a quick (and dirty) way of extracting average cell data for export.
    mises_values = zeros(getncells(grid))
    κ_values = zeros(getncells(grid))
    for (el, cell_states) in enumerate(states)
        for state in cell_states
            mises_values[el] += vonMises(state.σ)
            κ_values[el] += state.k*material.H
        end
        mises_values[el] /= length(cell_states) # average von Mises stress
        κ_values[el] /= length(cell_states)     # average drag stress
    end
    return u, dh, u_history, mises_values, κ_values
end

ext: jl


text:
using Neurodevo
using CGP

CGP.Config.init("cfg/cgp/base.yaml")
CGP.Config.init("cfg/cgp/functions.yaml")

function make_chromos!(cfg::Dict; cinds::Array{Int64}=collect(1:10))
    chromos = Array{CGP.PCGPChromo}(undef, length(cinds))
    nins = Neurodevo.input_lengths(cfg)
    nouts = Neurodevo.output_lengths(cfg)
    for i in eachindex(chromos)
        chromos[i] = CGP.PCGPChromo(nins[cinds[i]], nouts[cinds[i]])
    end
    cfg["chromosomes"] = chromos
    chromos
end

function make_chromos!(cfg::Dict, genes::Array{Array{Float64}};
                       cinds::Array{Int64}=collect(1:10))
    chromos = Array{CGP.PCGPChromo}(undef, length(cinds))
    nins = Neurodevo.input_lengths(cfg)
    nouts = Neurodevo.output_lengths(cfg)
    for i in eachindex(cinds)
        chromos[i] = CGP.PCGPChromo(genes[i], nins[cinds[i]], nouts[cinds[i]])
    end
    cfg["chromosomes"] = chromos
    chromos
end

function make_controller(cfg::Dict)
    chromos = cfg["chromosomes"]
    nouts = Neurodevo.output_lengths(cfg)
    cell_division(x::Array{Float64}) = false
    new_cell_params(x::Array{Float64}) = zeros(nouts[2])
    cell_death(x::Array{Float64}) = false
    cell_state_update(x::Array{Float64}) = process(chromos[1], x)
    cell_param_update(x::Array{Float64}) = process(chromos[2], x)
    connect(x::Array{Float64}) = false
    new_conn_params(x::Array{Float64}) = zeros(nouts[7])
    disconnect(x::Array{Float64}) = false
    conn_state_update(x::Array{Float64}) = process(chromos[3], x)
    conn_param_update(x::Array{Float64}) = process(chromos[4], x)

    Controller(cell_division, new_cell_params, cell_death,
               cell_state_update, cell_param_update,
               connect, new_conn_params, disconnect,
               conn_state_update, conn_param_update)
end

function cgp_cfg(cfg::Dict, genes::Array{Float64})
    cfg["T_learn"] = max(1, round(Int64, genes[1] * 5) * 5)
    cfg["T_devo"] = max(1, round(Int64, genes[2] * 5) * 5)
    cfg["n_hidden"] = round(Int64, genes[3] * 10) * 2
    cfg["n_step"] = max(1, round(Int64, genes[4] * 5))
    cfg
end

function cgp_controller(cfg::Dict; cinds::Array{Int64}=collect(1:10))
    make_chromos!(cfg; cinds=cinds)
    make_controller(cfg)
end

function cgp_controller(cfg::Dict, genes::Array{Array{Float64}};
                        cinds::Array{Int64}=collect(1:10))
    make_chromos!(cfg, genes; cinds=cinds)
    make_controller(cfg)
end


text:
"""
a Double if statement
"""

function midCriterion(a, b) #ok
    m = sumInterval(a, b) # | b - a |   definition
    condition = isEven(m) #even (divisible by 2)
    check = nothing
    if condition == true if iseven #ERROR double if!  #<----------------
        check = condition
    elseif condition == false
        check = condition
    else #if faulty input or Unexpected ERROR Occured
        check
    end
    return check #whether check is true, false, nothing
end

function midCriterion(a, b) #ok
    m = sumInterval(a, b) # | b - a |   definition
    condition = isEven2(m) #even (divisible by 2)
    check = nothing
    if condition == true #ERROR double if!  #<----------------
        check = condition
    elseif condition == false
        check = condition
    else #if faulty input or Unexpected ERROR Occured
        check
    end
    return check #whether check is true, false, nothing
end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
sqlitetype{T<:Integer}(::Type{T}) = "INT"
sqlitetype{T<:AbstractFloat}(::Type{T}) = "REAL"
sqlitetype{T<:AbstractString}(::Type{T}) = "TEXT"
sqlitetype(::Type{NullType}) = "NULL"
sqlitetype(x) = "BLOB"

"""
independent SQLite.Sink constructor to create a new or wrap an existing SQLite table with name `name`.
must provide a `Data.Schema` through the `schema` argument
can optionally provide an existing SQLite table name or new name that a created SQLite table will be called through the `name` argument
`temp=true` will create a temporary SQLite table that will be destroyed automatically when the database is closed
`ifnotexists=false` will throw an error if `name` already exists in `db`
"""
function Sink(db::DB, schema::Data.Schema; name::AbstractString="julia_"*randstring(), temp::Bool=false, ifnotexists::Bool=true, append::Bool=false)
    rows, cols = size(schema)
    temp = temp ? "TEMP" : ""
    ifnotexists = ifnotexists ? "IF NOT EXISTS" : ""
    columns = [string(esc_id(schema.header[i]), ' ', sqlitetype(schema.types[i])) for i = 1:cols]
    SQLite.execute!(db, "CREATE $temp TABLE $ifnotexists $(esc_id(name)) ($(join(columns, ',')))")
    !append && execute!(db, "delete from $(esc_id(name))")
    params = chop(repeat("?,", cols))
    stmt = SQLite.Stmt(db, "INSERT INTO $(esc_id(name)) VALUES ($params)")
    return Sink(schema, db, name, stmt, "")
end

"constructs a new SQLite.Sink from the given `SQLite.Source`; uses `source` schema to create the SQLite table"
function Sink(source::SQLite.Source; name::AbstractString="julia_"*randstring(), temp::Bool=false, ifnotexists::Bool=true)
    return Sink(source.db, source.schema; name=name, temp=temp, ifnotexists=ifnotexists)
end
"constructs a new SQLite.Sink from the given `Data.Source`; uses `source` schema to create the SQLite table"
function Sink(db::DB, source; name::AbstractString="julia_"*randstring(), temp::Bool=false, ifnotexists::Bool=true)
    return Sink(db, Data.schema(source); name=name, temp=temp, ifnotexists=ifnotexists)
end

# DataStreams interface
Data.streamtypes{T<:SQLite.Sink}(::Type{T}) = [Data.Field]

function Sink{T}(source, ::Type{T}, append::Bool, db::DB, name::AbstractString="julia_" * randstring())
    sink = Sink(db, Data.schema(source); name=name, append=append)
    return sink
end
function Sink{T}(sink, source, ::Type{T}, append::Bool)
    !append && execute!(sink.db, "delete from $(esc_id(sink.tablename))")
    return sink
end

function getbind!{T}(val::Nullable{T}, col, stmt)
    if isnull(val)
        SQLite.bind!(stmt, col, NULL)
    else
        SQLite.bind!(stmt, col, get(val))
    end
    return
end
if isdefined(:isna)
    function getbind!{T}(val::T, col, stmt)
        if isna(val)
            SQLite.bind!(stmt, col, NULL)
        else
            SQLite.bind!(stmt, col, val)
        end
        return
    end
else
    getbind!{T}(val::T, col, stmt) = SQLite.bind!(stmt, col, val)
end

function Data.open!(sink::SQLite.Sink, source)
    execute!(sink.db,"PRAGMA synchronous = OFF;")
    sink.transaction = string("SQLITE",randstring(10))
    transaction(sink.db, sink.transaction)
    return nothing
end

function Data.streamfield!{T}(sink::SQLite.Sink, source, ::Type{T}, row, col, cols)
    val = Data.getfield(source, T, row ,col)
    getbind!(val, col, sink.stmt)
    if col == cols
        SQLite.sqlite3_step(sink.stmt.handle)
        SQLite.sqlite3_reset(sink.stmt.handle)
    end
    return nothing
end

function Data.cleanup!(sink::SQLite.Sink)
    rollback(sink.db, sink.transaction)
    execute!(sink.db, "PRAGMA synchronous = ON;")
    return nothing
end

function Data.close!(sink::SQLite.Sink)
    commit(sink.db, sink.transaction)
    execute!(sink.db, "PRAGMA synchronous = ON;")
    SQLite.execute!(sink.db, "ANALYZE $(esc_id(sink.tablename))")
    return nothing
end

"""
Load a Data.Source `source` into an SQLite table that will be named `tablename` (will be auto-generated if not specified).

`temp=true` will create a temporary SQLite table that will be destroyed automatically when the database is closed
`ifnotexists=false` will throw an error if `tablename` already exists in `db`
"""
function load{T}(db, name, ::Type{T}, args...;
              temp::Bool=false,
              ifnotexists::Bool=true,
              append::Bool=false)
    source = T(args...)
    schema = Data.schema(source)
    sink = Sink(db, schema; name=name, temp=temp, ifnotexists=ifnotexists, append=append)
    Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end
function load{T}(db, name, source::T;
              temp::Bool=false,
              ifnotexists::Bool=true,
              append::Bool=false)
    sink = Sink(db, Data.schema(source); name=name, temp=temp, ifnotexists=ifnotexists, append=append)
    Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end

load{T}(sink::Sink, ::Type{T}, args...; append::Bool=false) = (sink = Data.stream!(T(args...), sink, append); Data.close!(sink); return sink)
load(sink::Sink, source; append::Bool=false) = (sink = Data.stream!(source, sink, append); Data.close!(sink); return sink)


text:
include("newton.jl")

include("cubic_poly.jl")

@time matrix = paranewton(roots3,colors3,cubic,20,0.01,(-1.0,1.0,-1.0,1.0),500,500);
@time matrix = paranewton(roots3,colors3,cubic,20,0.01,(-1.0,1.0,-1.0,1.0),500,500);

imwrite(convert(Image,getimage(matrix)), "img_cubic.png")

include("sin.jl")

@time matrix = paranewton(roots_sin,colors_sin,fsin,20,0.01,(-3.0,3.0,-3.0,3.0),500,500);
@time matrix = paranewton(roots_sin,colors_sin,fsin,20,0.01,(-3.0,3.0,-3.0,3.0),500,500);

imwrite(convert(Image,getimage(matrix)), "img_sin.png")

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 2d89ea8d-697c-4278-b467-d5ed7e123dd7
using Plots, Latexify, LaTeXStrings, PlutoUI

# ╔═╡ 751e5150-a132-11eb-2a66-5d81e11face8
md""" # Mínimos quadrados e ortogonalidade

Vamos ilustrar o processo algébrico da aproximação por mínimos quadrados """

# ╔═╡ 6c41c486-376e-42b1-9d0f-89178e1020a2
md" Primeiro vamos construir os vetores de dados, vamos dividir o intervalo $[0,10]$ em $n$ partes "

# ╔═╡ 3fc129c6-b277-4af5-b114-1b33a9f6ebcf
md" n = $(@bind n Slider(1:100, default=10, show_value=true))"

# ╔═╡ ec6d7196-ef3a-4c6b-a673-127c22a44aeb
x=range(0,10,step=10/n)

# ╔═╡ 0e86ddce-1b8c-4158-9713-9bc209089091
y=(10*(0.5 .- rand(n+1))) + 2*x + ones(n+1) + x.^2

# ╔═╡ d1982bed-31dc-4cbc-a11b-8bafa9295c9c

# ╔═╡ 34ed10fa-7241-475e-9611-402931dd3784
scatter(x,y, label=false)

# ╔═╡ b2e7caf3-3e6d-4fe7-a063-5820e86de6e6
md""" Agora quero saber qual é a melhor reta $y=ax+b$ que melhor se ajusta à essa massa de dados que minimiza $\sum_{i=0}^n (ax_i-b -y_i)^2= \| A\begin{bmatrix} a \\ b \end{bmatrix} - y\|^2$.

Pela teoria que acabamos de ver:
$\begin{bmatrix} a \\ b \end{bmatrix} = (A^tA)^{-1}A^ty$

"""

# ╔═╡ ce16bd9e-0753-43bd-af54-c2494b09eb46

# ╔═╡ 50886dfa-68a6-4c2a-8583-f5e6240c7364
begin
	um=ones(n+1);
	A=[x um];
	latexstring(latexify(A))
end

# ╔═╡ c5570f33-d189-43e2-b1af-2e892067e9d3
md"""
```math
A^tA
=
``` $(latexstring(latexify(A'*A)))
"""

# ╔═╡ 6566ec1b-b9b8-4210-8e43-efa816db54cc
md"""
```math
A^ty
=
``` $(latexstring(latexify(A'*y)))
"""

# ╔═╡ 42230d2b-cc08-414d-bf0f-610344cd4890
md""" ```math
\begin{bmatrix} a \\ b \end{bmatrix} =
``` $(latexstring(latexify(inv(A'*A)*(A'*y))))
"""

# ╔═╡ bd08d6d3-726d-47c5-aab4-d4a86606e99d
begin
	m=inv(A'*A)*(A'*y);
	r= m[1]*x+ m[2]*um;
end

# ╔═╡ 9b683dd1-cc8f-4129-ad44-2f94817e632e
begin
	plot(x,r,label=false)
	scatter!(x,y,label=false)
end

# ╔═╡ 82eb86ba-47aa-43bd-8b9f-7a42e0f944bf
md"""

 **y** = $(m[1])**x** + $(m[2])

é a equação da reta que melhor se ajusta aos dados
"""

# ╔═╡ 2a1ef996-4516-4fbc-8b6a-8b5fdc12f0ba
md" Para o ajuste de uma parábola: "

# ╔═╡ 30a7bba5-9105-4564-8f61-0ff49c21801e
x2=x.^2

# ╔═╡ 48382d7f-a66e-43a2-8cee-3ba1f79d2690
begin
	A2 =[x2 x um]
	N=A2'*A2
	p=inv(N)*A2'*y
	q=p[1]*x2 + p[2]*x+p[3]*um
end

# ╔═╡ 41271278-dc7f-42ee-a565-eb48afdf3dda
begin
	plot(x,q,label=false)
	plot!(x,r,label=false)
	scatter!(x,y,label=false)
end

# ╔═╡ Cell order:
# ╟─751e5150-a132-11eb-2a66-5d81e11face8
# ╟─2d89ea8d-697c-4278-b467-d5ed7e123dd7
# ╟─6c41c486-376e-42b1-9d0f-89178e1020a2
# ╟─3fc129c6-b277-4af5-b114-1b33a9f6ebcf
# ╟─ec6d7196-ef3a-4c6b-a673-127c22a44aeb
# ╟─0e86ddce-1b8c-4158-9713-9bc209089091
# ╟─d1982bed-31dc-4cbc-a11b-8bafa9295c9c
# ╟─34ed10fa-7241-475e-9611-402931dd3784
# ╟─b2e7caf3-3e6d-4fe7-a063-5820e86de6e6
# ╠═ce16bd9e-0753-43bd-af54-c2494b09eb46
# ╟─50886dfa-68a6-4c2a-8583-f5e6240c7364
# ╟─c5570f33-d189-43e2-b1af-2e892067e9d3
# ╟─6566ec1b-b9b8-4210-8e43-efa816db54cc
# ╟─42230d2b-cc08-414d-bf0f-610344cd4890
# ╟─bd08d6d3-726d-47c5-aab4-d4a86606e99d
# ╟─9b683dd1-cc8f-4129-ad44-2f94817e632e
# ╟─82eb86ba-47aa-43bd-8b9f-7a42e0f944bf
# ╟─2a1ef996-4516-4fbc-8b6a-8b5fdc12f0ba
# ╟─30a7bba5-9105-4564-8f61-0ff49c21801e
# ╟─48382d7f-a66e-43a2-8cee-3ba1f79d2690
# ╟─41271278-dc7f-42ee-a565-eb48afdf3dda


text:
if :recursion_relation in fieldnames(Method)

first(methods(Diffractor.∂⃖recurse{1}())).recursion_relation = function(method1, method2, parent_sig, new_sig)
	# Recursion from a higher to a lower order is always allowed
	parent_order = parent_sig.parameters[1].parameters[1]
	child_order = new_sig.parameters[1].parameters[1]
	#@Core.Main.Base.show (parent_order, child_order)
	if parent_order > child_order
		return true
	end
	wrapped_parent_sig = Tuple{parent_sig.parameters[2:end]...}
	wrapped_new_sig = Tuple{parent_sig.parameters[2:end]...}
	if method2 !== nothing && isdefined(method2, :recursion_relation)
		# TODO: What if method2 is itself a generated function.
		return method2.recursion_relation(method2, nothing, wrapped_parent_sig, wrapped_new_sig)
	end
	return Core.Compiler.type_more_complex(new_sig, parent_sig, Core.svec(parent_sig), 1, 3, length(method1.sig.parameters)+1)
end

first(methods(PrimeDerivativeBack(sin))).recursion_relation = function(method1, method2, parent_sig, new_sig)
	# Recursion from a higher to a lower order is always allowed
	parent_order = parent_sig.parameters[1].parameters[1]
	child_order = new_sig.parameters[1].parameters[1]
	#@Core.Main.Base.show (parent_order, child_order)
	if parent_order > child_order
		return true
	end
	wrapped_parent_sig = Tuple{parent_sig.parameters[2:end]...}
	wrapped_new_sig = Tuple{parent_sig.parameters[2:end]...}
	return Core.Compiler.type_more_complex(new_sig, parent_sig, Core.svec(parent_sig), 1, 3, length(method1.sig.parameters)+1)
end

which(Tuple{∂⃖{N}, T, Vararg{Any}} where {T,N}).recursion_relation = function(_, _, parent_sig, new_sig)
	# Any actual recursion will always be caught be one of the functions we're
	# recursing into.
	return isa(Base.unwrap_unionall(parent_sig.parameters[1].parameters[1]), Int) &&
           isa(Base.unwrap_unionall(new_sig.parameters[1].parameters[1]), Int)
end

which(Tuple{∂⃖{N}, ∂⃖{1}, Vararg{Any}} where {N}).recursion_relation = function(_, _, parent_sig, new_sig)
	# Allowed as long as both parent and new sig have concrete integers. In that
	# case, actual recursion will be caught elsewhere.
	return isa(Base.unwrap_unionall(parent_sig.parameters[1].parameters[1]), Int) &&
	       isa(Base.unwrap_unionall(new_sig.parameters[1].parameters[1]), Int)
end

end


text:
# ---
# author: "[Johnny Chen](https://github.com/johnnychen94); Jane Doe"
# ---

# author with markdown url


text:
using Statistics
print("Hello")
"blabla"


text:
using ManifoldHMC

using SimpleTraits
using Random
using LinearAlgebra
using Manifolds: Manifold, Euclidean, Sphere, IsDecoratorManifold
using Manifolds: Manifold, IsDecoratorManifold
import Manifolds: base_manifold, exp!
using StaticArrays
using ForwardDiff
using Test

include("flow.jl")
include("manifolds.jl")

using Test


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE

using Pkg, Documenter, HeatTransfer, Literate, Dates

# automatically generate documentation from tests

"""
    update_datetime(content)

Replace DATETIMEOFTODAY with current datetime.
"""
function update_datetime(content)
    content = replace(content, "DATETIMEOFTODAY" => DateTime(now()))
    return content
end

"""
    add_datetime(content)

Add page generation time to the end of the content.
"""
function add_datetime(content)
    line =  "\n# Page generated at " * string(DateTime(now())) * "."
    content = content * line
    println(content)
    return content
end

"""
    remove_license(content)

Remove licence strings from source file.
"""
function remove_license(content)
    # return replace(content, r"# This file is a part of JuliaFEM.*?\n\n"sm, "")
    lines = split(content, '\n')
    function islicense(line)
        occursin("# This file is a part of JuliaFEM.", line) && return false
        occursin("# License is MIT:", line) && return false
        return true
    end
    content = join(filter(islicense, lines), '\n')
    return content
end

function generate_docs()

    function preprocess(content)
        content = add_datetime(content)
        content = remove_license(content)
    end

    pkg_dir = dirname(dirname(pathof(HeatTransfer)))
    testdir = joinpath(pkg_dir, "test")
    outdir = joinpath(pkg_dir, "docs", "src", "tests")
    test_pages = []
    for test_file in readdir(testdir)
        startswith(test_file, "test_") || continue
        Literate.markdown(joinpath(testdir, test_file), outdir; documenter=true, preprocess=preprocess)
        generated_test_file = joinpath("tests", first(splitext(test_file)) * ".md")
        push!(test_pages, generated_test_file)
    end
    return test_pages

end

test_pages = generate_docs()

makedocs(modules=[HeatTransfer],
         format = Documenter.HTML(),
         checkdocs = :all,
         sitename = "HeatTransfer.jl",
         pages = [
                  "index.md",
                  "Examples" => test_pages
                 ]
        )


text:
using Genie.Router
using TasksController

route("/") do
  serve_static_file("welcome.html")
end

route("/hello") do
    "
Welcome to Genie!
"
end

route("/tasks", TasksController.index, named = :get_tasks)
route("/tasks/create", TasksController.create, method = POST, named = :create_task)
route("/tasks/edit", TasksController.edit, named = :edit_task)
route("/tasks/update", TasksController.update, method = POST, named = :update_task)
route("/tasks/ajax_create", TasksController.ajax_create, method = POST, named = :ajax_create_task)
route("/tasks/update_status", TasksController.ajax_update_status, method = POST, named = :ajax_update_status_task)
route("/tasks/destroy", TasksController.ajax_destroy, method = POST, named = :ajax_destroy_task)
route("/tasks/change_show_done", TasksController.ajax_change_show_done, method = POST, named = :ajax_change_show_done)

max_stars_repo_name: ky-s/MyGenieApp

max_issues_repo_name: ky-s/MyGenieApp

max_forks_repo_name: ky-s/MyGenieApp


text:
using Base.Test

include("../../build_out/generated/Shape.jl")
using FFIG

@test area(Square(3.0)) == 9.0

@test perimeter(Square(3.0)) == 12.0

@test name(Circle(3.0)) == "Circle"

max_stars_repo_name: FFIG/ffig

max_issues_repo_name: FFIG/ffig

max_forks_repo_name: FFIG/ffig


text:
# PROJECT_GRAPH    Project a point onto the graph of a linear operator.
#
#   If A is an m x n matrix and v = [v1;v2] is a vector such that v1 is of
#   length n and v2 is of length m, then
#
#     z = project_graph(v,A);
#
#   gives the projection of v onto the graph { (x,y) | y = Ax } of A. The
#   vector z can then be sliced into its first n components and the remaining
#   m components.
#
#   If A is dense and m <= n, then the projection involves computing
#   the Cholesky factor L of I + A*A'. One can obtain A*A' and L by calling
#
#     [z AA L] = project_graph(v,A)
#
#   and can pass in AA and L to avoid them being recomputed via
#
#     z = project_graph(v,A,AA,L)
#
#   This allows for factorization caching when calling project_graph
#   multiple times for the same A.
#
#   This applies similarly when m >= n, only AA is A'*A and L is the
#   Cholesky factor of I + A'*A.
#
#   This also applies similarly when A is sparse, in which case we
#   use a permuted LDL factorization of some matrix. In this case,
#   one can obtain P, L, and D via
#
#     [z L D P] = project_graph(v,A);
#
#   and can pass these back in via
#
#     z = project_graph(v,A,[],L,D,P);
#
#   Note that L has different meanings depending on whether or not A is sparse.
function project_graph(
    v,
    A,
    AA,
    L,
    D,
    P
)
    m, n = size(A)
    c = v[1:n]
    d = v[(n + 1):end]

    if issparse(A)
        if isempty(P) || isempty(L) || isempty(D)
            K = [speye(n) A';
                 A - speye(m)]
            L, D, P = ldl(K)
        end

        z = P * (L' \ (D \ (L \ (P' * sparse([c + A' * d; zeros(m, 1)])))))

        varargout[1] = {z}
        varargout[2] = {L}
        varargout[3] = {D}
        varargout[4] = {P}
    else
        if m <= n
            if isempty(AA)
                AA = A * A'
            end
            if isempty(L)
                L = chol(eye(m) + AA)
            end
            y = L \ (L' \ (A * c + AA * d))
            x = c + A' * (d - y)
        else
            if isempty(AA)
                AA = A' * A
            end
            if isempty(L)
                L = chol(eye(n) + AA)
            end
            x = L \ (L' \ (c + A' * d))
            y = A * x
        end

        varargout[1] = {[x;y]}
        varargout[2] = {AA}
        varargout[3] = {L}
    end

    return varargout
end


text:
@testset "AdaptiveParticleSwarm Tuning Hyperparameters" begin
    warning = "AdaptiveParticleSwarm requires at least 3 particles. " *
              "Resetting n_particles=3. " *
              "AdaptiveParticleSwarm requires 1.5 ≤ c1 ≤ 2.5, 1.5 ≤ c2 ≤ 2.5, and " *
              "c1 + c2 ≤ 4. Resetting coefficients c1=2.0, c2=2.0. " *
              "AdaptiveParticleSwarm requires 0 ≤ prob_shift < 1. " *
              "Resetting prob_shift=0.25. "
    ps = @test_logs (:warn, warning) AdaptiveParticleSwarm(
        n_particles=2, c1=3, c2=3, prob_shift=2
    )
    @test ps.n_particles == 3
    @test ps.c1 == 2.0
    @test ps.c2 == 2.0
    @test ps.prob_shift == 0.25
end

@testset "Evolutionary Algorithm" begin
    # Evolutionary factor
    X = [1 1 1; 2 2 2; 10 10 10]
    # If the global best particle is the closest to all particles, then f is 0
    @test PSO._evolutionary_factor(X, 2) == 0
    # If the global best particle is the furthest from all particles, then f is 1
    @test PSO._evolutionary_factor(X, 3) == 1

    # Evolutionary phase
    # Initially pick most likely phase
    @test PSO._evolutionary_phase(0.55, nothing) == 1 # marginal probs [0.75, 0.25, 0, 0]
    @test PSO._evolutionary_phase(0.25, nothing) == 2 # marginal probs [0, 0.5, 0.25, 0]
    @test PSO._evolutionary_phase(0.225, nothing) == 3 # marginal probs [0, 0.25, 0.375, 0]
    @test PSO._evolutionary_phase(0.775, nothing) == 4 # marginal probs [0.25, 0, 0, 0.375]
    # Move to the next phase if possible
    @test PSO._evolutionary_phase(0.55, 1) == 2
    @test PSO._evolutionary_phase(0.25, 2) == 3
    @test PSO._evolutionary_phase(0.775, 4) == 1
    # Stay in the current phase if possible and moving to the next is impossible
    @test PSO._evolutionary_phase(0.55, 2) == 2
    @test PSO._evolutionary_phase(0.25, 3) == 3
    @test PSO._evolutionary_phase(0.775, 1) == 1
    # Pick the most likely phase otherwise
    @test PSO._evolutionary_phase(0.55, 3) == 1
    @test PSO._evolutionary_phase(0.25, 4) == 2
    @test PSO._evolutionary_phase(0.225, 4) == 3
    @test PSO._evolutionary_phase(0.775, 2) == 4

    # Cognitive and social coefficients clamping
    @test PSO._clamp_coefficients(1.0, 1.0) == (1.5, 1.5) # lower bound is 1.5
    @test PSO._clamp_coefficients(3.0, 1.0) == (2.5, 1.5) # upper bound is 2.5
    @test PSO._clamp_coefficients(2.25, 2.25) == (2.0, 2.0) # sum cannot be larger than 4

    # Coefficient adaptive control
    rng = StableRNG(1234)
    # Exploration state
    f, phase = 0.65, 1
    w, c1, c2 = PSO._adapt_parameters(rng, 1.75, 1.75, f, phase)
    @test w == 1 / (1 + 1.5*exp(-2.6*f))
    @test c1 ≥ 1.8 # increase cognitive
    @test c2 ≤ 1.7 # decrease social
    # Exploitation state
    f, phase = 0.35, 2
    w, c1, c2 = PSO._adapt_parameters(rng, 1.75, 1.75, f, phase)
    @test w == 1 / (1 + 1.5*exp(-2.6*f))
    @test 1.775 ≤ c1 ≤ 1.8 # slightly increase cognitive
    @test 1.7 ≤ c2 ≤ 1.725 # slightly decrease social
    # Convergence state
    f, phase = 0.0, 3
    w, c1, c2 = PSO._adapt_parameters(rng, 1.75, 1.75, f, phase)
    @test w == 1 / (1 + 1.5*exp(-2.6*f))
    @test 1.775 ≤ c1 ≤ 1.8 # slightly increase cognitive
    @test 1.775 ≤ c2 ≤ 1.8 # slightly increase social
    # Jumping out state
    f, phase = 1.0, 4
    w, c1, c2 = PSO._adapt_parameters(rng, 1.75, 1.75, f, phase)
    @test w == 1 / (1 + 1.5*exp(-2.6*f))
    @test c1 ≤ 1.7 # decrease cognitive
    @test c2 ≥ 1.8 # increase social
end

for acceleration in (CPU1(), CPUProcesses(), CPUThreads())
    @testset "EvoTree Tuning with AdaptiveParticleSwarm and $(typeof(acceleration))" begin
        rng = StableRNG(123)
        features = rand(rng, 10_000) .* 5 .- 2
        X = MLJBase.table(reshape(features, (size(features)[1], 1)))
        y = sin.(features) .* 0.5 .+ 0.5
        y = EvoTrees.logit(y) + randn(rng, size(y))
        y = EvoTrees.sigmoid(y)

        tree = EvoTreeRegressor(rng=rng)
        r1 = range(tree, :max_depth; values=[3:7;])
        r2 = range(tree, :η; lower=-2, upper=0, scale=exp10)

        baseline_self_tuning_tree = TunedModel(
            model=tree,
            tuning=RandomSearch(rng=StableRNG(1234)),
            resampling=CV(nfolds=5, rng=StableRNG(8888)),
            range=[r1, r2],
            measure=mae,
            n=15,
            acceleration=acceleration
        )
        baseline_mach = machine(baseline_self_tuning_tree, X, y)
        fit!(baseline_mach, verbosity=0)
        baseline_rep = report(baseline_mach)
        baseline_best_loss = baseline_rep.best_history_entry.measurement[1]

        self_tuning_tree = TunedModel(
            model=tree,
            tuning=AdaptiveParticleSwarm(rng=StableRNG(1234)),
            resampling=CV(nfolds=5, rng=StableRNG(8888)),
            range=[r1, r2],
            measure=mae,
            n=15,
            acceleration=acceleration
        )
        mach = machine(self_tuning_tree, X, y)
        fit!(mach, verbosity=0)
        rep = report(mach)
        best_loss = rep.best_history_entry.measurement[1]

        # Compare with random search result with the same settings
        @test best_loss < baseline_best_loss ||
              isapprox(best_loss, baseline_best_loss; atol=1e-3)
    end
end


text:
"""
fit(type::UnfoldModel,f::FormulaTerm,tbl::DataFrame,data::Array{T,3},times)
fit(type::UnfoldModel,f::FormulaTerm,tbl::DataFrame,data::Array{T,2},basisfunction::BasisFunction)
fit(type::UnfoldModel,d::Dict,tbl::DataFrame,data::Array)

Generates Designmatrix & fits model, either mass-univariate (one model per epoched-timepoint) or time-expanded (modeling linear overlap).

# Examples
Mass Univariate Linear
```julia-repl
julia> data,evts = loadtestdata("testCase1")
julia> data_r = reshape(data,(1,:))
julia> data_e,times = Unfold.epoch(data=data_r,tbl=evts,τ=(-1.,1.9),sfreq=10) # cut the data into epochs. data_e is now ch x times x epoch

julia> f  = @formula 0~1+continuousA+continuousB # 1
julia> model = fit(UnfoldModel,f,evts,data_e,times)
```
Timexpanded Univariate Linear
```julia-repl
julia> basisfunction = firbasis(τ=(-1,1),sfreq=10,name="A")
julia> model = fit(UnfoldModel,Dict(Any=>(f,basisfunction),evts,data_r)
```

"""
function StatsModels.fit(
    UnfoldModelType::Type{T},
    f::FormulaTerm,
    tbl::DataFrame,
    data::AbstractArray,
    basisOrTimes::Union{BasisFunction,AbstractArray};
    kwargs...,
) where {T<:Union{<:UnfoldModel}}
    fit(UnfoldModelType, Dict(Any => (f, basisOrTimes)), tbl, data; kwargs...)
end

function StatsModels.fit(
    UnfoldModelType::Type{T},
    design::Dict,
    tbl::DataFrame,
    data::AbstractArray;
    kwargs...,
) where {T<:Union{<:UnfoldModel}}
    to = get_timer("Shared")
    if UnfoldModelType == UnfoldModel
        UnfoldModelType = designToModeltype(design)
    end
    uf = UnfoldModelType(design)

    @timeit to "designmatrix" designmatrix!(uf, tbl; kwargs...)
    @timeit to "fit" fit!(uf, data; kwargs...)

    return uf
end

function StatsModels.fit(
    UnfoldModelType::Type{T},
    X::DesignMatrix,
    data::AbstractArray;
    kwargs...,
) where {T<:Union{<:UnfoldModel}}
    if UnfoldModelType == UnfoldModel
        error(
            "Can't infer model automatically, specify with e.g. fit(UnfoldLinearModel...) instead of fit(UnfoldModel...)",
        )
    end
    uf = UnfoldModelType(Dict(), X)

    fit!(uf, data; kwargs...)

    return uf
end

isMixedModelFormula(f::ConstantTerm) = false
isMixedModelFormula(f::FormulaTerm) = isMixedModelFormula(f.rhs)

function isMixedModelFormula(f::Tuple)
    ix = [isa(t, FunctionTerm) for t in f]
    return any([isa(t.forig, typeof(|))|isa(t.forig,typeof(MixedModels.zerocorr)) for t in f[ix]])
end
function designToModeltype(design)
    # autoDetect
    tmp = collect(values(design))[1]
    f = tmp[1] # formula
    t = tmp[2] # Vector or BasisFunction

    isMixedModel = isMixedModelFormula(f)

    if typeof(t) <: BasisFunction
        if isMixedModel
            UnfoldModelType = UnfoldLinearMixedModelContinuousTime
        else
            UnfoldModelType = UnfoldLinearModelContinuousTime
        end
    else
        if isMixedModel
            UnfoldModelType = UnfoldLinearMixedModel
        else
            UnfoldModelType = UnfoldLinearModel
        end
    end
    return UnfoldModelType
end

# helper function for 1 channel data
function StatsModels.fit(
    UnfoldModelType::Type{T},
    design::Dict,
    tbl::DataFrame,
    data::AbstractVector,
    args...;
    kwargs...,
) where {T<:Union{<:UnfoldModel}}

    @debug("data array is size (X,), reshaping to (1,X)")
    data = reshape(data, 1, :)
    return fit(UnfoldModelType, design, tbl, data, args...; kwargs...)
end

"""
fit!(uf::UnfoldModel,data::Union{<:AbstractArray{T,2},<:AbstractArray{T,3}}) where {T<:Union{Missing, <:Number}}

Fit a DesignMatrix against a 2D/3D Array data along its last dimension
Data is typically interpreted as channel x time (with basisfunctions) or channel x time x epoch (for mass univariate)

Returns an UnfoldModel object

# Examples
```julia-repl
```

"""
function StatsModels.fit!(
    uf::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime},
    data::AbstractArray;
    kwargs...,
)
    # function content partially taken from MixedModels.jl bootstrap.jl
    df = Array{NamedTuple,1}()
    dataDim = length(size(data)) # surely there is a nicer way to get this but I dont know it

    Xs = modelmatrix(uf)
    # If we have3 dimension, we have a massive univariate linear mixed model for each timepoint
    if dataDim == 3
        firstData = data[1, 1, :]
        ntime = size(data, 2)
    else
        # with only 2 dimension, we run a single time-expanded linear mixed model per channel/voxel
        firstData = data[1, :]
        ntime = 1
    end
    nchan = size(data, 1)

    _,data = zeropad(Xs[1],data)
    # get a un-fitted mixed model object

    mm = LinearMixedModel_wrapper(formula(uf), firstData, Xs)

    # prepare some variables to be used
    βsc, θsc = similar(MixedModels.coef(mm)), similar(mm.θ) # pre allocate
    p, k = length(βsc), length(θsc)
    #β_names = (Symbol.(fixefnames(mm))..., )

    β_names = (Symbol.(vcat(fixefnames(mm)...))...,)
    β_names = (unique(β_names)...,)

    @assert(
        length(β_names) == length(βsc),
        "Beta-Names & coefficient length do not match. Did you provide two identical basis functions?"
    )

    @debug println("beta_names $β_names")
    @debug println("uniquelength: $(length(unique(β_names))) / $(length(β_names))")
    # for each channel
    prog = Progress(nchan * ntime, 0.1)
    #@showprogress .1
    for ch in range(1, stop = nchan)
        # for each time
        for t in range(1, stop = ntime)

            #@debug "ch:$ch/$nchan, t:$t/$ntime"
            @debug "data-size: $(size(data))"
            #@debug println("mixedModel: $(mm.feterms)")
            if ndims(data) == 3
                MixedModels.refit!(mm, data[ch, t, :])
            else
                MixedModels.refit!(mm, data[ch, :])
            end
            #@debug println(MixedModels.fixef!(βsc,mm))

            β = NamedTuple{β_names}(MixedModels.fixef!(βsc, mm))

            out = (
                objective = mm.objective,
                σ = mm.σ,
                β = NamedTuple{β_names}(MixedModels.fixef!(βsc, mm)),
                se = SVector{p,Float64}(MixedModels.stderror!(βsc, mm)), #SVector not necessary afaik, took over from MixedModels.jl
                θ = SVector{k,Float64}(MixedModels.getθ!(θsc, mm)),
                channel = ch,
                timeIX = ifelse(dataDim == 2, NaN, t),
            )
            push!(df, out)
            ProgressMeter.next!(prog; showvalues = [(:channel, ch), (:time, t)])
        end
    end

    uf.modelfit = UnfoldMixedModelFitCollection(
        df,
        deepcopy(mm.λ),
        getfield.(mm.reterms, :inds),
        copy(mm.optsum.lowerbd),
        NamedTuple{Symbol.(fnames(mm))}(map(t -> (t.cnames...,), mm.reterms)),
    )

    return uf.modelfit
end

function StatsModels.coef(
    uf::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime},
)
    beta = [x.β for x in MixedModels.tidyβ(modelfit(uf))]
    return reshape_lmm(uf, beta)
end

function MixedModels.ranef(
    uf::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime},
)
    sigma = [x.σ for x in MixedModels.tidyσs(modelfit(uf))]
    return reshape_lmm(uf, sigma)
end

function reshape_lmm(uf::UnfoldLinearMixedModel, est)
    ntime = length(collect(values(design(uf)))[1][2])
    nchan = modelfit(uf).fits[end].channel
    return permutedims(reshape(est, :, ntime, nchan), [3 2 1])
end
function reshape_lmm(uf::UnfoldLinearMixedModelContinuousTime, est)
    nchan = modelfit(uf).fits[end].channel
    return reshape(est, :, nchan)'

end

function StatsModels.fit!(
    uf::Union{UnfoldLinearModelContinuousTime,UnfoldLinearModel},
    data;
    solver = (x, y) -> solver_default(x, y),
    kwargs...,
)

    @assert ~isempty(designmatrix(uf))

    to = get_timer("Shared")

    @timeit to "modelmatrix" X = modelmatrix(uf)
    # mass univariate, data = ch x times x epochs
    X, data = zeropad(X, data)

    @debug "UnfoldLinearModel, datasize: $(size(data))"

    @timeit to "solver" uf.modelfit = solver(X, data)

end

LinearMixedModel_wrapper(form,data::Array{<:Union{TData},1},Xs;wts = []) where {TData<:Union{Missing,Number}}= @error("currently no support for missing values in MixedModels.jl")

"""
$(SIGNATURES)

Wrapper to generate a LinearMixedModel. Code taken from MixedModels.jl and slightly adapted.

"""
function LinearMixedModel_wrapper(
    form,
    data::Array{<:Union{TData},1},
    Xs;
    wts = [],
) where {TData<:Number}
    #    function LinearMixedModel_wrapper(form,data::Array{<:Union{Missing,TData},1},Xs;wts = []) where {TData<:Number}

    # XXX Push this to utilities zeropad
    # Make sure X & y are the same size
    m = size(Xs[1])[1]

    if m != size(data)[1]
        fe,data = zeropad(Xs[1],data)

        Xs = changeMatSize!(size(data)[1], fe, Xs[2:end])
    end

    y = (reshape(float(data), (:, 1)))

    MixedModels.LinearMixedModel(y, Xs, form, wts)
end

function MixedModels.LinearMixedModel(y, Xs, form::Array, wts)

    form_combined = form[1]
    for f in form[2:end]

        form_combined =
            form_combined.lhs ~
                MatrixTerm(form_combined.rhs[1] + f.rhs[1]) +
                form_combined.rhs[2:end] +
                f.rhs[2:end]
    end
    MixedModels.LinearMixedModel(y, Xs, form_combined, wts)
end

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
module NonNegLeastSquaresMLJInterface

using MLJModelInterface
using NonNegLeastSquares: nonneg_lsq

################################
########### Structure ##########
################################

MLJModelInterface.@mlj_model mutable struct NonNegativeLeastSquareRegressor <: MLJModelInterface.Deterministic
    alg::Symbol = :pivot::(_ in (:pivot, :nnls, :fnnls, :admm))
    variant::Symbol = :none::(_ in (:comb, :cache, :none))
    gram::Bool = false
    fit_intercept::Bool = true
end

################################
########### Methods ############
################################

"""
    augment_X(X, b)
Augment the matrix `X` with a column of ones if the intercept is to be fitted (`b=true`), return
`X` otherwise.
"""
function augment_X(X::Matrix, b::Bool)::Matrix
    b && return hcat(X, ones(eltype(X), size(X, 1), 1))
    return X
end

function MLJModelInterface.fit(m::NonNegativeLeastSquareRegressor, verbosity::Int, X, y)
    X = augment_X(MLJModelInterface.matrix(X), m.fit_intercept)
    fitresult = vec(nonneg_lsq(X, y;
                                alg=m.alg,
                                variant=m.variant,
                                gram=m.gram))
    cache = nothing
    report = NamedTuple{}()
    return (fitresult, cache, report)
end

function MLJModelInterface.fitted_params(model::NonNegativeLeastSquareRegressor, fitresult)
    return (coefs      = fitresult[1:end-Int(model.fit_intercept)],
            intercept = ifelse(model.fit_intercept, fitresult[end], nothing))
end

function MLJModelInterface.predict(m::NonNegativeLeastSquareRegressor, fitresult, Xnew)
    Xnew = augment_X(MLJModelInterface.matrix(Xnew), m.fit_intercept)
    ypred = Xnew * fitresult
    return ypred
end

################################
########### METADATA ###########
################################

const ALL_MODELS = Union{NonNegativeLeastSquareRegressor}

MLJModelInterface.metadata_pkg.(ALL_MODELS,
    name       = "NonNegLeastSquares",
    uuid       = "b7351bd1-99d9-5c5d-8786-f205a815c4d",
    url        = "https://github.com/ahwillia/NonNegLeastSquares.jl",
    julia      = true,
    license    = "MIT",
    is_wrapper = false,
)

MLJModelInterface.metadata_model(NonNegativeLeastSquareRegressor,
    input_scitype    = MLJModelInterface.Table(MLJModelInterface.Continuous),
    target_scitype   = AbstractVector{MLJModelInterface.Continuous},
    output_scitype   = MLJModelInterface.Table(MLJModelInterface.Continuous),
    supports_weights = false,
    descr            = "A non negative least square regressor model",
	load_path        = "NonNegLeastSquaresMLJInterface.NonNegativeLeastSquareRegressor"
    )

end


text:
# Implementation of Atomistic API with Molly.jl

include("molly/molly_base.jl")
include("molly/molly_simulator.jl")
include("molly/molly_result.jl")


text:
#!TODO: Yet to write more tests here

#=
Disclaimer: Many of the tests in this file taken from:
    https://github.com/tpapp/DynamicHMC.jl/blob/master/test/test-NUTS.jl
=#

############################################################################################
for iter in eachindex(objectives)
    _obj = objectives[iter]
    _flattentype = _obj.model.info.reconstruct.default.output
    @testset "Kernel construction and propagation, all models" begin
        mcmcdefault = MCMCDefault(;
        stepsize = ConfigStepsize(;
            stepsizeadaption=UpdateFalse()),
            GradientBackend = backends[1],
            generated = UpdateTrue()
        )
        mcmckernel = MCMC(
            _rng,
            kernel,
            _obj,
            mcmcdefault
        )
    end
end

############################################################################################
@testset "random booleans" begin
    @test rand_bool_logprob(_rng, 0)
    @test !rand_bool_logprob(_rng, log(0))
end


text:
export CUFFTError

struct CUFFTError <: Exception
    code::cufftResult
    msg::AbstractString
end
Base.show(io::IO, err::CUFFTError) = print(io, "CUFFTError(code $(err.code), $(err.msg))")

function CUFFTError(code::cufftResult)
    msg = status_message(code)
    return CUFFTError(code, msg)
end

function status_message(status)
    if status == CUFFT_SUCCESS
        return "the operation completed successfully"
    elseif status == CUFFT_INVALID_PLAN
        return "cuFFT was passed an invalid plan handle"
    elseif status == CUFFT_ALLOC_FAILED
        return "cuFFT failed to allocate GPU or CPU memory"
    elseif status == CUFFT_INVALID_TYPE
        return "cuFFT invalid type " # No longer used
    elseif status == CUFFT_INVALID_VALUE
        return "User specified an invalid pointer or parameter"
    elseif status == CUFFT_INTERNAL_ERROR
        return "Driver or internal cuFFT library error"
    elseif status == CUFFT_EXEC_FAILED
        return "Failed to execute an FFT on the GPU"
    elseif status == CUFFT_SETUP_FAILED
        return "The cuFFT library failed to initialize"
    elseif status == CUFFT_INVALID_SIZE
        return "User specified an invalid transform size"
    elseif status == CUFFT_UNALIGNED_DATA
        return "cuFFT unaligned data" # No longer used
    elseif status == CUFFT_INCOMPLETE_PARAMETER_LIST
        return "Missing parameters in call"
    elseif status == CUFFT_INVALID_DEVICE
        return "Execution of a plan was on different GPU than plan creation"
    elseif status == CUFFT_PARSE_ERROR
        return "Internal plan database error"
    elseif status == CUFFT_NO_WORKSPACE
        return "No workspace has been provided prior to plan execution"
    elseif status == CUFFT_NOT_IMPLEMENTED
        return "Function does not implement functionality for parameters given."
    elseif status == CUFFT_LICENSE_ERROR
        return "cuFFT license error" # Used in previous versions.
    elseif status == CUFFT_NOT_SUPPORTED
        return "Operation is not supported for parameters given."
    else
        return "unknown status"
    end
end

macro check(fft_func)
    quote
        local err::cufftResult
        err = $(esc(fft_func::Expr))
        if err != CUFFT_SUCCESS
            throw(CUFFTError(err))
        end
        err
    end
end


text:
# copy-free regular-expression searches on DenseStringViews, implemented
# by copying the low-level PCRE calls from julia/base/regex.jl

import Base.PCRE

function Base.occursin(r::Regex, s::DenseStringViewAndSub; offset::Integer=0)
    Base.compile(r)
    return PCRE.exec_r(r.regex, s, offset, r.match_options)
end

function Base.startswith(s::DenseStringViewAndSub, r::Regex)
    Base.compile(r)
    return PCRE.exec_r(r.regex, s, 0, r.match_options | PCRE.ANCHORED)
end

function Base.endswith(s::DenseStringViewAndSub, r::Regex)
    Base.compile(r)
    return PCRE.exec_r(r.regex, s, 0, r.match_options | PCRE.ENDANCHORED)
end

function Base.match(re::Regex, str::DenseStringViewAndSub, idx::Integer, add_opts::UInt32=UInt32(0))
    Base.compile(re)
    opts = re.match_options | add_opts
    matched, data = PCRE.exec_r_data(re.regex, str, idx-1, opts)
    if !matched
        PCRE.free_match_data(data)
        return nothing
    end
    n = div(PCRE.ovec_length(data), 2) - 1
    p = PCRE.ovec_ptr(data)
    mat = SubString(str, unsafe_load(p, 1)+1, prevind(str, unsafe_load(p, 2)+1))
    cap = Union{Nothing,SubString{String}}[unsafe_load(p,2i+1) == PCRE.UNSET ? nothing :
                                           SubString(str, unsafe_load(p,2i+1)+1,
                                           prevind(str, unsafe_load(p,2i+2)+1)) for i=1:n]
    off = Int[ unsafe_load(p,2i+1)+1 for i=1:n ]
    result = RegexMatch(mat, cap, unsafe_load(p,1)+1, off, re)
    PCRE.free_match_data(data)
    return result
end

Base.findnext(re::Regex, str::DenseStringViewAndSub, idx::Integer) = Base._findnext_re(re, str, idx, C_NULL)

function Base._findnext_re(re::Regex, str::DenseStringViewAndSub, idx::Integer, match_data::Ptr{Cvoid})
    if idx > nextind(str,lastindex(str))
        throw(BoundsError())
    end
    opts = re.match_options
    Base.compile(re)
    alloc = match_data == C_NULL
    if alloc
        matched, data = PCRE.exec_r_data(re.regex, str, idx-1, opts)
    else
        matched = PCRE.exec(re.regex, str, idx-1, opts, match_data)
        data = match_data
    end
    if matched
        p = PCRE.ovec_ptr(data)
        ans = (Int(unsafe_load(p,1))+1):prevind(str,Int(unsafe_load(p,2))+1)
    else
        ans = nothing
    end
    alloc && PCRE.free_match_data(data)
    return ans
end

# copied from Base.RegexMatchIterator
struct RegexMatchIterator{T<:DenseStringViewAndSub}
    regex::Regex
    string::T
    overlap::Bool
end
Base.compile(itr::RegexMatchIterator) = (compile(itr.regex); itr)
Base.eltype(::Type{<:RegexMatchIterator}) = RegexMatch
Base.IteratorSize(::Type{<:RegexMatchIterator}) = Base.SizeUnknown()

function Base.iterate(itr::RegexMatchIterator, (offset,prevempty)=(1,false))
    opts_nonempty = UInt32(PCRE.ANCHORED | PCRE.NOTEMPTY_ATSTART)
    while true
        mat = match(itr.regex, itr.string, offset,
                    prevempty ? opts_nonempty : UInt32(0))

        if mat === nothing
            if prevempty && offset <= sizeof(itr.string)
                offset = nextind(itr.string, offset)
                prevempty = false
                continue
            else
                break
            end
        else
            if itr.overlap
                if !isempty(mat.match)
                    offset = nextind(itr.string, mat.offset)
                else
                    offset = mat.offset
                end
            else
                offset = mat.offset + ncodeunits(mat.match)
            end
            return (mat, (offset, isempty(mat.match)))
        end
    end
    nothing
end

Base.eachmatch(re::Regex, str::DenseStringViewAndSub; overlap = false) =
    RegexMatchIterator(re, str, overlap)

# copied from julia/base/pcre.jl:
function PCRE.exec(re, subject::DenseStringViewAndSub, offset, options, match_data)
    rc = ccall((:pcre2_match_8, PCRE.PCRE_LIB), Cint,
                (Ptr{Cvoid}, Ptr{UInt8}, Csize_t, Csize_t, UInt32, Ptr{Cvoid}, Ptr{Cvoid}),
                re, subject, ncodeunits(subject), offset, options, match_data, PCRE.get_local_match_context())
    # rc == -1 means no match, -2 means partial match.
    rc < -2 && error("PCRE.exec error: $(err_message(rc))")
    return rc >= 0
end


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

function iterativealign(parameters::Parameters, cascades::AbstractVector{GammaCascade};
                        kwargs...)
    ungated = QuickReader(parameters);
    gated = QuickReader(parameters);
    iterativealign(ungated, gated, cascades; kwargs...)
end

function iterativealign(ungated::QuickReader, gated::QuickReader,
                        cascades::AbstractVector{GammaCascade};
                        nbins=5000, nbins2=3000, kwargs...)
    # TODO This is complicated with a lot of parts. Think more carefully
    # For each γ cascade, gate on the corresponding excitation energy window
    # And fit each γ decay
    # TODO Needs to perform some initial calibration so that the channel
    #      numbers aren't random
    default = GammaCalibrator([-913.5, 5.6])
    println("Reading ungated")
    read(ungated)
    ungateddetector = ChannelDetector(ungated)
    default(ungateddetector)
    signal = pyimport("scipy.signal")
    smoother(x) = signal.savgol_filter(x, 51, 4)
    chns = goodchannels(ungateddetector)
    # ungatedhist, edges = histogram(ungateddetector, nbins=nbins)
    # Contains the peaks for all gamma cascades for all spectra
    allpeaks = Vector{Vector{Float64}}[]
    for cascade in cascades
        println("Reading gated for $(cascade.gate)")
        read(gated, cascade.gate)
        gateddetector = ChannelDetector(gated)
        default(gateddetector)
        # gatedhist = histogram(gateddetector, edges)
        # Contains the peaks for all decays in a cascade
        # for all spectra
        peaks = Vector{Float64}[]
        individual = true
        for (e, region) in zip(cascade.e, cascade.regions)
            nbins = get(cascade.keywords, :bins1, 5000)
            nbins2 = get(cascade.keywords, :bins2, 3000)

            region = first(region):last(region)
            calib = localalign(ungateddetector, region; nbins=nbins, cascade.keywords[:align]...)
            aligned = calib(gateddetector, inplace=false)
            ax = plot(aligned, nbins=nbins)
            @show region
            plt.axvline(x=first(region), c="r")
            plt.axvline(x=last(region), c="r")

            total = sum(aligned)
            plot(total, nbins=nbins2)
            plt.axvline(x=first(region), c="k")
            plt.axvline(x=last(region), c="k")

            inv = invert(calib)
            peaks_ = if !individual
                peaks_ = localpeak(total, region, nbins=nbins2, smoother=smoother)
                peaks_ = inv(peaks_)[chns]
            else
                peaks_ = []
                for (detector, chn) in zip(aligned.data, 1:33)
                    chn ∉ chns && continue
                    peak = localpeak(detector, region, nbins=nbins2, smoother=smoother)
                    peak = inv[chn](peak)
                    push!(peaks_, peak)
                end
                peaks_
            end

            ax = plot(ungateddetector)
            ax.scatter(peaks_, 1.5:length(peaks_)+1, marker="x", color="r")

            push!(peaks, peaks_)
        end
        push!(allpeaks, peaks)
    end
    ax = plot(ungateddetector)
    for peaks in allpeaks
        for peaks_ in peaks
            ax.scatter(peaks_, 1.5:length(peaks_)+1, marker="x", color="r")
        end
    end
    allpeaks
end

function localpeak(detector::GammaDetector, region; nbins=3000, smoother=nothing)
    hist, edges = histogram(detector, nbins=nbins)
    index2e, e2index, e2indexint = transformations(edges)

    if !isnothing(smoother)
        hist = smoother(hist)
    end

    regioni = e2indexint(first(region)):e2indexint(last(region))
    peaki = findpeaks(hist[regioni], numpeaks=1)[1]
    #@show hist[regioni][peaki]
    peaki += first(regioni) |> Int
    #@show hist[peaki]
    #fig, ax = plt.subplots()
    #ax.step(edges[regioni], hist[regioni])
    #ax.plot(edges[regioni], smoothed[regioni])
    #ax.scatter(edges[peaki], hist[peaki])
    #plt.step(edges[regioni], hist[regioni])
    #plt.scatter(edges[peaki], hist[peaki])
    peak = index2e(peaki)
    #plt.scatter(peak, hist[peaki], marker="v", color="r", s=50)
    peak
end

function iterativealign(cascade::GammaCascade, ungatedhist,
                        edges, gated::QuickReader, peakwidth=50
                        ;kwargs...)
    println("Reading gated for $(cascade.gate)")
    read(gated, gate=cascade.gate)
    gateddetector = ChannelDetector(gated)
    gatedhist = histogram(gateddetector, edges)
    # Contains the peaks for all decays in a cascade
    # for all spectra
    peaks = Float64[]
    for e in cascade.e
        region = (e-peakwidth, e+peakwidth)
        peak = regionalign(ungatedhist, gatedhist,
                           edges, region;
                           kwargs...)
        push!(peaks, peak)
    end
    peaks
end

function transformations(range::AbstractRange)
    shift = first(range)
    gain = step(range)
    indextoe = Poly([shift, gain])
    etoindex = Poly([-shift/gain, 1/gain])
    etoindexint = Int∘round∘etoindex
    indextoe, etoindex, etoindexint
end

function transformations(detector::ChannelDetector, nbins=5000)
    hist, edges = histogram(detector, nbins=nbins)
    transformations(edges)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
# Stochastic programming constructs #
# ========================== #
"""
    WS(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenarioaData; optimizer = nothing)

Generate a **wait-and-see** (`WS`) model of the two-stage `stochasticprogram`, corresponding to `scenario`.

In other words, generate the first stage and the second stage of the `stochasticprogram` as if `scenario` is known to occur. Optionally, a capable `optimizer` can be supplied to `WS`.

See also: [`DEP`](@ref), [`EVP`](@ref)
"""
function WS(stochasticprogram::StochasticProgram{2}, scenario::AbstractScenario; optimizer = nothing)
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
    # Return WS model
    return _WS(generator(stochasticprogram,:stage_1),
               generator(stochasticprogram,:stage_2),
               stage_parameters(stochasticprogram, 1),
               stage_parameters(stochasticprogram, 2),
               scenario,
               Decisions(),
               Decisions(),
               optimizer)
end
function _WS(stage_one_generator::Function,
             stage_two_generator::Function,
             stage_one_params::Any,
             stage_two_params::Any,
             scenario::AbstractScenario,
             stage_one_decisions::Decisions,
             stage_two_decisions::Decisions,
             optimizer_constructor)
    ws_model = optimizer_constructor == nothing ? Model() : Model(optimizer_constructor)
    # Prepare decisions
    ws_model.ext[:stage_map] = Dict{MOI.VariableIndex, Int}()
    ws_model.ext[:decisions] = (stage_one_decisions, stage_two_decisions)
    add_decision_bridges!(ws_model)
    # Generate first stage and cache objective
    stage_one_generator(ws_model, stage_one_params)
    # Cache objective
    ws_obj = copy(objective_function(ws_model))
    ws_sense = objective_sense(ws_model)
    ws_sense = ws_sense == MOI.FEASIBILITY_SENSE ? MOI.MIN_SENSE : ws_sense
    # Generate second stage and finalize objective
    stage_two_generator(ws_model, stage_two_params, scenario)
    if ws_sense == objective_sense(ws_model)
        ws_obj += objective_function(ws_model)
    else
        ws_obj -= objective_function(ws_model)
    end
    set_objective_function(ws_model, ws_obj)
    set_objective_sense(ws_model, ws_sense)
    return ws_model
end
"""
    wait_and_see_decision(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenario, optimizer_constructor = nothing)

Calculate the optimizer of the **wait-and-see** (`WS`) model of the two-stage `stochasticprogram`, corresponding to `scenario`.

If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`WS`](@ref)
"""
function wait_and_see_decision(stochasticprogram::StochasticProgram{2}, scenario::AbstractScenario)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Solve WS model for supplied scenario
    ws_model = WS(stochasticprogram, scenario, optimizer = subproblem_optimizer(stochasticprogram))
    optimize!(ws_model)
    # Return WS decision
    return JuMP.value.(all_decision_variables(ws_model, 1))
end
"""
    EWS(stochasticprogram::StochasticProgram)

Calculate the **expected wait-and-see result** (`EWS`) of the `stochasticprogram`.

In other words, calculate the expectated result of all possible wait-and-see models, using the provided scenarios in `stochasticprogram`.

If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function EWS(stochasticprogram::StochasticProgram)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Solve all possible WS models and compute EWS
    return EWS(stochasticprogram, structure(stochasticprogram))
end
# Default implementation
function EWS(stochasticprogram::StochasticProgram, structure::AbstractStochasticStructure)
    return mapreduce(+, scenarios(stochasticprogram)) do scenario
        ws = _WS(stochasticprogram.generator[:stage_1],
                 stochasticprogram.generator[:stage_2],
                 stage_parameters(stochasticprogram, 1),
                 stage_parameters(stochasticprogram, 2),
                 scenario,
                 Decisions(),
                 Decisions(),
                 subproblem_optimizer(stochasticprogram))
        optimize!(ws)
        probability(scenario)*objective_value(ws)
    end
end
"""
    EWS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)

Approximately calculate the **expected wait-and-see result** (`EWS`) of the two-stage `stochasticmodel` to the current confidence level, over the scenario distribution induced by `sampler`.

The attribute [`NumEWSSamples`](@ref) is the size of the sampled models used to generate the interval and generally governs how tight it is. The confidence level can be set through the [`Confidence`](@ref) attribute.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function EWS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    N = MOI.get(stochasticmodel, NumEWSSamples())
    # Generate a sample model and statistically evaluate EWS
    let eval_model = instantiate(stochasticmodel, sampler, N; optimizer = optimizer)
        𝔼WS, σ² = statistical_EWS(eval_model, structure(eval_model))
        σ = sqrt(σ²)
        z = quantile(Normal(0,1), confidence)
        L = 𝔼WS - z * σ / sqrt(N)
        U = 𝔼WS + z * σ / sqrt(N)
        return ConfidenceInterval(L, U, confidence)
    end
end
# Default implementation
function statistical_EWS(stochasticprogram::StochasticProgram, structure::AbstractStochasticStructure)
    ws_models = map(scenarios(stochasticprogram)) do scenario
        ws = _WS(stochasticprogram.generator[:stage_1],
                 stochasticprogram.generator[:stage_2],
                 stage_parameters(stochasticprogram, 1),
                 stage_parameters(stochasticprogram, 2),
                 scenario,
                 Decisions(),
                 Decisions(),
                 subproblem_optimizer(stochasticprogram))
    end
    return welford(ws_models, probability.(scenarios(stochasticprogram)))
end
"""
    DEP(stochasticprogram::TwoStageStochasticProgram; optimizer = nothing)

Generate the **deterministically equivalent problem** (`DEP`) of the two-stage `stochasticprogram`, unless a cached version already exists.

In other words, generate the extended form the `stochasticprogram` as a single JuMP model. Optionally, a capable `optimizer` can be supplied to `DEP`.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function DEP(stochasticprogram::StochasticProgram{2}; optimizer = nothing)
    return DEP(stochasticprogram, structure(stochasticprogram); optimizer = optimizer)
end
function DEP(stochasticprogram::StochasticProgram{2}, ::AbstractStochasticStructure; optimizer = nothing)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache, :dep)
        dep = cache[:dep]
        optimizer != nothing && set_optimizer(dep, optimizer)
        return dep
    end
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
    # Generate and cache deterministic equivalent
    dep = StochasticStructure(scenario_types(stochasticprogram), Deterministic())
    generate!(stochasticprogram, dep)
    cache[:dep] = dep.model
    # Return DEP
    return dep
end
"""
    VRP(stochasticprogram::StochasticProgram)

Calculate the **value of the recouse problem** (`VRP`) in `stochasticprogram`.

In other words, optimize the stochastic program and return the optimal value.

If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`EVPI`](@ref), [`EWS`](@ref)
"""
function VRP(stochasticprogram::StochasticProgram)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Check if cached solution is available
    cache = solutioncache(stochasticprogram)
    if haskey(cache, :solution)
        # Return cached value
        return MOI.get(cache[:solution], MOI.ObjectiveValue())
    end
    # Solve DEP
    optimize!(stochasticprogram)
    # Return optimal value
    return objective_value(stochasticprogram)
end
"""
    VRP(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; confidence = 0.95)

Return a confidence interval around the **value of the recouse problem** (`VRP`) of `stochasticmodel` to the given `confidence` level.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`EVPI`](@ref), [`VSS`](@ref), [`EWS`](@ref)
"""
function VRP(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Optimize stochastic model using sample-based method
    optimize!(stochasticmodel, sampler)
    return objective_value(stochasticmodel)
end
"""
    EVPI(stochasticprogram::TwoStageStochasticProgram)

Calculate the **expected value of perfect information** (`EVPI`) of the two-stage `stochasticprogram`.

In other words, calculate the gap between `VRP` and `EWS`. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`VRP`](@ref), [`EWS`](@ref), [`VSS`](@ref)
"""
function EVPI(stochasticprogram::StochasticProgram{2})
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Calculate VRP
    vrp = VRP(stochasticprogram)
    # Solve all possible WS models and calculate EWS
    ews = EWS(stochasticprogram)
    # Return EVPI = EWS - VRP
    return abs(ews-vrp)
end
"""
    EVPI(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)

Approximately calculate the **expected value of perfect information** (`EVPI`) of the two-stage `stochasticmodel` to the current confidence level, over the scenario distribution induced by `sampler`.

In other words, calculate confidence intervals around `VRP` and `EWS`. If they do not overlap, the EVPI is statistically significant, and a confidence interval is calculated and returned.

The attribute [`NumSamples`](@ref) is the size of the sampled models used to generate the interval and generally governs how tight it is. The attribute [`NumLowerTrials`](@ref) is the number of sampled models used in the lower bound calculation and the attribute [`NumUpperTrials`](@ref) is the number of sampled models used in the upper bound calculation. The attribute [`NumEvalSamples`](@ref) is the size of the sampled models used in the upper bound calculation. The confidence level can be set through the [`Confidence`](@ref) attribute.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`VRP`](@ref), [`EWS`](@ref), [`VSS`](@ref)
"""
function EVPI(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    # Modify confidence for two-sided interval
    α = (1 - confidence) / 2
    MOI.set(stochasticmodel, Confidence(), 1 - α)
    # Calculate confidence interval around VRP
    optimize!(stochasticmodel, sampler)
    # Check return status
    status = termination_status(stochasticmodel)
    if !(status in AcceptableTermination)
        error("Stochastic model could not be solved to optimality, returned status $status.")
    end
    vrp = objective_value(stochasticmodel)
    # EWS solution of the corresponding size
    ews = EWS(stochasticmodel, sampler)
    # Restore confidence level
    MOI.set(stochasticmodel, Confidence(), confidence)
    # Check overlap
    if (upper(vrp) >= lower(ews) + eps() && upper(ews) >= lower(vrp) + eps()) ||
       (lower(vrp) <= upper(ews) - eps() && upper(vrp) >= lower(ews) + eps())
        @warn "EVPI is not statistically significant to the chosen confidence level and tolerance"
        # Return confidence interval as tolerance around zero
        tolerance = MOI.get(stochasticmodel, RelativeTolerance())
        return ConfidenceInterval(-tolerance, tolerance, confidence)
    end
    # Switch on sign
    if lower(ews) >= upper(vrp)
        # Return confidence interval around EVPI
        return ConfidenceInterval(lower(ews) - upper(vrp), upper(ews) - lower(vrp), confidence)
    else
        # Return confidence interval around EVPI
        return ConfidenceInterval(lower(vrp) - upper(ews), upper(vrp) - lower(ews), confidence)
    end
end
"""
    EVP(stochasticprogram::TwoStageStochasticProgram; optimizer = nothing)

Generate the **expected value problem** (`EVP`) of the two-stage `stochasticprogram`.

In other words, generate a wait-and-see model corresponding to the expected scenario over all available scenarios in `stochasticprogram`. Optionally, a capable `optimizer` can be supplied to `WS`.

See also: [`expected_value_decision`](@ref), [`EEV`](@ref), [`EV`](@ref), [`WS`](@ref)
"""
function EVP(stochasticprogram::StochasticProgram{2}; optimizer = nothing)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache,:evp)
        evp = cache[:evp]
        optimizer != nothing && set_optimizer(evp, optimizer)
        return evp
    end
    # Create EVP as a wait-and-see model of the expected scenario
    ev_model = WS(stochasticprogram, expected(stochasticprogram).scenario, optimizer = optimizer)
    # Cache EVP
    cache[:evp] = ev_model
    # Return EVP
    return ev_model
end
"""
    expected_value_decision(stochasticprogram::TwoStageStochasticProgram)

Calculate the optimizer of the `EVP` of the two-stage `stochasticprogram`.

If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`EVP`](@ref), [`EV`](@ref), [`EEV`](@ref)
"""
function expected_value_decision(stochasticprogram::StochasticProgram{2})
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Solve EVP
    evp = EVP(stochasticprogram, optimizer = master_optimizer(stochasticprogram))
    optimize!(evp)
    # Return EVP decision
    return JuMP.value.(all_decision_variables(evp, 1))
end
"""
    EV(stochasticprogram::TwoStageStochasticProgram)

Calculate the optimal value of the `EVP` of the two-stage `stochasticprogram`.

If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`EVP`](@ref), [`expected_value_decision`](@ref), [`EEV`](@ref)
"""
function EV(stochasticprogram::StochasticProgram{2})
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Solve EVP model
    evp = EVP(stochasticprogram, optimizer = master_optimizer(stochasticprogram))
    optimize!(evp)
    # Return optimal value
    return objective_value(evp)
end
"""
    EEV(stochasticprogram::TwoStageStochasticProgram)

Calculate the **expected value of the expected value solution** (`EEV`) of the two-stage `stochasticprogram`.

In other words, evaluate the `EVP` decision. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`EVP`](@ref), [`EV`](@ref)
"""
function EEV(stochasticprogram::StochasticProgram{2})
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Solve EVP model
    evp_decision = expected_value_decision(stochasticprogram)
    # Calculate EEV by evaluating the EVP decision
    eev = evaluate_decision(stochasticprogram, evp_decision)
    # Return EEV
    return eev
end
"""
    EEV(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)

Approximately calculate the **expected value of the expected value decision** (`EEV`) of the two-stage `stochasticmodel` to the current confidence level, over the scenario distribution induced by `sampler`.

The attribute [`NumEEVSamples`](@ref) is the size of the sampled models used in the eev calculation. The confidence level can be set through the [`Confidence`](@ref) attribute.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`EVP`](@ref), [`EV`](@ref)
"""
function EEV(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get the instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    N = MOI.get(stochasticmodel, NumSamples())
    # Generate expected value decision
    sp = instantiate(stochasticmodel, sampler, N, optimizer = optimizer)
    x̄ = expected_value_decision(sp)
    # Evaluate expected value decision
    M = MOI.get(stochasticmodel, NumEvalSamples())
    MOI.set(stochasticmodel, NumEvalSamples(), MOI.get(stochasticmodel, NumEEVSamples()))
    eev = evaluate_decision(stochasticmodel, x̄, sampler)
    # Restore NumEvalSamples()
    MOI.set(stochasticmodel, NumEvalSamples(), M)
    # Return EEV
    return eev
end
"""
    VSS(stochasticprogram::TwoStageStochasticProgram)

Calculate the **value of the stochastic solution** (`VSS`) of the two-stage `stochasticprogram`.

In other words, calculate the gap between `EEV` and `VRP`. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function VSS(stochasticprogram::StochasticProgram{2})
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Solve EVP and determine EEV
    eev = EEV(stochasticprogram)
    # Calculate VRP
    vrp = VRP(stochasticprogram)
    # Return VSS = VRP-EEV
    return abs(vrp-eev)
end
"""
    VSS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)

Approximately calculate the **value of the stochastic solution** (`VSS`) of the two-stage `stochasticmodel` to the current confidence level, over the scenario distribution induced by `sampler`.

In other words, calculate confidence intervals around `EEV` and `VRP`. If they do not overlap, the VSS is statistically significant, and a confidence interval is calculated and returned. `Ñ` is the number of samples in the out-of-sample evaluation of EEV.

The attribute [`NumSamples`](@ref) is the size of the sampled models used to generate the interval and generally governs how tight it is. The same size is used to generate the expected value decision. The attribute [`NumLowerTrials`](@ref) is the number of sampled models used in the lower bound calculation and the attribute [`NumUpperTrials`](@ref) is the number of sampled models used in the upper bound calculation. The attribute [`NumEvalSamples`](@ref) is the size of the sampled models used in the upper bound calculation and the attribute [`NumEEVSamples`] is the size of the sampled models used in the `EEV` calculation. The confidence level can be set through the [`Confidence`](@ref) attribute.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`VRP`](@ref), [`EEV`](@ref), [`EVPI`](@ref)
"""
function VSS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; confidence::AbstractFloat = 0.95, Ñ::Integer = 1000, tol::AbstractFloat = 1e-1, kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    # Modify confidence for two-sided interval
    α = (1-confidence)/2
    MOI.set(stochasticmodel, Confidence(), 1 - α)
    # Calculate confidence interval around VRP
    optimize!(stochasticmodel, sampler)
    # Check return status
    status = termination_status(stochasticmodel)
    if !(status in AcceptableTermination)
        error("Stochastic model could not be solved to optimality, returned status $status.")
    end
    vrp = objective_value(stochasticmodel)
    # Calculate confidence interval around EEV
    eev = EEV(stochasticmodel, sampler)
    # Restore confidence level
    MOI.set(stochasticmodel, Confidence(), confidence)
    # Check overlap
    if (upper(eev) >= lower(vrp) + eps() && upper(vrp) >= lower(eev) + eps()) ||
       (lower(eev) <= upper(vrp) - eps() && upper(eev) >= lower(vrp) + eps())
        @warn "VSS is not statistically significant to the chosen confidence level and tolerance"
        # Return confidence interval as tolerance around zero
        tolerance = MOI.get(stochasticmodel, RelativeTolerance())
        return ConfidenceInterval(-tolerance, tolerance, confidence)
    end
    # Switch on sign
    if lower(vrp) >= upper(eev)
        # Return confidence interval around VSS
        return ConfidenceInterval(lower(vrp) - upper(eev), upper(vrp) - lower(eev), confidence)
    else
        # Return confidence interval around VSS
        return ConfidenceInterval(lower(eev) - upper(vrp), upper(eev) - lower(vrp), confidence)
    end
end
# ========================== #


text:
# N-gons

#Makie

#Pkg.add("Makie")

using CairoMakie, Makie.GeometryBasics
CairoMakie.activate!()
let
    n = 20
    function ngonShape(h, k, r, n)
        Polygon([Point2f(h .+ r * sin.(m * 2π / n), k .+ r * cos.(m * 2π / n)) for m in 1:n])
    end
    polysCentric = [ngonShape(0, 0, 3 / i^1.5, i) for i in 3:n]
    polysCircular = [ngonShape(√2 / 2 * sin(θ), √2 / 2 * cos(θ), 0.15 / √idx, idx + 2)
                     for (idx, θ) in enumerate(LinRange(0, 2π * (1 - 1 / (n - 2)), n - 2))]
    cmap = to_colormap(:gist_ncar)[3:end]

    with_theme(theme_light()) do
        fig, ax, = poly(polysCentric; color = 1:n-2, colormap = cmap,
            axis = (; aspect = DataAspect()), figure = (; resolution = (600, 400)))
        poly!(polysCircular; color = 1:n-2, colormap = cmap)
        hidedecorations!(ax; grid = false)
        hidespines!(ax)
        display(fig)
    end
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
function can_solve_ut(A::SMat{T}, g::SRow{T}) where T <: Union{FieldElem, nmod}
  # Works also for non-square matrices
  #@hassert :HNF 1  ncols(A) == nrows(A)
  @hassert :HNF 2  isupper_triangular(A)
  # assumes A is upper triangular, reduces g modulo A to zero and collects
  # the transformation
  # supposed to be over a field...

  sol = sparse_row(base_ring(g))

  while length(g)>0
    s = g.pos[1]
    j = 1
    while j<= nrows(A) && A.rows[j].pos[1] < s
      j += 1
    end
    if j > nrows(A) || A.rows[j].pos[1] > s
      break
    end
    @hassert :HNF 2  A.rows[j].pos[1] == g.pos[1]
    p = divexact(g.values[1], A.rows[j].values[1])
    push!(sol.pos, j)
    push!(sol.values, p)
    _g = Hecke.add_scaled_row(A[j], g, -p)
    @assert _g != g
    g = _g
    @hassert :HNF 2  length(g)==0 || g.pos[1] > A[j].pos[1]
  end
  if length(g) == 0
    return true, sol
  else
    return false, sol
  end

  if length(g) > 0
    li = inv(g.values[1])
    for i=1:length(g)
      g.values[i] *= li
    end
  end
  return sol
end

function solve_ut(A::SMat{T}, g::SRow{T}) where T <: Union{FieldElem, nmod}
  fl, sol = can_solve_ut(A, g)
  @assert fl
  return sol
end

#TODO: write vector reconstruction and use it here.
@doc Markdown.doc"""
    rational_reconstruction(A::SRow{fmpz}, M::fmpz) -> Bool, SRow{fmpz}, fmpz

Apply rational reconstruction to the entries of $A$. Returns true iff
successful. In this case, the numerator is returned as a matrix and the
common denominator in the third value.
"""
function rational_reconstruction(A::SRow{fmpz}, M::fmpz)
  B = sparse_row(FlintZZ)
  de = fmpz(1)
  M2 = div(M, 2)
  nbM = div(nbits(M), 2)
  fl, d,n = true, fmpz(1), fmpz(1)
  for (p,v) = A
    vv = mod_sym(d*v, M)
    if nbits(vv) < nbM
      fl, n = true, vv
    else
      fl, n, d = rational_reconstruction(v, M)
    end
    if !fl
      return false, B, de
    end
    if de % d == 0
      push!(B.pos, p)
      push!(B.values, n*(div(de, d)))
    else
      l = lcm(d, de)
      B = div(l, de) * B
      push!(B.pos, p)
      push!(B.values, div(l, d)*n)
      de = l
    end
  end
  return true, B, de
end

function solve_ut(A::SMat{fmpz}, b::SRow{fmpz})
  @hassert :HNF 1  isupper_triangular(A)
  #still assuming A to be upper-triag

  sol = sparse_row(FlintZZ)
  den = fmpz(1)
  while length(b) > 0
    p = b.pos[1]
    v = b.values[1]
    d = gcd(A[p].values[1], v)
    if d != A[p].values[1]  # A does not divide b, the RHSS
      f = div(A[p].values[1], d)
      sol = f * sol
      b = f * b
      v *= f
      den *= f
    end
    push!(sol.pos, p)
    @hassert :HNF 2  mod(v, A[p].values[1]) == 0
    push!(sol.values, div(v, A[p].values[1]))
    b = b - sol.values[end]*A[p]
    @hassert :HNF 2  length(b) == 0 || b.pos[1] > p
  end
  return sol, den
end

function solve_ut(A::SMat{fmpz}, b::SMat{fmpz})
  @hassert :HNF 1  isupper_triangular(A)
  #still assuming A to be upper-triag
  d = fmpz(1)
  r = sparse_matrix(FlintZZ)
  for i = b
    x, dx = solve_ut(A, i)
    nd = lcm(d, dx)
    if nd != d
      dd = div(nd, d)
      for j = 1:r.r
        Hecke.scale_row!(r, j, dd)
      end
      d = nd
    end
    push!(r, div(d, dx)*x)
  end
  return r, d
end

################################################################################
#
#  Determinant
#
################################################################################

@doc Markdown.doc"""
    det_mc(A::SMat{fmpz})

Computes the determinant of $A$ using a LasVegas style algorithm,
i.e. the result is not proven to be correct.
Uses the dense (nmod_mat) determinant on $A$ for various primes $p$.
"""
function det_mc(A::SMat{fmpz})

  @hassert :HNF 1  A.r == A.c
  if isupper_triangular(A)
    z = fmpz[ A[i, i] for i in 1:A.r]
    return prod(z)
  end

  b = sparse_matrix(matrix(FlintZZ, 1, A.c, rand(1:10, A.c)))
  _, qq = solve_dixon_sf(A, b)

  q = p_start # global prime
  first = true
  dd = fmpz(1)
  mm = fmpz(1)
  last = fmpz(0)
  while true
    R = ResidueRing(FlintZZ, q, cached = false)
    d = det(matrix(change_base_ring(R, A)))*inv(R(qq))
    if first
      dd = fmpz(d)
      mm = fmpz(q)
      first = false
    else
      dd = crt(dd, mm, fmpz(d), fmpz(q), true)
      mm *= q
    end
    q = next_prime(q)
    if dd == last
      return dd*qq
    else
      last = dd
    end
  end
end

@doc Markdown.doc"""
    det(A::SMat{fmpz})

The determinant of $A$ using a modular algorithm.
Uses the dense (nmod_mat) determinant on $A$ for various primes $p$.
"""
function det(A::SMat{fmpz})
  @hassert :HNF 1  A.r == A.c
  if isupper_triangular(A)
    return prod(fmpz[A[i,i] for i=1:A.r]) end

  b = div(nbits(hadamard_bound2(A)), 2)
  lp = fmpz[p_start]
  while b > 0
    push!(lp, next_prime(lp[end]))
    b -= nbits(lp[end])
  end

  #TODO: re-use the nmod_mat....
  ld = fmpz[]
  for q in lp
    R = ResidueRing(FlintZZ, Int(q), cached = false)
    push!(ld, fmpz(det(matrix(change_base_ring(R, A)))))
  end
  #ld = [fmpz(det(matrix(sparse_matrix(A, Int(q))))) for q = lp]
  return crt_signed(ld, crt_env(lp))
end

@doc Markdown.doc"""
    echelon_with_transform(A::SMat{nmod}) -> SMat, SMat

Find a unimodular matrix $T$ and an upper-triangular $E$ s.th.
$TA = E$ holds.
"""
function echelon_with_transform(A::SMat{nmod})
  z = hcat(A, identity_matrix(SMat, base_ring(A), A.r))
  M = Hecke.ModuleCtxNmod(base_ring(A), z.c)
  for i=z
    Hecke.add_gen!(M, i)
  end
  return sub(M.basis, 1:A.r, 1:A.c), sub(M.basis, 1:A.r, A.c+1:z.c)
end

@doc Markdown.doc"""
    solve_dixon_sf(A::SMat{fmpz}, b::SRow{fmpz}, is_int::Bool = false) -> SRow{fmpz}, fmpz
    solve_dixon_sf(A::SMat{fmpz}, B::SMat{fmpz}, is_int::Bool = false) -> SMat{fmpz}, fmpz

For a sparse square matrix $A$ of full rank and a sparse matrix (row), find
a sparse matrix (row) $x$ and an integer $d$ s.th.
$$x A = bd$$
holds.
The algorithm is a Dixon-based linear p-adic lifting method.
If \code{is_int} is given, then $d$ is assumed to be $1$. In this case
rational reconstruction is avoided.
"""
function solve_dixon_sf(A::SMat{fmpz}, b::SRow{fmpz}, is_int::Bool = false)
  B = sparse_matrix(FlintZZ)
  push!(B, b)
  s, d = solve_dixon_sf(A, B, is_int)
  return s[1], d
end

function solve_dixon_sf(A::SMat{fmpz}, B::SMat{fmpz}, is_int::Bool = false)
  #for square matrices (s) of full rank (f) only.
  p = next_prime(2^20)
  R = ResidueRing(FlintZZ, p, cached = false)

  Ap = change_base_ring(R, A)

  #want AT = upper_triag.
  #Let J = anti-identity, so JA inverts the rows of A and AJ the columns
  #
  # B := (JA)' = A'J and TB = E in upper triag
  # J (TB)' J = JE'J is still upper triag and
  # J (TB)' J = J B' T' J = J (JA) T' J = A (JT)' = JE'J still upp-triag

  Bp = copy(Ap)
  reverse_rows!(Bp)
  Bp = transpose(Bp)
  Ep, Tp = echelon_with_transform(Bp)
  @hassert :HNF 1  Ep.c == Ep.r
#  @hassert :HNF 1  nmod_mat(Tp) * nmod_mat(Bp) == nmod_mat(Ep)

  reverse_rows!(Ep)
  Ep = transpose(Ep)
  reverse_rows!(Ep)
#  @hassert :HNF 1  Hecke.isupper_triangular(Ep)

  reverse_rows!(Tp)
  Tp = transpose(Tp)
#  @hassert :HNF 1  nmod_mat(Ap)*nmod_mat(Tp) == nmod_mat(Ep)

  #now, to solve xA = b, we do
  #              xAT = bT since AT is upper-triag, we can do this!

  sol_all = sparse_matrix(FlintZZ)
  den_all = fmpz(1)

  for b in B
    pp = fmpz(1)
    b_orig = b

    bp = change_base_ring(R, b)

    sol = sparse_row(FlintZZ)
    last = (sol, fmpz(1))

    while true
      bp = mul(bp, Tp)
      zp = solve_ut(Ep, bp)
      z = lift(zp)

      sol += pp*z

      pp *= fmpz(p)

  #    @hassert :HNF 1  iszero(SRow(b_orig - Hecke.mul(sol, A), pp))

      if is_int
        fl = true
        nu = copy(sol)
        mod_sym!(nu, pp)
        de = fmpz(1)
      else
        fl, nu, de = rational_reconstruction(sol, pp)
      end
      if fl
  #      @hassert :HNF 1  SRow(de*sol, pp) == SRow(nu, pp)
  #      @hassert :HNF 1  SRow(mul(nu, A), pp) == SRow(de*b_orig, pp)
        if last == (nu, de)
          if mul(nu, A) == de*b_orig
            l = lcm(den_all, de)
            if l == den_all
              push!(sol_all, div(l, de)*nu)
            else
              sol_all = div(l, den_all)*sol_all
              push!(sol_all, div(l, de)*nu)
              den_all = l
            end
            break
          end
          println("bad")
        else
          last = (nu, de)
        end
      end

  #    @hassert :HNF 1  SRow(Hecke.mul(z, A), p) == bp
      b = b - mul(z, A)

      for i=1:length(b.values)
  #      @hassert :HNF 1  b.values[i] % p == 0
        b.values[i] = div(b.values[i], p)
      end
      bp = change_base_ring(R, b)
    end
  end
  return sol_all, den_all
end

function echelon!(S::SMat{T}; complete::Bool = false) where T <: FieldElem
  i = 1
  while i <= nrows(S)
    m = ncols(S)+1
    mp = 0
    for j=i:nrows(S)
      if m > S[j].pos[1]
        m = S[j].pos[1]
        mp = j
      end
    end
    if mp == 0
      return
    end
    if mp != i
      swap_rows!(S, i, mp)
    end
    Si = -inv(S[i].values[1])
    j = i+1
    while j <= nrows(S)
      if S[j].pos[1] == m
        add_scaled_row!(S, i, j, S[j].values[1]*Si)
        if length(S[j].values) == 0
          deleteat!(S.rows, j)
          S.r -= 1
        else
          j += 1
        end
      else
        j += 1
      end
    end
    i += 1
  end
  if complete
    for i = nrows(S):-1:2
      p = S[i].pos[1]
      c = S[i,p]
      if !isone(c)
        scale_row!(S, i, inv(c))
      end
      for j=i-1:-1:1
        v = S[j, p]
        if !iszero(v)
          add_scaled_row!(S, i, j, -v)
        end
      end
    end
  end
end

function solve(a::SMat{T}, b::SRow{T}) where T <: FieldElem
  fl, sol = can_solve_with_solution(a, b)
  @assert fl
  return sol
end

function can_solve_with_solution(a::SMat{T}, b::SRow{T}) where T <: FieldElem
  c = hcat(a, identity_matrix(SMat, base_ring(a), a.r))
  echelon!(c)
  fl, sol = can_solve_ut(sub(c, 1:nrows(c), 1:a.c), b)
  if fl
    return fl, mul(sol, sub(c, 1:nrows(c), a.c+1:c.c))
  else
    return fl, sol
  end
end

#TODO: can_solve using Dixon for Q, NF
#      for SMat rather than SRow only