text:
using Distributed

#addprocs before including LAIprocessing.jl !
# 9 processors ideal because typical size of image set
addprocs(min(Base.Sys.CPU_THREADS, 9) - nprocs())

using Logging, Dates

# Setup Logging
isdir("logs") || mkdir("logs")
logger_io = open(joinpath("logs", "logjulia_$(Dates.format(Dates.now(), "dd-mm-yyyy_HHhMM")).log"), "a")
logger = SimpleLogger(logger_io, Logging.Debug)
global_logger(logger)
TEMPSETLOG = joinpath("logs", "tempsetlog.log")
DATALOG = joinpath("logs", "data.txt")

#get version of LeafAreaIndex.jl
# using Pkg
# import LibGit2
# LAICOMMIT = LibGit2.readchomp(`rev-parse HEAD`, dir=Pkg.dir("LeafAreaIndex"))

using PyCall, DataFrames
# use Conda.add("pyodbc") if not yet installed
@pyimport pyodbc

include("LAIprocessing.jl")

# check available memory
#import Humanize
#mem() = Humanize.datasize(1000*parse(Int, split(readall(`wmic os get FreePhysicalMemory`))[2]))

# convenience functions for quering the database
function selectcolnames(cursor::PyObject, tablename)
    sql = "SELECT column_name FROM information_schema.columns WHERE table_name = '$tablename'"
    cex = cursor[:execute](sql)
    pytable = cex[:fetchall]()
    String[collect(obj)[1] for obj in pytable]
end
function selecttable(cursor::PyObject, tablename, where::String, justone::Bool, ascdesc::String)
    # if justone: only select the first valid row
    if justone
        cex = cursor[:execute]("SELECT top 1 * FROM $tablename WHERE $where ORDER BY id $ascdesc")
    else
        cex = cursor[:execute]("SELECT * FROM $tablename WHERE $where ORDER BY id $ascdesc")
    end
    pytable = cex[:fetchall]()
    res = map(collect, pytable)
    df = DataFrame()
    # convert string to Symbol for DataFrame indexing
    colnames = [Symbol(c) for c in selectcolnames(cursor, tablename)]
    for col in eachindex(colnames)
        colvals = [res[row][col] for row in eachindex(res)]
        df[!, colnames[col]] = colvals
    end
    df
end
function selectimages(cursor::PyObject,where::String)
    cex = cursor[:execute]("select path, plotLocations.slope, plotLocations.slopeAspect, images.ID from images left join [dbo].[plotLocations] on images.plotLocationID = [dbo].[plotLocations].ID WHERE $where  order by images.ID ASC")
    pytable = cex[:fetchall]()
    res = map(collect, pytable)
    df = DataFrame()
    # convert string to Symbol for DataFrame indexing
    colnames = ["path","slope","slopeAspect", "ID"]
    for col in eachindex(colnames)
        colvals = [res[row][col] for row in eachindex(res)]
        df[!, colnames[col]] = colvals
    end
    df
end
function updatetable(conn::PyObject, tablename, ID::Int, columnname::Symbol,newvalue)
    cursor = conn[:cursor]()
    sql = "UPDATE $tablename SET $columnname = '$(newvalue)' WHERE ID = $ID"
    cex = cursor[:execute](sql)
    conn[:commit]()
    nothing
end

function process_calibration(conn)

    cursor = conn[:cursor]()
    cameraSetup = selecttable(cursor, :cameraSetup, " processed = 0 and pathCenter is not null ", true, "ASC")

    size(cameraSetup, 1) != 0 || return
    setupID = cameraSetup[1, :ID]
    try
        @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - detected new processed=false in cameraSetup table"

        #setupID    = cameraSetup[1, :ID]
        pathCenter = cameraSetup[1, :pathCenter]
        pathProj   = cameraSetup[1, :pathProj]
        width      = cameraSetup[1, :width]
        height     = cameraSetup[1, :height]

        df = readtable(pathCenter, names=[:x, :y, :circle])
        lensx, lensy = processcenterfile(df, height, width, TEMPSETLOG)
        @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - result x: $lensx y:$lensy"
        updatetable(conn, "cameraSetup", setupID, :x, lensx)
        updatetable(conn, "cameraSetup", setupID, :y, lensy)

        df = readtable(pathProj, names=[:cm, :px, :H, :pos])
        lensa, lensb = processprojfile(df, height, width, TEMPSETLOG)
        @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - result a: $lensa b:$lensb"
        updatetable(conn, "cameraSetup", setupID, :a, lensa)
        updatetable(conn, "cameraSetup", setupID, :b, lensb)

        updatetable(conn, "cameraSetup", setupID, :processed, 1)
    catch y
        error("Could not process center calibration: $setupID with error $y")
    finally
        updatetable(conn, "cameraSetup", setupID, :processed, 1)
    end
end

function gettabledata(cursor, results)

    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - reading LAI_App database tables"

    plotSetID = results[1, :plotSetID];
    @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - plotSetID = $plotSetID"
    plotSet = selecttable(cursor, :plotSets, "ID = $plotSetID", true, "ASC")
    if size(plotSet)[1] == 0
        error("Could not find plotSetID $plotSetID from results table in plotSets table.")
        set_processed()
        return
    end

    uploadSetID = plotSet[1, :uploadSetID]
    @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - uploadSetID = $uploadSetID"
    uploadSet = selecttable(cursor, :uploadSet, "ID = $uploadSetID", true, "ASC")
    if size(uploadSet)[1] == 0
        error("Could not find uploadSetID $uploadSetID from results table in uploadSet table.")
        set_processed()
        return
    end

    plotID = plotSet[1, :plotID]
    @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - plotID = $plotID"
    plot = selecttable(cursor, :plots, "ID = $plotID", true, "ASC")
    if size(plot)[1] == 0
        err("Could not find plotID $plotID from plots table in uploadSet table.")
        set_processed()
        return
    end

    camSetupID = uploadSet[1, :camSetupID]
    @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - camSetupID = $camSetupID"
    cameraSetup = selecttable(cursor, :cameraSetup, "ID = $camSetupID", true,"ASC")
    if size(cameraSetup)[1] == 0
        error("Could not find camSetupID $camSetupID from uploadSet table in cameraSetup table.")
        set_processed()
        return
    end

    return uploadSet, uploadSetID, plot, plotSetID, cameraSetup
end

function process_images(conn)
    cursor = conn[:cursor]()
    results = selecttable(cursor, "results", "processed = 0", true,"DESC")
    size(results)[1] != 0 || return
    resultsID = results[1, :ID]
    set_processed() = updatetable(conn, "results", resultsID, :processed, 1)
    # TODO make macro for all these `if not found return` statements
    try
        @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - detected new processed=false in results table"

        uploadSet, uploadSetID, plot, plotSetID, cameraSetup = gettabledata(cursor, results)

        lensx = cameraSetup[1, :x]
        lensy = cameraSetup[1, :y]
        lensa = cameraSetup[1, :a]
        lensb = cameraSetup[1, :b]
        lensρ = cameraSetup[1, :maxRadius]
        lensparams = (lensx, lensy, lensa, lensb, lensρ)
        @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - lens parameters: $lensparams"

        images = selectimages(cursor, "plotSetID = $plotSetID")
        @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - start images processing"
        success = false
        LAIres = Dict()
        try
            LAIres = processimages(images,lensparams,TEMPSETLOG,DATALOG)
            success = LAIres["success"]
            @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - uploadset $uploadSetID process completed with success: $success"
        catch y
            error("$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - Could not process uploadset: $uploadSetID with error $y")
        end
        println("$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - uploadset $uploadSetID process completed with success: $success")
        updatetable(conn, "results", resultsID, :processed, 1)
        updatetable(conn, "results", resultsID, :succes, ifelse(success,1,0))
        updatetable(conn, "results", resultsID, :resultLog, string(read(open(TEMPSETLOG), String)))
        datafile = open(DATALOG)
        updatetable(conn, "results", resultsID, :data, read(datafile, String))
        close(datafile)
        # todo
        LAICOMMIT = "739b1af9f496b6fd919ee09088e5210287cd9e56"
        updatetable(conn, "results", resultsID, :scriptVersion, LAICOMMIT)
        if success
            LAIvalue = LAIres["LAI"]
            LAIsd    = LAIres["LAIsd"]
            try
                println("LAIvalue: $LAIvalue")
                if !(isnan(LAIvalue))
                    updatetable(conn, "results", resultsID, :LAI, LAIvalue)
                end
                @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - added LAI to results table for ID $resultsID"
                println("added LAI to results table for ID $resultsID")
                if !(isnan(LAIsd))
                    updatetable(conn, "results", resultsID, :LAI_SD, LAIsd)
                end
                @info "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - added LAI_SD to results table for ID $resultsID"
                for (reskey, col) in [
                    ("csv_gapfraction", :gapfraction), ("csv_exif", :exif),
                    ("csv_histogram", :histogram), ("csv_stats", :stats),
                    ("jpgpath", :jpgPath), ("binpath", :binPath), ("LAIs", :LAI),
                    ("LAIe", :LAIe), ("threshold", :threshold), ("clumping", :clumping),
                    ("overexposure", :overexposure),("ALIA", :ALIA)]
                    for (imgp, csv) in LAIres[reskey]
                        imgp in images.path || continue
                        imageID = images.ID[images.path.==imgp][1]
                        if (!( col == :LAI && isnan(csv)) && !( col == :LAIe && isnan(csv)))
                            updatetable(conn, "images", imageID , col, csv)
                        end
                    end
                end
                for im in eachrow(images)
                    imID = im.ID
                    slope = im.slope
                    slopeaspect = im.slopeAspect
                    if (slope == nothing || slope == zero(slope) )
                        updatetable(conn, "images", imID , :slope, 0)
                        updatetable(conn, "images", imID , :slopeAspect, 360)
                    else
                        updatetable(conn, "images", imID , :slope, slope)
                        updatetable(conn, "images", imID , :slopeAspect, slopeaspect)
                    end
                end
            catch y
                error("$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - could not add LAI to results table, error: $y")
            end
        end
    catch y
        error("$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - caugth general error in interior loop: $y")
    finally
        updatetable(conn, "results", resultsID, :processed, 1)
    end
end

# Main loop
function mainloop(conn)
    println("Started Leaf Area Index Service")
    while true
        try
            #process_calibration(conn)
            process_images(conn)
        catch y
            error("$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - caugth general error in interior loop: $y")
        end
        flush(logger_io)
        sleep(1)
    end
end

# Connect to the database

cnxn  = pyodbc.connect("DSN=LAI")
mainloop(cnxn)

max_stars_repo_name: ETC-UA/LAIscript

max_issues_repo_name: ETC-UA/LAIscript

max_forks_repo_name: ETC-UA/LAIscript


text:
using DFTK
using Aqua
Aqua.test_all(DFTK, ambiguities=false, stale_deps=(ignore=[:Primes, ], ))

ext: jl


text:
module omegasim  #=

is the module definition of the Julia package "omegasim".

omegasim aims at practical simulations of composite industrial systems. A composite system consists of two or more subsystems characterized by different response time scales. In omegasim, different simulation methods are used for subsystems with different response time scales, and inputs and outputs of those subsystems are interconnected for the whole system simulation. This feature is quite reasonable for the simulation of industrial systems which are inherently composed of two or more subsystems with different response time scales. Currently, omegasim has the following two time-domain simulation methods.

1. Transient simulation method.
2. Phasor simulation method.

Before talking about these methods, let us define scientific terms. When a system does not include distributed-parameter elements, then it can be fully described by a set of differential and algebraic equations, or DAEs for short. Differential equations describe the behavior of dynamic elements, whereas algebraic equations describe that of time-independent elements. The equations also describe the relationships among those elements. Since a set of DAEs fully describes a system, it is also called a descriptor system.

The transient simulation method solves the DAEs of a subsystem in a straightforward way with respect to time. Each nonlinearity which exists in the subsystem is represented by a piecewise linear curve, since a nonlinearity is usually given as point-list data for most industrial systems. And the point-list data immediately give a piecewise linear representation. Since a subsystem in an industrial system is often stiff, it must be solved by an L-stable integration scheme which does not produce spurious or non-existing numerical oscillation. An L-stable integration scheme is an implicit scheme, and a set of simultaneous nonlinear equations therefore has to be solved at each time step. The piecewise linear representation of nonlinearity mentioned above is also suitable for solving the set of nonlinear equations. This is because its convergence in the solution process is detected by a fairly simple condition. After all, the transient simulation method is summarized as follows. First, all nonlinearities in the subsystem are represented by piecewise linear curves, and then the resultant set of DAEs is integrated using an L-stable integration scheme of your choice with respect to time by solving the set of piecewise linear equations obtained at each time step.

Considering the overall response time scale, the responses of some subsystems are quite fast and negligible and can be substituted by their steady-state ones. In this case, the steady-state response can be obtained by solving a set of algebraic equations describing the subsystem. The subsystem may operate on a sinusoidal carrier. If this is the case, the solution of the algebraic equations is complex-valued, where its absolute value and argument respectively give the magnitude and phase angle of the sinusoidal carrier. The algebraic equations are solved at each time step and repeated as time evolves, and the solution follows the variation of the sinusoidal carrier with respect to time. This method is often called the phasor simulation method. If the subsystem operates on a dc signal basis, then the solution of its algebraic equations becomes real-valued. This can be considered as a special case of the phasor simulation method.

omegasim uses a constant time step size for all subsystems, since this strategy is reasonable for the simulation of industrial systems. Consider this awkward but real situation. The designers of industrial systems are engineers. When those engineers implement, for instance, a digital PI controller, they use fairly simple numerical integration algorithm such as the forward Euler method, the trapezoidal method and so on with a constant time step. Do we have to use an elaborate numerical integration algorithm with variable time steps for the simulation of such a system, as suggested by applied mathematicians? The answer is of course not. This is the primary reason why I, an engineer, have started creating omegasim.

At each time step, the responses of all subsystems are calculated as their outputs. Then, considering the interconnected topology, those outputs are considered as the inputs to the subsystems with allowing one time step delay. The whole system is divided into subsystems so that the one time step delay does not affect the solution. For industrial systems, this process should not be difficult, since existing subsystems are actually designed to allow one time step delay. So, once the responses of all subsystems have been calculated, the output-to-input interconnection process can be carried out by simple substitutions, and this process is repeated with respect to time so that the response of the whole system is obtained.

An industrial system often has one or more electrical circuits as part of the system. To deal with electrical circuits, omegasim has methods to derive the DAEs of a given electrical circuit.

omegasim also has stability analysis functions ...

Author: Taku Noda
Started on Oct. 31, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

# -- Packages used in this modeule -------------------------------------

using EzXML
using LinearAlgebra
using SparseArrays
using Plots

# -- Structure definitions ---------------------------------------------

struct Component  #=
is used to store component definitions when reading an iSDF file.
=#
  name::String  # name of the component.
  type::String  # type of the component.
  calc::String  # content of the calc section of the component.
  root::EzXML.Node  # EzXML pointer to this component definition.
end

struct Subsystem  #=
is used to store subsystem definitions when reading an iSDF file.
=#
  name::String  # name of the subsystem.
  type::String  # type of the subsystem.
  calc::String  # content of the calc section of the subsystem.
  root::EzXML.Node  # EzXML pointer to this subsystem definition.
end

struct Use  #=
is used to store use declarations when reading an iSDF file.
=#
  cmpnnt_id::Int  # id of the component to be used.
  as::String  # name of the instance generated by this use declaration.
  use_elm::EzXML.Node  # EzXML pointer to this use declaration.
  subst::String  # content of the use declaration (variable substitutions).
end

struct Branch  #=
holds the information of a branch, or an electrical component. It is assumed that an array of the struct branch will be created for each circuit and the array holds information of all branches in the circuit.
=#
  # fields used by all branches:
  name::String  # name of the branch
  type::String  # type of the branch
  nph::Int  # number of phases
  val::Expr  # val(s) of the branch; given as an Expr to accept different formats defined by different types
  lhs::Vector{ Int }  # left-hand-side node numbers
  rhs::Vector{ Int }  # right-hand-side node numbers
  p_vars::Int  # pointer to the state variables of the branch
  # fields used by dynamic branches:
  p_de::Int  # pointer to the 1st differential equation
  # fields used by nonlinear and source branches:
  p_ns::Int  # pointer to nonlinear or source variable in u
  plx::Array{ Float64 }  # point list data (x, y) to represent the
  ply::Array{ Float64 }  # piecewise linear curve of this branch.
  op::Array{ Int }  # operating point of this branch.
end

struct Records  #=
is used to store calculation results .
=#
end

# -- Global variables --------------------------------------------------

t = 0.0  # time.

# -- Type definition ---------------------------------------------------

BranchTypes = (
  "R",  # linear resistor (single- or multiphase)
  "L",  # linear inductor (single- or multiphase)
  "G",  # linear conductance (single- or multiphase)
  "C",  # linear capacitor (single- or multiphase)
  "NR",  # nonlinear resistor (single-phase only)
  "NL",  # nonlinear inductor (single-phase only)
  "NC",  # nonlinear capacitor (single-phase only)
  "E",  # voltage source (single-phase only)
  "J",  # current source (single-phase only)
  "Sw",  # time controlled switch (single-phase only)
  "CR",  # controlled resistor (single-phase only)
  "CSw",  # controlled switch (single-phase only)
  "CE",  # controlled voltage source (single-phase only)
  "CJ",  # controlled current source (single-phase only)
  "VP",  # voltage probe (single-phase only)
  "CP" )  # current probe (single-phase only)

# -- Files to be included ----------------------------------------------

# iSDF-related source files:
include( "isdf_error.jl" )
include( "ordinal_string.jl" )
include( "check_component.jl" )
include( "read_isdf.jl" )
include( "expand_uses.jl" )
include( "isdf_pp1.jl" )
# omegasim source files:
include( "omegasim_error.jl" )
include( "collect_nodes.jl" )
include( "collect_branches.jl" )
include( "build_Ac.jl" )
include( "build_ZYB!.jl" )
include( "discretize_dyns!.jl" )
include( "update_y!.jl" )
# omegasim
include( "extract_nodes_branches.jl" )
include( "calc_submatrices.jl" )
include( "build_F0_tr.jl" )
include( "transient.jl" )

# -- Test code ---------------------------------------------------------

function test()
  filename = "test01.xml"
  h = 0.05e-3
  Tmax = 20.0e-3
  #
  nbset = extract_nodes_branches( filename )
  smset = calc_submatrices( nbset )
  F = build_F0_tr( nbset, smset )
  transient( F, nbset, h, Tmax )
  #
  return nothing
end

# -- This is the end of the module omegasim ----------------------------

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
module GraphLasso
using LinearAlgebra, SparseArrays
using DataFrames, Lasso, Graphs

export graphlasso

function graphlasso(S::Matrix{Float64}, α::Float64; tol::Float64=1e-5,
                    maxit::Int=1000, penalize_diag::Bool=true, verbose::Bool=false)
    p = size(S,1)
    adj = convert(SparseMatrixCSC{Int64,Int64}, abs.(S) .> α)
    blocks = connected_components(sparse2adjacencylist2(max.(adj, sparse(I, p, p))))
    if(verbose)
        print(blocks)
    end
    W = zeros(p,p)
    Θ = zeros(p,p)
    for block = blocks
        W[block,block], Θ[block,block] = fit_block(S[block,block], α, tol, maxit, penalize_diag)
    end
    return W, Θ
end

function fit_block(S::Matrix{Float64}, α::Float64, tol::Float64,
                   maxit::Int, penalize_diag::Bool)
    p = size(S,1)
    W = copy(S)
    if penalize_diag
        W += α*I
    end
    if size(S,1) == 1
        return W, inv(W)
    end
    Θ = zeros(p,p)
    i = 0
    β = zeros(p-1,p)
    while i < maxit
        i += 1
        W_old = copy(W)
        for j = 1:p
            inds = collect(1:p)
            splice!(inds, j)
            W11 = W[inds,inds]
            sqrtW11 = sqrt(W11)
            β[:,j] = fit(LassoPath, sqrtW11, sqrtW11 \ S[inds,j], λ=[α/(p-1)],
                         standardize=false, intercept=false).coefs
            W[inds,j] = W11 * β[:,j]
            W[j,inds] = W[inds,j]
        end
        if norm(W - W_old) < tol
            break
        end
    end
    if i == maxit
        warn("Maximum number of iterations reached, graphlasso failed to converge")
    end
    for j in 1:p
        inds = collect(1:p)
        splice!(inds, j)
        Θ[j,j] = 1/(W[j,j] - dot(W[inds,j],β[:,j]))
        Θ[inds,j] = -Θ[j,j] * β[:,j]
        Θ[j,inds] = Θ[inds,j]
    end
    return W, Θ
end

# Modify code from unmaintained Graph.jl
function sparse2adjacencylist2(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
    colptr = A.colptr
    rowval = A.rowval
    n = size(A, 1)
    adjlist = Array{Array{Ti,1}}(undef, n)
    s = 0
    for j in 1:n
        adjj = Ti[]
        sizehint!(adjj, colptr[j+1] - colptr[j] - 1)
        for k in colptr[j]:(colptr[j+1] - 1)
            rvk = A.rowval[k]
            if rvk != j push!(adjj, rvk) end
        end
        s += length(adjj)
        adjlist[j] = adjj
    end
    GenericAdjacencyList{Ti, UnitRange{Ti}, Vector{Vector{Ti}}}(!ishermitian(A),
                                                             one(Ti):convert(Ti,n),
                                                             s, adjlist)
end

end


text:
# NEEDED MODULES -- This can take time at the first launch
# Code main be needed to run twice to plot correctly (shift-enter)
using DifferentialEquations
using Plots
using LaTeXStrings
# Load the backend (engine) for plotting based on javascript
plotlyjs()

function weighted_leastsquare(O,R,E)
    denom  = sum(O.^2)
    weights = 1. / (E.^2)
    return sqrt(sum(weights.*(O-R).^2)/denom)
end

function leastsquare(O,R)
    denom  = sum(O.^2)
    return sqrt(sum((O-R).^2)/denom)
end

# Saturation function to positive values
# (concentrations are positive function)
# This is just used as a safety for the parameter search phase
@inline function sat(x)
    if x < 0.0
        return 0.0
    else
        x
    end
end

# Gaussian model of plasma Glucose appearance rate
function GlcGauss(t)
    σ = 35.0 # width
    μ = 30.0 # peak time
    dose = 3729 # parameter K in paper
    return dose*1/(σ*sqrt(2*pi))*exp((-(t-μ)^2/(2*σ^2)))
end

# Parametric ODEs of te Glucose response model
# t : time
# u : variables
# p : parameters
# Here I call the guassian model of glucose input GlcGauss(t)
# And I saturate the variables to the positive cone
function f_param!(du,u,p,t)
 du[1] = -p[1]*(sat(u[1])-p[2])- p[10]*sat(u[2])*sat(u[1])+GlcGauss(t)
 du[2] = -p[3]*sat(u[2])+p[4]*(sat(u[3])-p[5])
 du[3] = -p[6]*sat(u[3])+p[7]*(sat(u[1])-p[8])*t + p[9]*u[4]
 du[4] = -p[9]*u[4]
end

# Initial conditions at 21j
u0 = [110.0,0.0,16.0,5950]

tspan = (0.0,120.0)

# #  A Parameter set from the MINMOD model as a reference
# #       p1  Gb    p2    p3     Ib  n   γ    h     |added parameters => |     p4    r
# p0 = [0.01,100.0,0.01,0.00001,8.0,0.1,0.01,100.0,                          0.07, 1.0]

# Fitted parameter set to the control group at 21 Days ===> err: 0.00311
#         p1    Gb    p2     p3     Ib     n      γ      h     p4     r
  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0] # err = 0.00311
###################################################
## Simulation of the Control experiment at 21Days
prob_param = ODEProblem(f_param!,u0,tspan,p0)
sctrl = solve(prob_param, alg_hints = [:stiff])

## Experimental results
texp = [0,10,20,30,45,60,90,120]
glcexp = [110.1,178.2,209.3,203.6,187.4,161.9,131.4,114.5]
uppererror = [1.4,4.0,3.3,2.6,2.8,2.4,1.9,1.5]
lowererror = [1.4,4.0,3.3,2.6,2.8,2.4,1.9,1.5]
errs = [lowererror;uppererror]
res_tot = [sctrl(ti) for ti in texp]
res = [ri[1] for ri in res_tot]

## Computation of the weighted least square error and least square
wls_err_Ctrl = weighted_leastsquare(glcexp,res,uppererror)
print("CONTROL wls = ", wls_err_Ctrl, "\n")
ls_err_Ctrl = leastsquare(glcexp,res)

## Generation of the figure and first plots
plot(sctrl,vars=[(0,1)],label="Glucose (control) simulations")
scatter!(texp,glcexp,yerr = errs,marker = stroke(2, :blue),
       xlabel="time (min)", ylabel= "Concentration: mg/dL",
    label="Control exp results", legend=:top, markerstrokecolor = :auto)

# Parameters sets for the different hypothesis for the group CD1 at 21 Days
########################################################################
# Control parameter set  ===>   err = 0.0188
#         p1    Gb    p2     p3     Ib     n      γ      h     p4     r
  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0]
# hypothese 1.1: sensibility to interstitial insuline (r) ===> err : 0.0144
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 0.88]
# hypothese 1.2: disponibility of insuline for X(t) (p3) ==> err: 0.0144
#  p0 = [0.01, 100.0, 0.56, 0.0136, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0]
# hypothese 2: insuline degradation rate (n) ==> err : 0.0145
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 11.10, 0.0310, 85.0, 0.033, 1.0]
# hypothese 3: interstitial insuline degradation (p2) ==> err : 0.0140
# p0 = [0.01, 100.0, 0.65, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0]
# hypothesis 4: Slow phase 'acceleration' γ decreased  ==> err = 0.00566
 p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0258, 85.0, 0.033, 1.0]
# hypothesis 5:  change p4 (stored insuline release rate)  ==> err = 0.0188
#   p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0]
# hypothese 6: reponse Glucose threshold (h)  ==> err: 0.0102
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 100.0, 0.033, 1.0]
##############################################################################
# Simulation CD1 21 Days
prob_param = ODEProblem(f_param!,u0,tspan,p0)
scd1 = solve(prob_param, alg_hints = [:stiff])
# Experimental results
glcexpCd1 = [117.5,174.4,222.0,220.6,204.2,177.1,135.1,119.7]
uppererror = [1.8,2.7,3.8,4.1,3.7,2.9,2.2,1.5]
lowererror = [1.8,2.7,3.8,4.1,3.7,2.9,2.2,1.5]
errs = [lowererror;uppererror]
res_tot = [scd1(ti) for ti in texp]
res = [ri[1] for ri in res_tot]
# Error computation
wls_err_Cd1 = weighted_leastsquare(glcexpCd1,res,uppererror)
print("CD1 wls = ", wls_err_Cd1, "\n")
ls_err_Cd1 = leastsquare(glcexpCd1,res)
# Plotting
plot!(scd1,vars=[(0,1)],label="Glucose(Cd1) hyp 3 simu")
scatter!(texp,glcexpCd1,yerr = errs,marker = stroke(2, :green),
      xlabel="time (min)", ylabel= "Concentration: mg/dL",
      label="Cd1 exp results", legend=:top, markerstrokecolor = :auto)

# Parameters sets for the different hypothesis for the group CD2 at 21 Days
########################################################################
# Control parameter set  ===> 0.0144
#         p1    Gb    p2     p3     Ib     n      γ      h     p4     r
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0]
# hypothese 1.1: sensibility to interstitial insuline (r) ==> err: 0.0082
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 0.70]
# hypothese 1.2: disponibility of insuline for X(t) (p3) ==> err : 0.00821
#  p0 = [0.01, 100.0, 0.56, 0.0110, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0]
# hypothese 2: degradation  insuline (n) ==> err : 0.00852
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 12.05, 0.0310, 85.0, 0.033, 1.0]
# hypothese 3: interstitial insuline degradation (p2) ==> err: 0.00789
#  p0 = [0.01, 100.0, 0.81, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.033, 1.0]
# hypothesis 4: Slow phase 'acceleration' γ decreased ==> err =  0.00433
 p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0215, 85.0, 0.033, 1.0]
# hypothesis 5:  change p4 (stored insuline release rate)  ==> err = 0.01370
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 85.0, 0.029, 1.0]
# hypothese 6: reponse Glucose threshold (h)  ==> err: 0.00740
#  p0 = [0.01, 100.0, 0.56, 0.0155, 10.0, 10.53, 0.0310, 115.0, 0.033, 1.0]
##############################################################################
# Simulation CD2 21 Days
prob_param = ODEProblem(f_param!,u0,tspan,p0)
scd2 = solve(prob_param, alg_hints = [:stiff])
# Experimental results
glcexpCd2 = [105.0,196.9,226.8,243.7,220.7,191.7,151.6,119.2]
uppererror = [2.9,10.2,10.1,9.6,7.7,5.4,5.0,3.9]
lowererror = [2.9,10.2,10.1,9.6,7.7,5.4,5.0,3.9]
errs = [lowererror;uppererror]
res_tot = [scd2(ti) for ti in texp]
res = [ri[1] for ri in res_tot]
# Error computation
wls_err_Cd2 = weighted_leastsquare(glcexpCd2,res,uppererror)
print("CD2 wls = ", wls_err_Cd2, "\n")
ls_err_Cd2 = leastsquare(glcexpCd2,res)
# Plotting
plot!(scd2,vars=[(0,1)],label="Glucose (Cd2) simulations")
cur_colors = get_color_palette(:auto, plot_color(:white), 17)
scatter!(texp,glcexpCd2,yerr = errs,marker = stroke(2, cur_colors[5]),
       xlabel="time (min)", ylabel= "Concentration: mg/dL",
       label="Cd2 exp results", legend=:top, markerstrokecolor = :auto)

# Additional Plot options
plot!(size=(1280,960),legendfontsize=15, xtickfont=font(13), ytickfont=font(13), guidefont=font(18))


text:
module TypeTraits

export @trait, @mixin

# global registry of all abstract classes defined
traits_declarations = Dict{Symbol, Array}()

# @trait type X
# ...
# end
macro trait(typedef)
  # TODO: add checks that no inheritance occurs here

  sym = typedef.args[2]
  block = typedef.args[3]

  traits_declarations[sym] = [Expr(:(::), var.args[1], var.args[2]) for var in block.args[2:2:end]]
end

# @mixin type Foo(X) <: Y end
macro mixin(typedef)
  decl = typedef.args[2]
  block = typedef.args[3]

  if decl.head === :call
    sym = decl.args[1]
    mixin = decl.args[2]
    parent = None
  elseif decl.head === :<:
    traits = decl.args[1]
    sym = traits.args[1]
    mixin = traits.args[2]
    parent = decl.args[2]
  end

  trait_declaration = traits_declarations[mixin]
  child_declarations = [Expr(:(::), var.args[1], var.args[2]) for var in block.args[2:2:end]]
  declarations = vcat(trait_declaration, child_declarations)

  if parent !== None
    Expr(:type, true, Expr(:<:, esc(sym), esc(parent)), Expr(:block, declarations...))
  else
    Expr(:type, true, esc(sym), Expr(:block, declarations...))
  end
end

end # module


text:
using aoc20
using Test, aoc20

@testset "aoc20.jl" begin
    @test_nowarn aoc20.benchmarkAll(onlyOnce = true)
    @test aoc20.formatTime(0.1) == "100.000 ms"
end

@testset "day1" begin
    @testset "part1" begin
        @test part1(day = 1) == 1006875
    end

    @testset "part2" begin
        @test part2(day = 1) == 165026160
    end
end

@testset "day2" begin
    @testset "part1" begin
        @test aoc20.day2.part1(["1-3 a: abaca"]) == 1
        @test aoc20.day2.part1(["1-3 a: abaaca"]) == 0
        @test aoc20.day2.part1(["1-3 a: bc"]) == 0
        @test aoc20.day2.part1(["0-0 a: abaaca"]) == 0
        @test part1(day = 2) == 506
    end

    @testset "part2" begin
        @test aoc20.day2.part2(["1-3 a: abc"]) == 1
        @test aoc20.day2.part2(["1-3 a: cba"]) == 1
        @test aoc20.day2.part2(["1-3 a: aba"]) == 0
        @test aoc20.day2.part2(["1-3 a: bac"]) == 0
        @test part2(day = 2) == 443
    end
end

@testset "day3" begin
    @testset "part1" begin
        x = [ 0 0 0 0 0 ;
              0 0 0 1 0 ;
              0 1 0 0 0 ;
              0 0 0 0 1]'
        @test aoc20.day3.part1(x) == 3
        @test part1(day = 3) == 257
    end

    @testset "part2" begin
        @test part2(day = 3) == 1744787392
    end
end

@testset "day3" begin
    @testset "part1" begin
        x = [ 0 0 0 0 0 ;
              0 0 0 1 0 ;
              0 1 0 0 0 ;
              0 0 0 0 1]'
        @test aoc20.day3.part1(x) == 3
        @test part1(day = 3) == 257
    end

    @testset "part2" begin
        @test part2(day = 3) == 1744787392
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using SpecialFunctions
using Random123
import AdaptiveRejectionSampling

"""
    struct BAT.GaussianShell <: Distribution{Multivariate,Continuous}

*Experimental feature, not part of stable public API.*

Gaussian Shell (see
[Caldwell et al.](https://arxiv.org/abs/1808.08051) for definition).

Constructors:

* ```BAT.GaussianShell(; r::Real=5, w::Real=2, n::Integer=2)```

Fields:

$(TYPEDFIELDS)

!!! note

    Fields `c` and `lognorm` are considered internal and subject to
    change without deprecation.
"""
struct GaussianShell{T<:Real, V<:AbstractVector{<:Real}, F<:AbstractFloat} <: Distribution{Multivariate,Continuous}
    "The radius of the Gaussian shell distribution."
    r::T

    "Variance of the Gaussian shell distribution"
    w::T

    "Number of dimensions"
    n::Int

    c::V

    lognorm::F
end

function GaussianShell(;r::Real=5, w::Real=2, n::Integer=2)
    c = zeros(n)
    r,w = promote(r, w)
    radial_integral(ρ) = ρ^(n-1)*exp(-(ρ - r)^2 / (2*w^2))
    density_norm = sqrt(2*π*w^2)
    # Normalization for coordinate transform
    radial_norm = (sqrt(2)*π^((n-1)/2)) / (gamma(n/2)*w)*QuadGK.quadgk(radial_integral, 0, r+w*20)[1]
    lognorm = log(density_norm) + log(radial_norm)
    GaussianShell(r, w, n, c, lognorm)
end

function _find_gs_search_range(ndims::Integer, r0::Real, w::Real)
    f(r) = gs_pdf_r(r, ndims, r0, w)
    o = AdaptiveRejectionSampling.Objective(f)

    r_grid = Array(range(maximum([r0/2, r0-w]), stop=sqrt(2*ndims)*w+r0/(w^2) + sqrt(ndims)*r0/w, length=10^6))
    y_grid = f.(r_grid)
    y_grads = o.grad.(r_grid)

    # Both sides of a concave peak
    idx1 = argmax(y_grads)
    idx2 = argmin(y_grads)

    return (r_grid[idx1], r_grid[idx2]+sqrt(ndims)*r0)
end

nball_surf_area(r, ndims) = 2 * π^(ndims/2) / gamma(ndims/2) * r^(ndims-1)
gs_pdf_r(r, ndims, r0, w) = nball_surf_area(r, ndims) * 1/(2π * w^2) * exp(-(abs(r) - r0)^2 / (2 * w^2))

function gauss_shell_radial_samples(rng::AbstractRNG, ndims::Integer, r0::Real, w::Real, n::Integer)
    f(r) = gs_pdf_r(r, ndims, r0, w)
    result = Vector{Float64}()

    # ToDo: AdaptiveRejectionSampling sometimes fails sporadically here,
    # find a fix to replace this workaround:
    ntries::Int = 0
    while length(result) < n
        ntries += 1
        try
            sampler = AdaptiveRejectionSampling.RejectionSampler(f, (10^-10, Inf), search_range=_find_gs_search_range(ndims, r0, w))
            smpls = AdaptiveRejectionSampling.run_sampler!(rng, sampler, n)
            append!(result, smpls)
        catch err
            if !(err isa AssertionError || err isa DomainError) || ntries > 4
                rethrow()
            end
        end
    end

    result
end

function rand_nball_surf_samples!(rng::AbstractRNG, X::AbstractMatrix)
    randn!(rng, X)
    X ./= sqrt.(sum(X .* X, dims = 1))
    return X
end

Distributions.rand(rng::AbstractRNG, d::GaussianShell) = vec(rand(rng, d, 1))

function Distributions.rand(rng::AbstractRNG, d::GaussianShell, n::Int)
    X = Matrix{eltype(d)}(undef, length(d), n)
    Distributions._rand!(rng, d, X)
end

function Distributions._rand!(rng::AbstractRNG, d::GaussianShell, X::AbstractMatrix)
    rand_nball_surf_samples!(rng, X)
    R = gauss_shell_radial_samples(rng, d.n, d.r, d.w, size(X,2))
    X .*= R'
    return X
end

function Distributions._logpdf(d::GaussianShell, x::AbstractArray)
    integral_result = -(sqrt(sum((x .- d.c).^2)) - d.r)^2 / (2*d.w^2)
    result = integral_result - d.lognorm
    return result
end

function Statistics.cov(d::GaussianShell)
    cov(nestedview(rand(bat_determ_rng(), bat_sampler(d), 10^5)))
end

Base.length(d::GaussianShell) = length(d.c)

Base.eltype(d::GaussianShell) = eltype(d.c)

Distributions.params(d::GaussianShell) = (d.r, d.w, d.c)


text:
__precompile__(true)
module XLSX

import ZipFile, EzXML, Missings

# https://github.com/fhs/ZipFile.jl/issues/39
if !method_exists(Base.nb_available, Tuple{ZipFile.ReadableFile})
    Base.nb_available(f::ZipFile.ReadableFile) = f.uncompressedsize - f._pos
end

const SPREADSHEET_NAMESPACE_XPATH_ARG = [ "xpath" => "http://schemas.openxmlformats.org/spreadsheetml/2006/main" ]

include("structs.jl")
include("cellref.jl")
include("sst.jl")
include("stream.jl")
include("table.jl")
include("relationship.jl")
include("read.jl")
include("workbook.jl")
include("worksheet.jl")
include("cell.jl")
include("styles.jl")
include("helper.jl")
include("deprecated.jl")
include("write.jl")

end # module XLSX


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

if Sys.iswindows()
    ttyhascolor(term_type = nothing) = true
else
    function ttyhascolor(term_type = get(ENV, "TERM", ""))
        startswith(term_type, "xterm") && return true
        try
            @static if Sys.KERNEL === :FreeBSD
                return success(`tput AF 0`)
            else
                return success(`tput setaf 0`)
            end
        catch e
            return false
        end
    end
end
function get_have_color()
    global have_color
    have_color === nothing && (have_color = ttyhascolor())
    return have_color::Bool
end
in(key_value::Pair{Symbol,Bool}, ::TTY) = key_value.first === :color && key_value.second === get_have_color()
haskey(::TTY, key::Symbol) = key === :color
getindex(::TTY, key::Symbol) = key === :color ? get_have_color() : throw(KeyError(key))
get(::TTY, key::Symbol, default) = key === :color ? get_have_color() : default

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
module MiniMax

export play_bot_turn!, State, MARK_P1, MARK_P2, MARK_TIE, MARK_EMPTY

using Random

const MARK_P1 = Int8(1)
const MARK_P2 = Int8(-1)
const MARK_EMPTY = Int8(0)
const MARK_TIE = Int8(2)

const ONE = Int8(1)

const IdxTable = Int8[1, 2, 3, 4, 5, 6, 7, 8, 9]

mutable struct State
    player::Int8
    winner::Int8
    last_entry::Int8
    board::Array{Int8,2}
end

function play_bot_turn!(board::Array{Int8,2}, game_type::String = "normal")::State
    state = State(MARK_P1, MARK_EMPTY, MARK_EMPTY, copy(board))

    find_optimal_move!(state, game_type)

    state
end

function find_optimal_move!(state::State, game_type::String)
    init_depth::Int8 = sum(==(MARK_EMPTY), state.board)

    if init_depth == length(state.board)
        mid_cell::Int8 = ceil(length(state.board) / 2)

        state.board[mid_cell] = state.player
        state.last_entry = mid_cell
        return
    end

    best_move::Int8 = MARK_EMPTY

    if game_type == "normal"
        α::Int8 = typemin(Int8) # min value that the maximizing player get certainly
        β::Int8 = typemax(Int8) # max value that the minimizing player get certainly

        _, best_move = minimax!(state, init_depth, α, β, true)
    else
        best_move = novice_move(state)
    end

    state.last_entry = best_move

    if state.last_entry > MARK_EMPTY
        state.board[state.last_entry] = state.player
    end

    win_state::Int8 = compute_win_state(state)

    if win_state == 0 && sum(==(MARK_EMPTY), state.board) == 0
        state.winner = MARK_TIE
    else
        state.winner = win_state
    end

end

function minimax!(
    state::State,
    depth::Int8,
    α::Int8,
    β::Int8,
    maximize::Bool,
)::Tuple{Int8,Int8}
    win_state::Int8 = compute_win_state(state)

    if win_state == MARK_P1 || win_state == MARK_P2 || depth == 0
        return win_state * (depth + 1), Int8(0)
    end

    best_entry::Int8 = 0
    best_value::Int8 = (maximize == true) ? typemin(Int8) : typemax(Int8)
    current_player::Int8 = state.player

    allowed_entries::Array{Int8,1} = allowed_board_entries(state)

    for entry in allowed_entries
        state.board[entry] = current_player
        state.player = next_player(current_player)

        if maximize == true
            value, _ = minimax!(state, depth - ONE, α, β, false)

            if value > best_value
                best_value = value
                best_entry = entry
            end

            state.board[entry] = MARK_EMPTY
            state.player = current_player

            if β <= α # prune rest branches
                break
            end
            α = max(α, best_value)
        else
            value, _ = minimax!(state, depth - ONE, α, β, true)

            if value < best_value
                best_value = value
                best_entry = entry
            end

            state.board[entry] = MARK_EMPTY
            state.player = current_player

            if β <= α # prune rest branches
                break
            end
            β = min(β, best_value)
        end
    end

    best_value, best_entry
end

function allowed_board_entries(state::State)::Array{Int8,1}
    entries::Array{Int8,1} = []

    for cell::Int8 = 1:length(state.board)
        if state.board[cell] == MARK_EMPTY
            push!(entries, cell)
        end
    end

    shuffle(entries)
end

function next_player(player::Int8)::Int8
    (player == MARK_P1) ? MARK_P2 : MARK_P1
end

function compute_win_state(state::State)::Int8
    rows = size(state.board)[begin]
    cols = size(state.board)[end]

    row_states = [
        all(state.board[r, begin] .== state.board[r, :]) ? state.board[r, begin] :
        MARK_EMPTY for r = 1:rows
    ]
    col_states = [
        all(state.board[begin, c] .== state.board[:, c]) ? state.board[begin, c] :
        MARK_EMPTY for c = 1:cols
    ]

    diag_state =
        (
            state.board[begin, begin] == state.board[2, 2] &&
            state.board[2, 2] == state.board[end, end]
        ) ? state.board[begin, begin] : MARK_EMPTY
    antidiag_state =
        (
            state.board[begin, end] == state.board[2, 2] &&
            state.board[2, 2] == state.board[end, begin]
        ) ? state.board[begin, end] : MARK_EMPTY

    max_state::Int8 = max(row_states..., col_states..., diag_state, antidiag_state)
    min_state::Int8 = min(row_states..., col_states..., diag_state, antidiag_state)

    (max_state > MARK_EMPTY) ? max_state : min_state
end

function novice_move(state::State)::Int8
    win_state::Int8 = compute_win_state(state)

    if win_state == MARK_P1 || win_state == MARK_P2
        return Int8(0)
    end

    free_indices::Array{Int8,1} =
        filter(x -> x > 0, reshape(state.board .== 0, (length(state.board),)) .* IdxTable)

    (length(free_indices) == 0) ? Int8(0) : shuffle(free_indices)[begin]
end

end

max_stars_repo_name: elmomoilanen/Tic-tac-toe

max_issues_repo_name: elmomoilanen/Tic-tac-toe

max_forks_repo_name: elmomoilanen/Tic-tac-toe


text:
"""
    precompile_pather(packageName::String)

To get the path of precompile_packageName.jl file

Written exclusively for SnoopCompile Github actions.
# Examples
```julia
precompilePath, precompileFolder = precompile_pather("MatLang")
```
"""
function precompile_pather(packageName::String)
    return "\"../deps/SnoopCompile/precompile/precompile_$packageName.jl\"",
    "$(pwd())/deps/SnoopCompile/precompile/"
end

precompile_pather(packageName::Symbol) = precompile_pather(string(packageName))
precompile_pather(packageName::Module) = precompile_pather(string(packageName))

################################################################

function precompile_regex(precompilePath)
    # https://stackoverflow.com/questions/3469080/match-whitespace-but-not-newlines
    # {1,} for any number of spaces
    c1 = Regex("#[^\\S\\r\\n]{0,}include\\($(precompilePath)\\)")
    c2 = r"#\s{0,}_precompile_\(\)"
    a1 = "include($precompilePath)"
    a2 = "_precompile_()"
    return c1, c2, a1, a2
end
################################################################

"""
    precompile_activator(packagePath, precompilePath)

Activates precompile of a package by adding or uncommenting include() of *.jl file generated by SnoopCompile and _precompile_().

packagePath is the same as `pathof`. However, `pathof(module)` isn't used to prevent loadnig the package.

Written exclusively for SnoopCompile Github actions.
"""
function precompile_activator(packagePath::String, precompilePath::String)

    packageText = Base.read(packagePath, String)

    c1, c2, a1, a2 = precompile_regex(precompilePath)

    # Checking availability of _precompile_ code
    commented = occursin(c1, packageText) && occursin(c2, packageText)
    available = occursin(a1, packageText) && occursin(a2, packageText)

    if commented
        packageEdited = foldl(replace,
                     (
                      c1 => a1,
                      c2 => a2,
                     ),
                     init = packageText)

                     Base.write(packagePath, packageEdited)

        println("precompile is activated")
    elseif available
        # do nothing
        println("precompile is already activated")
    else
        # TODO: add code automatiaclly
        error(""" add the following codes into your PackageName.jl file under src folder:
         #include($precompilePath)
         #_precompile_()
         """)
    end

end

"""
    precompile_deactivator(packagePath, precompilePath)

Deactivates precompile of a package by commenting include() of *.jl file generated by SnoopCompile and _precompile_().

packagePath is the same as `pathof`. However, `pathof(module)` isn't used to prevent loadnig the package.

Written exclusively for SnoopCompile Github actions.
"""
function precompile_deactivator(packagePath::String, precompilePath::String)

    packageText = Base.read(packagePath, String)

    c1, c2, a1, a2 = precompile_regex(precompilePath)

    # Checking availability of _precompile_ code
    commented = occursin(c1, packageText) && occursin(c2, packageText)
    available = occursin(a1, packageText) && occursin(a2, packageText)

    if available && !commented
        packageEdited = foldl(replace,
                     (
                      a1 => "#"*a1,
                      a2 => "#"*a2,
                     ),
                     init = packageText)

                     Base.write(packagePath, packageEdited)

        println("precompile is deactivated")
    elseif commented
        # do nothing
        println("precompile is already deactivated")
    else
        # TODO: add code automatiaclly
        error(""" add the following codes into your PackageName.jl file under src folder:
         #include($precompilePath)
         #_precompile_()
         """)
    end

end


text:
using CompScienceMeshes
using BEAST
using Test

p1 = point(0,0,0)
p2 = point(1/3,0,0)
p3 = point(0,1/3,0)

q1 = point(0,0,0)
q2 = point(0.125,0,0)

m = Mesh([p1,p2,p3], [index(1,2,3)])
n = Mesh([q1,q2], [index(1,2)])
translate!(n, point(0,0,20))

X = lagrangec0d1(m, boundary(m))
@test numfunctions(X) == 3

x = refspace(X)
s = chart(m, m.faces[X.fns[1][1].cellid])
c = neighborhood(s, [1,1]/3)      # get the barycenter of that patch
v = x(c, Val{:withcurl})       # evaluate the Lagrange elements in c, together with their curls

@test (s[3] - s[2]) / (2 * volume(s)) ≈ v[1][2]
@test (s[1] - s[3]) / (2 * volume(s)) ≈ v[2][2]
@test (s[2] - s[1]) / (2 * volume(s)) ≈ v[3][2]

Y = lagrangec0d1(n, boundary(n))
@test numfunctions(Y) == 2

y = refspace(Y)
t = chart(n, n.faces[Y.fns[1][1].cellid])
d = neighborhood(t, [1]/2)
tg = normalize(tangents(d,1))
w = y(d)

@test w[1][1] ≈ 0.5
@test w[2][1] ≈ 0.5

κ = 0.0
T = BEAST.NitscheHH3(κ)
Tyx = assemble(T, Y, X)

@test size(Tyx) == (numfunctions(Y), numfunctions(X))

R = norm(cartesian(c)-cartesian(d))
estimate = volume(t) * volume(s) * dot(v[1][2], w[1][1] * tg) / (4π*R)
actual = Tyx[1,1]

@test (estimate - actual) / actual < 0.005

## test the value of the skeleton gram matrix
I = BEAST.Identity()
Iyy = assemble(I, Y, Y)

@test size(Iyy) == (2,2)

qps = BEAST.quadpoints(y, [t], (10,))[1,1]
estimate = 0.0
for qp in qps
    igd = qp.value[1][1] * qp.value[2][1]
    global estimate += qp.weight * igd
end
actual = Iyy[1,2]
@test norm(estimate - actual) < 1e-6


text:
export is_zero_row, modular_hnf, submat, howell_form, _hnf_modular, kernel_mod, matrix, zeromatrix

#

function matrix{T <: RingElem}(A::Array{T, 2})
  r, c = size(A)
  (r < 0 || c < 0) && error("Array must be non-empty")
  m = MatrixSpace(parent(A[1, 1]), size(A)...)(A)
  return m
end

function matrix{T}(R::Ring, n::Int, m::Int, A::Array{T, 2})
  m = MatrixSpace(R, n, m)(A)
  return m
end

function zeromatrix(R::Ring, n::Int, m::Int)
  return MatrixSpace(R, n, m)()
end

#

function Array{T}(a::fmpz_mat; S::Type{T} = fmpz)
  A = Array(T, rows(a), cols(a))
  for i = 1:rows(a)
    for j = 1:cols(a)
      A[i,j] = T(a[i,j])
    end
  end
  return A
end

function is_zero_row(M::fmpz_mat, i::Int)
  for j = 1:cols(M)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function is_zero_row(M::nmod_mat, i::Int)
  zero = UInt(0)
  for j in 1:cols(M)
    t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &M, i - 1, j - 1)
    if t != zero
      return false
    end
  end
  return true
end

function is_zero_row{T}(M::MatElem{T}, i::Int)
  for j in 1:cols(M)
    if !iszero(M[i,j])
      return false
    end
  end
  return true
end

function is_zero_row{T <: Integer}(M::Array{T, 2}, i::Int)
  for j = 1:Base.size(M, 2)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function is_zero_row(M::Array{fmpz, 2}, i::Int)
  for j = 1:Base.size(M, 2)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function is_zero_row{T <: RingElem}(M::Array{T, 2}, i::Int)
  for j in 1:Base.size(M, 2)
    if !iszero(M[i,j])
      return false
    end
  end
  return true
end

function divexact!(a::fmpz_mat, b::fmpz_mat, d::fmpz)
  ccall((:fmpz_mat_scalar_divexact_fmpz, :libflint), Void,
               (Ptr{fmpz_mat}, Ptr{fmpz_mat}, Ptr{fmpz}), &a, &a, &d)
end

function mul!(a::fmpz_mat, b::fmpz_mat, c::fmpz)
  ccall((:fmpz_mat_scalar_mul_fmpz, :libflint), Void,
                  (Ptr{fmpz_mat}, Ptr{fmpz_mat}, Ptr{fmpz}), &a, &b, &c)
end

#computes (hopefully) the hnf for vcat(a*I, m) and returns ONLY the
#non-singular part. By definition, the result wil have full rank
#
#Should be rewritten to use Howell and lifting rather the big HNF
#
function modular_hnf(m::fmpz, a::fmpz_mat, shape::Symbol = :upperright)
  c = vcat(parent(a)(m), a)
  n = cols(a)
  w = window(c, n+1, 1, 2*n, n)
  ccall((:fmpz_mat_scalar_mod_fmpz, :libflint), Void, (Ptr{fmpz_mat}, Ptr{fmpz_mat}, Ptr{fmpz}), &w, &w, &m)
  if shape == :lowerleft
    c = _hnf(c, shape)
    return sub(c, n+1:2*n, 1:n)
  else
    c = hnf(c)
    c = sub(c, 1:n, 1:n)
  end
end

#function _lift_howell_to_hnf(x::nmod_mat)
## Assume that x is square, in howell normal form and all non-zero rows are at the bottom
## NOTE: _OUR_ Howell normal form algorithm always puts the rows at the right position
## If row i is non-zero then i is the rightmost non-zero entry
## Thus lifting is just replacing zero diagonal entries
#  !issquare(x) && error("Matrix has to be square")
#  y = lift_unsigned(x)
#  for i in cols(y):-1:1
#    z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &y, i - 1, i - 1)
#    if Bool(ccall((:fmpz_is_zero, :libflint), Int, (Ptr{fmpz}, ), z))
##      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &y, 0, i - 1)
#      ccall((:fmpz_set_ui, :libflint), Void, (Ptr{fmpz}, UInt), z, x._n)
##      for k in 1:i-1
##        _swaprows!(y, k, k+1)
##      end
#    end
#  end
#  return y
#end

function submat{T <: Integer}(x::nmod_mat, r::UnitRange{T}, c::UnitRange{T})
  z = deepcopy(window(x, r, c))
  return z
end

function submat{T <: Integer}(x::fmpz_mat, r::UnitRange{T}, c::UnitRange{T})
  z = deepcopy(window(x, r, c))
  return z
end

#function _hnf_modular(x::fmpz_mat, m::fmpz, shape::Symbol = :lowerleft)
#  if abs(m) < fmpz(typemax(UInt))
#    y = MatrixSpace(ResidueRing(FlintZZ, m), rows(x), cols(x))(x)
#    howell_form!(y, shape)
#    y = submat(y, rows(y) - cols(y) + 1:rows(y), 1:cols(y))
#    return _lift_howell_to_hnf(y)
#  end
#  return __hnf_modular(x, m, shape)
#end
#
#function __hnf_modular(x::fmpz_mat, m::fmpz, shape::Symbol = :lowerleft)
## See remarks above
#  y = deepcopy(x)
#  howell_form!(y, m, shape)
#  y = submat(y, rows(y) - cols(y) + 1:rows(y), 1:cols(y))
#  for i in cols(y):-1:1
#    z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &y, i - 1, i - 1)
#    if Bool(ccall((:fmpz_is_zero, :libflint), Int, (Ptr{fmpz}, ), z))
#    #if ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, i - 1, i - 1) == 0
##      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &y, 0, i - 1)
#      ccall((:fmpz_set, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), z, &m)
##      for k in 1:i-1
##        _swaprows!(y, k, k+1)
##      end
#    end
#  end
#  return y
#end

function _hnf(x::fmpz_mat, shape::Symbol = :upperright)
  if shape == :lowerleft
    h = hnf(_swapcols(x))
    _swapcols!(h)
    _swaprows!(h)
    return h::fmpz_mat
  end
  return hnf(x)::fmpz_mat
end

function _hnf_modular_eldiv(x::fmpz_mat, m::fmpz, shape::Symbol = :upperright)
  if shape == :lowerleft
    h = hnf_modular_eldiv(_swapcols(x), m)
    _swapcols!(h)
    _swaprows!(h)
    return h
  elseif shape == :upperright
    return hnf_modular_eldiv(x, m)
  end
end

function hnf_modular_eldiv!(x::fmpz_mat, d::fmpz)
   (rows(x) < cols(x)) &&
                error("Matrix must have at least as many rows as columns")
   ccall((:fmpz_mat_hnf_modular_eldiv, :libflint), Void,
                (Ptr{fmpz_mat}, Ptr{fmpz}), &x, &d)
   return x
end

#function howell_form!(x::fmpz_mat, m::fmpz, shape::Symbol = :upperright)
#  if shape == :lowerleft
#    _swapcols!(x)
#    ccall((:_fmpz_mat_howell, :libflint), Int, (Ptr{fmpz_mat}, Ptr{fmpz}), &x, &m)
#    _swapcols!(x)
#    _swaprows!(x)
#  else
#    ccall((:_fmpz_mat_howell, :libflint), Int, (Ptr{fmpz_mat}, Ptr{fmpz}), &x, &m)
#  end
#end
#
#function howell_form(x::fmpz_mat, m::fmpz, shape::Symbol = :upperright)
#  y = deepcopy(x)
#  howell_form!(y, m, shape)
#  return y
#end
#
#function howell_form!(x::nmod_mat, shape::Symbol = :upperright)
#  if shape == :lowerleft
#    _swapcols!(x)
#    ccall((:_nmod_mat_howell, :libflint), Int, (Ptr{nmod_mat}, ), &x)
#    _swapcols!(x)
#    _swaprows!(x)
#  else
#    ccall((:_nmod_mat_howell, :libflint), Int, (Ptr{nmod_mat}, ), &x)
#  end
#end
#
#function howell_form(x::nmod_mat, shape::Symbol = :upperright)
#  y = deepcopy(x)
#  howell_form!(y, shape)
#  return y
#end

function _swaprows(x::fmpz_mat)
  y = deepcopy(x)
  _swaprows!(y)
  return y
end

function _swapcols(x::fmpz_mat)
  y = deepcopy(x)
  _swapcols!(y)
  return y
end

function _swaprows!(x::fmpz_mat)
  r = rows(x)
  c = cols(x)

  if r == 1
    return x
  end

  if r % 2 == 0
    for i in 1:div(r,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, i - 1, j - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, (r - i + 1) - 1, j - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  else
    for i in 1:div(r-1,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, i - 1, j - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, (r - i + 1) - 1, j - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  end
  nothing
end

function _swaprows!(x::fmpz_mat, i::Int, j::Int)
  ccall((:_fmpz_mat_swap_rows, :libflint), Void, (Ptr{fmpz_mat}, Int, Int), &x, i-1, j-1)
  nothing
end

function _swaprows!(x::nmod_mat, i::Int, j::Int)
  ccall((:_nmod_mat_swap_rows, :libflint), Void, (Ptr{nmod_mat}, Int, Int), &x, i-1, j-1)
#  for k in 1:rows(x)
#    s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Lim, (Ptr{nmod_mat}, Int, Int), &x, i, k)
#    t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Lim, (Ptr{nmod_mat}, Int, Int), &x, j, k)
#    set_entry!(x, i, k, t)
#    set_entry!(y, j, k, s)
#  end
  nothing
end

function _swaprows!(x::nmod_mat)
  r = rows(x)
  c = cols(x)

  if r == 1
    return nothing
  end

  if r % 2 == 0
    for i in 1:div(r,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, i - 1, j - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, (r - i + 1) - 1, j - 1)
        set_entry!(x, r - i + 1, j, s)
        set_entry!(x, i, j, t)
      end
    end
  else
    for i in 1:div(r-1,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, i - 1, j - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, (r - i + 1) - 1, j - 1)
        set_entry!(x, i, j, t)
        set_entry!(x, r - i + 1, j, s)
      end
    end
  end
  nothing
end

function _swapcols!(x::nmod_mat)
  r = rows(x)
  c = cols(x)

  if c == 1
    return nothing
  end

  if c % 2 == 0
    for i in 1:div(c,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, j - 1, i - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, j - 1, (c - i + 1 ) - 1)
        set_entry!(x, j, i, t)
        set_entry!(x, j, c - i + 1, s)
      end
    end
  else
    for i in 1:div(c-1,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, j - 1, i - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &x, j - 1, (c - i + 1 ) - 1)
        set_entry!(x, j, i, t)
        set_entry!(x, j, c - i + 1, s)
      end
    end
  end
  nothing
end

function _swapcols!(x::fmpz_mat)
  r = rows(x)
  c = cols(x)

  if c == 1
    return x
  end

  if c % 2 == 0
    for i in 1:div(c,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, j - 1, i - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, j - 1, (c - i + 1 ) - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  else
    for i in 1:div(c-1,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, j - 1, i - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &x, j - 1, (c - i + 1 ) - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  end
  nothing
end

function _swapcols(x::GenMat)
  z = deepcopy(x)
  _swapcols!(z)
  return z
end

function _swapcols!(x::GenMat)
  r = rows(x)
  c = cols(x)
  t = base_ring(x)(0)

  if c == 1
    return x
  end

  if c % 2 == 0
    for i in 1:div(c,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        x[j, i], x[j, c - i + 1] = x[j, c - i + 1], x[j, i]
      end
    end
  else
    for i in 1:div(c-1,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        x[j, i], x[j, c - i + 1] = x[j, c - i + 1], x[j, i]
      end
    end
  end
  nothing
end

function _swaprows(x::GenMat)
  z = deepcopy(x)
  _swaprows(z)
  return z
end

function _swaprows!(x::GenMat)
  r = rows(x)
  c = cols(x)

  if r == 1
    return x
  end

  if r % 2 == 0
    for i in 1:div(r,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        x[i, j], x[r - i + 1, j] = x[r - i + 1, j], x[i, j]
      end
    end
  else
    for i in 1:div(r-1,2)
      for j = 1:c
        x[i, j], x[r - i + 1, j] = x[r - i + 1, j], x[i, j]
        # we swap x[i,j] <-> x[r-i+1,j]
      end
    end
  end
  nothing
end
################################################################################
#
################################################################################

function abs_max(a::fmpz_mat)
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &a, 0,0)
  for i=1:rows(a)
    for j=1:cols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &a, i-1, j-1)
      if ccall((:fmpz_cmpabs, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) < 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_abs, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), &r, m)
  return r
end

function max(a::fmpz_mat)  #should be maximum in julia
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &a, 0,0)
  for i=1:rows(a)
    for j=1:cols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &a, i-1, j-1)
      if ccall((:fmpz_cmp, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) < 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_set, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), &r, m)
  return r
end

function min(a::fmpz_mat)  #should be minimum in julia
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &a, 0,0)
  for i=1:rows(a)
    for j=1:cols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &a, i-1, j-1)
      if ccall((:fmpz_cmp, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) > 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_set, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), &r, m)
  return r
end

function lift_unsigned(a::nmod_mat)
  z = MatrixSpace(FlintZZ, rows(a), cols(a))()
  ccall((:fmpz_mat_set_nmod_mat_unsigned, :libflint), Void,
          (Ptr{fmpz_mat}, Ptr{nmod_mat}), &z, &a)
  return z
end

################################################################################
# possibly a slice or window in fmpz_mat?
# the nr x nc matrix starting in (a,b)
################################################################################

function submat(A::fmpz_mat, a::Int, b::Int, nr::Int, nc::Int)
  @assert nr >= 0 && nc >= 0
  M = MatrixSpace(FlintZZ, nr, nc)()::fmpz_mat
  t = ZZ()
  for i = 1:nr
    for j = 1:nc
      getindex!(t, A, a+i-1, b+j-1)
      M[i,j] = t
    end
  end
  return M
end

function submat{T <: Integer}(A::fmpz_mat, r::UnitRange{T}, c::UnitRange)
  @assert !isdefined(r, :step) || r.step==1
  @assert !isdefined(c, :step) || c.step==1
  return submat(A, r.start, c.start, r.stop-r.start+1, c.stop-c.start+1)::fmpz_mat
end

function sub(A::fmpz_mat, r::UnitRange, c::UnitRange)
  @assert !isdefined(r, :step) || r.step==1
  @assert !isdefined(c, :step) || c.step==1
  return submat(A, r.start, c.start, r.stop-r.start+1, c.stop-c.start+1)::fmpz_mat
end

################################################################################
#
#  misc.jl : At some point this should migrate to Nemo
#
################################################################################

# compute basis for the right kernel space by calling flint
# look at flint documentation of nmod_mat_nullspace

function _right_kernel(x::nmod_mat)
  z = MatrixSpace(base_ring(x), cols(x), max(rows(x),cols(x)))()
  n = ccall((:nmod_mat_nullspace, :libflint), Clong, (Ptr{nmod_mat}, Ptr{nmod_mat}), &z, &x)
  return z,n
end

# compute basis for the left kernel space
# output is array of arrays, which span the kernel

function kernel(a::nmod_mat)
  x = transpose(a)
  z,n = _right_kernel(x)
  z = transpose(z)
  #println(z)
  ar = typeof(Array(GenRes{fmpz}, cols(z)))[]
  for i in 1:n
    t = Array(GenRes{fmpz}, cols(z))
    for j in 1:cols(z)
      t[j] = z[i,j]
    end
    push!(ar,t)
  end
  #println(ar)
  return ar
end

function kernel_mod(a::fmpz_mat, m::fmpz)
  b = deepcopy(a)
  for i in 1:rows(b)
    for j in 1:cols(b)
      b[i,j] = b[i,j] % m
    end
  end

  # fmpz_mat_rref_mod assumes that input is reduced modulo m
  r = ccall((:fmpz_mat_rref_mod, :libflint), Int, (Ptr{Void}, Ptr{fmpz_mat}, Ptr{fmpz}), C_NULL, &b, &m)
  pivots = Array(Int, r)
  nonpivots = Array(Int, cols(b) - r)
  X = zero(MatrixSpace(ZZ,cols(b),cols(b)))

  if r == 0
    return one(MatrixSpace(ZZ, cols(b), cols(b)))
  elseif !((cols(b) - r) == 0)
    i = 1
    j = 1
    k = 1
    for i in 1:r
      while b[i,j] == 0
        nonpivots[k] = j
        k += 1
        j += 1
      end
      pivots[i] = j
      j += 1
    end
    while k <= cols(b) - r
      nonpivots[k] = j
      k += 1
      j += 1
    end

    for i in 1:cols(b) - r
      for j in 1:r
        X[pivots[j],i] = - ZZ(b[j,nonpivots[i]])
      end
      X[nonpivots[i],i] = ZZ(1)
    end
  end
  return X, r
end

# Another kernel function
function _kernel(x::fmpz_mat)
  H, U = hnf_with_transform(x)
  i = 1
  for i in 1:rows(H)
    if is_zero_row(H, i)
      break
    end
  end
  return submat(U, i:rows(U), 1:cols(U))
end

################################################################################
#
#  Copy matrix into another matrix
#
################################################################################

# Copy B into A at position (i, j)
function _copy_matrix_into_matrix(A::fmpz_mat, i::Int, j::Int, B::fmpz_mat)
  for k in 0:rows(B) - 1
    for l in 0:cols(B) - 1
      d = ccall((:fmpz_mat_entry, :libflint),
                Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &B, k, l)
      t = ccall((:fmpz_mat_entry, :libflint),
                Ptr{fmpz}, (Ptr{fmpz_mat}, Int, Int), &A, i - 1 + k, j - 1 + l)
      ccall((:fmpz_set, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, d)
    end
  end
end

function swap_rows!(M::MatElem, i::Int, j::Int)
  for k in 1:cols(M)
    t = M[i, k]
    M[i, k] = M[j, k]
    M[j, k] = t
  end
end


text:
@testset "FFT of geopotential" begin

    P = Params(NF=Float64)
    G = GeoSpectral{P.NF}(P)
    B = Boundaries{P.NF}(P,G)

    @test all(B.ϕ0trunc .≈ fourier_inverse(fourier(B.ϕ0trunc,G),G))
end

@testset "Legendre transform of geopotential" begin

    P = Params(NF=Float64)
    G = GeoSpectral{P.NF}(P)
    B = Boundaries{P.NF}(P,G)

    # using the already spectrally truncated geopotential
    @test all(B.ϕ0trunc .≈ gridded(spectral(B.ϕ0trunc,G),G))
end

@testset "Spectral transform of spectral noise" begin
    # original speedy T30 with 96x48 grid
    P = Params(NF=Float64;trunc=30,nlat=48,nlon=96)
    G = GeoSpectral{P.NF}(P)
    B = Boundaries{P.NF}(P,G)

    mx = G.spectral.mx
    nx = G.spectral.nx

    A = rand(mx,nx) + im*rand(mx,nx)
    At = spectral(gridded(A,G),G)   # the first transform includes truncation

    # the next should be exact to machine precision
    At2 = spectral(gridded(At,G),G)
    At3 = spectral(gridded(At2,G),G)
    @test all(At .≈ At2)
    @test all(At2 .≈ At3)
end


text:
using Dates
using TimeZones

"""
Parses the date/time part of the request URL for cases such as querytimeseries where a start- and end date & time plus
a time-step is required. Calls parse_datetime (singular) to parse the start- and end date & time; calls parse_interval
to parse the time-step.
"""
function parse_datetimes(startdate, enddate, interval)
    return "$(parse_datetime(startdate))--$(parse_datetime(enddate)):$(parse_interval(interval))"
end

"""
Parses date & time information passed as a string. Checks that the format of the string is as expected.
"""
function parse_datetime(datetime::String)
    @assert occursin(r"^\d\d\d\d-\d\d-\d\dT\d\d:\d\d:\d\dZ$", datetime) "Make sure start- and end datetime strings match 'YYYY-MM-DDTHH:MM:SST' format, where T and Z are literal characters."
    return "$datetime"
end

"""
Parses date & time information passed as a Date.DateTime.
"""
function parse_datetime(datetime::DateTime)
    return "$(Dates.format(datetime, "yyyy-mm-ddTHH:MM:SS"))Z"
end

"""
Parses date & time information passed as a TimeZones.ZonedDateTime by converting to a Dates.DateTime in UTC
and calling the method for a Dates.DateTime argument.
"""
function parse_datetime(datetime::ZonedDateTime)
    return parse_datetime(DateTime(datetime, UTC))
end

"""
Parses the time-step part of the URL if the user has provided a string. Checks that the string has the correct format.
"""
function parse_interval(interval::String)
    @assert occursin(r"^P(\dD)?(T(\d[HMS])+)?$", interval) "Intervals should start with P. Days may be left blank, but sub-day periods should be separated from days by T."
    return "$interval"
end

"""
Parses the time-step part of the URL if the user has provided a (number of) Dates.Day
"""
function parse_interval(interval::Dates.Day)
    return "P$(interval.value)D"
end

"""
Parses the time-step part of the URL if the user has provided a (number of) Dates.Hour
"""
function parse_interval(interval::Dates.Hour)
    return "PT$(interval.value)H"
end

"""
Parses the time-step part of the URL if the user has provided a (number of) Dates.Minute
"""
function parse_interval(interval::Dates.Minute)
    return "PT$(interval.value)M"
end

"""
Parses the time-step part of the URL if the user has provided a (number of) Dates.Second
"""
function parse_interval(interval::Dates.Second)
    return "PT$(interval.value)S"
end

"""
Parses the time-step part of the URL if the user has provided a Dates.CompoundPeriod.
"""
function parse_interval(interval::Dates.CompoundPeriod)
    lookup = Dict(Dates.Day=>'D', Dates.Hour=>'H', Dates.Minute=>'M',Dates.Second=>'S')
    undelimited = join(["$(p.value)$(lookup[typeof(p)])" for p in interval.periods])
    ind = findfirst(c->c=='D', undelimited)
    if isnothing(ind)
        daypart=""
        timepart=undelimited
    else
        daypart="$(undelimited[1:ind])"
        timepart="$(undelimited[ind+1:end])"
    end
    return "P$(daypart)T$(timepart)"
end

"""
Parses the parameters part of the URL, assuming the user has provided a vector containing their parameters.
Concatenates the parameters into a string and passes to string version of the method.
"""
function parse_parameters(parameters::Vector{String})
    return parse_parameters(join(parameters, ','))
end

"""
Checks that the parameter string provided (either by the user or by the vector method) fits the expected pattern.
"""
function parse_parameters(parameters::String)
    @assert occursin(r"(^.+:.+)+$", parameters) "There may be a typo in one or more of your parameters. Parameters should feature a name and a unit, separated by ':'"
    return parameters
end

"""
Parses the location part of the URL if the user has provided a vector of coordinate tuples as required for e.g.
querytimeseries. Checks each coordinate pair to ensure that they are within an accepted range (wrapping longitudes
if they are in the 180-360 degree range).
"""
function parse_locations(coordinate_list::Vector{<:NTuple{2, Real}})
    coordinate_string = ""
    for (lat, lon) in coordinate_list
        if abs(lat) > 90 throw(DomainError(lat, "value literally out of this world, in a bad way")) end
        if -180 > lon || lon > 360 throw(DomainError(lon, "value outside of supported range"))
        elseif lon > 180 lon -= 360 end
        coordinate_string*="$lat,$lon+"
    end
    return first(coordinate_string, length(coordinate_string)-1)
end

"""
The user may feasibly forget to wrap a single location tuple in a vector, in which case this method saves the day.
"""
function parse_locations(coordinate::Tuple{Real, Real})
    lat = coordinate[1]
    lon = coordinate[2]
    if abs(lat) > 90 throw(DomainError(lat, "value literally out of this world, in a bad way")) end
    if -180 > lon || lon > 360 throw(DomainError(lon, "value outside of supported range"))
    elseif lon > 180 lon -= 360 end
    return "$lat,$lon"
end

"""
The user may feasibly enter a single location as a vector, rather than a vector of tuples. This method permits this mistake
"""
function parse_locations(coordinate::Vector{T} where T<:Real)
    if length(coordinate) != 2 throw(ArgumentError("Please wrap your coordinate pairs as Tuples")) end
    lat = coordinate[1]
    lon = coordinate[2]
    if abs(lat) > 90 throw(DomainError(lat, "value literally out of this world, in a bad way")) end
    if -180 > lon || lon > 360 throw(DomainError(lon, "value outside of supported range"))
    elseif lon > 180 lon -= 360 end
    return "$lat,$lon"
end

"""
Parses the location part of the URL if the user has provided bounding latitudes and longitudes as required for
e.g. querygrid. Checks each coordinate  to ensure that they are within an accepted range (wrapping longitudes
if they are in the 180-360 degree range).
"""
function parse_locations(north::Real, west::Real, south::Real, east::Real, reslat::Real, reslon::Real)
    if abs(north) > 90 throw(DomainError(north, "value literally out of this world, in a bad way")) end
    if abs(south) > 90 throw(DomainError(south, "value literally out of this world, in a bad way")) end
    if east < west throw(ArgumentError("`east` should be greater than `west`. Areas crossing the dateline are not currently supported")) end
    if -180 > west || west > 360 throw(DomainError(north, "value outside of supported range"))
    elseif west > 180 west -= 360 end
    if -180 > east || east > 360 throw(DomainError(east, "value outside of supported range"))
    elseif east > 180 east -= 360 end
    return "$(north),$(west)_$(south),$(east):$reslat,$reslon"
end

"""
Unpacks a DataTypes.data_url into the format anticipated for a request.
"""
function parse_url(username::String, password::String, url::data_url)
    request_url = "https://$username:$password@api.meteomatics.com/"
    request_url *= "$(url.validdatetime)/$(url.parameters)/$(url.location)/$(url.format)"
    request_url *= "?"
    for (key, value) in url.opts
        if value === nothing
        else
            request_url *= "$key=$value&"
        end
    end
    return first(request_url, length(request_url) - 1)
end

"""
Unpacks a DataTypes.meta_url into the format anticipated for a request.
"""
function parse_url(username::String, password::String, url::meta_url)
    base_url = "https://$username:$password@api.meteomatics.com/"
    base_url *= url.metaquery*"?"
    base_url *= "model=$(url.model)&"
    if !isnothing(url.validdatetime) base_url *= "valid_date=$(url.validdatetime)&" end
    base_url *= "parameters=$(url.parameters)"
    return base_url
end

max_stars_repo_name: meteomatics/julia-connector-api

max_issues_repo_name: meteomatics/julia-connector-api

max_forks_repo_name: meteomatics/julia-connector-api


text:
using Test, QBase

using LinearAlgebra

@testset "./src/types/brakets.jl" begin

@testset "is_braket()" begin
    @testset "simple valid wavefunction" begin
        @test is_braket([1;0])
        @test is_braket([1])
        @test is_braket([0,-im]')
        @test is_braket([1;1;-1]/sqrt(3))
        @test is_braket([1;im]/sqrt(2))
    end

    @testset "invalid wavefunctions" begin
        @test !is_braket([0;0])
        @test !is_braket([1;1])
        @test_throws MethodError is_braket([1 0;0 0])
        @test_throws MethodError is_braket([1 0;0 0]')
    end

    @testset "valid input vector types" begin
        @test is_braket([1,0])
        @test is_braket([1,0]')
        @test is_braket([1;0])
        @test is_braket(Ket([1,0]))
        @test is_braket(Bra([1,0]))
    end

    @testset "tuning atol" begin
        ϵ = 1e-7

        @test is_braket([1,1], atol=1)
        @test is_braket([1,1], atol=sqrt(2)-1+ϵ)
        @test !is_braket([1,1], atol=sqrt(2)-1-ϵ)
        @test !is_braket([1,1], atol=0)
    end
end

@testset "_braket_error()" begin
    @test_throws DomainError QBase._braket_error([1,0])

    try
        QBase._braket_error([1,0])
        @test false
    catch err
        @test err isa DomainError
        @test err.val == [1, 0]
        @test err.msg ==  "ψ must be a vector that satisfies `norm(ψ,2) != 1`."
    end
end

@testset "Ket()" begin
    @testset "valid input types" begin
        ψ = Ket([1,0])
        @test ψ isa Ket{Int64}
        @test ψ isa Ket{Int64}
        @test ψ == [1,0]

        ψ = Ket([1.,0.])
        @test ψ isa Ket{Float64}
        @test ψ isa Ket{Float64}
        @test ψ == [1,0]

        ψ = Ket([1//1,0//1])
        @test ψ isa Ket{Rational{Int64}}
        @test ψ isa Ket{Rational{Int64}}
        @test ψ == [1,0]

        ψ = Ket([im,0])
        @test ψ isa Ket{Complex{Int64}}
        @test ψ isa Ket{Complex{Int64}}
        @test ψ == [im,0]

        ψ = Ket([im,0.])
        @test ψ isa Ket{Complex{Float64}}
        @test ψ isa Ket{Complex{Float64}}
        @test ψ == [im,0]
    end

    @testset "simple test cases" begin
        valid_cases = [
            [1,0,0],
            [1,0,0,0],
            [1,1,im]/sqrt(3),
            [0.5im,0//1,-sqrt(3)/2]
        ]

        for case in valid_cases
            @test Ket(case) == case
        end
    end

    @testset "input vector formats" begin
        ψ_match = [1,0]

        ψ = Ket([1;0] :: Vector{Int64})
        @test ψ == ψ_match

        ψ = Ket([1,0]' :: Adjoint{Int64,Vector{Int64}})
        @test ψ == ψ_match

        ψ = Ket([1,0] :: Vector{Int64})
        @test ψ == ψ_match

        ψ = Ket(Bra([im,0]'))
        @test ψ == [im,0]
    end

    @testset "invalid inputs" begin
        @test_throws DomainError Ket([1,1,1])
        @test_throws MethodError Ket(["a","b","c"])
    end

    @testset "tuning atol" begin
        ϵ = 1e-7
        ψ = Ket([1,1], atol=sqrt(2)-1+ϵ)

        @test ψ.atol ≈ sqrt(2)-1+ϵ
        @test ψ isa Ket{Int64}
        @test_throws DomainError Ket([1,1], atol=sqrt(2)-1-ϵ)
    end
end

@testset "Bra()" begin
    @testset "valid input types" begin
        ψ = Bra([1,0])
        @test ψ isa Bra{Int64}
        @test ψ isa Bra{Int64}
        @test ψ == [1,0]'

        ψ = Bra([1.,0.])
        @test ψ isa Bra{Float64}
        @test ψ isa Bra{Float64}
        @test ψ == [1,0]'

        ψ = Bra([1//1,0//1])
        @test ψ isa Bra{Rational{Int64}}
        @test ψ isa Bra{Rational{Int64}}
        @test ψ == [1,0]'

        ψ = Bra([im,0])
        @test ψ isa Bra{Complex{Int64}}
        @test ψ isa Bra{Complex{Int64}}
        @test ψ == [im,0]'

        ψ = Bra([im,0.])
        @test ψ isa Bra{Complex{Float64}}
        @test ψ isa Bra{Complex{Float64}}
        @test ψ == [im,0]'
    end

    @testset "simple test cases" begin
        valid_cases = [
            [1,0,0]',
            [1,0,0,0]',
            [1,1,im]'/sqrt(3),
            [0.5im,0//1,-sqrt(3)/2]'
        ]

        for case in valid_cases
            @test Bra(case) == case
        end
    end

    @testset "input vector formats" begin
        ψ_match = [1 0]

        ψ = Bra([1;0] :: Vector{Int64})
        @test ψ == ψ_match

        ψ = Bra([1,0]' :: Adjoint{Int64,Vector{Int64}})
        @test ψ == ψ_match

        ψ = Bra([1,0] :: Vector{Int64})
        @test ψ == ψ_match

        ψ = Bra(Ket([im,0]))
        @test ψ == [-im 0]

        ψ = Bra([im,0])
        @test ψ == [im,0]'
    end

    @testset "invalid inputs" begin
        @test_throws DomainError Bra([1,1,1])
        @test_throws MethodError Bra(["a","b","c"])
    end

    @testset "tuning atol" begin
        ϵ = 1e-7
        ψ = Bra([1,1], atol=sqrt(2)-1+ϵ)

        @test ψ.atol ≈ sqrt(2)-1+ϵ
        @test ψ isa Bra{Int64}
        @test_throws DomainError Bra([1,1], atol=sqrt(2)-1-ϵ)
    end
end

@testset "bra and ket multiplication" begin
    @testset "product of Complex{Int64} types" begin
        ket = Ket([im,0])
        bra = Bra([im,0])

        inner_product = bra*ket
        @test inner_product isa Complex{Int64}
        @test inner_product == 1

        outer_product = ket*bra
        @test outer_product isa Matrix{Complex{Int64}}
        @test outer_product == [1 0;0 0]
    end

    @testset "product of Complex{Int64} types" begin
        ket = Ket([1,0])
        bra = Bra([1,im]/sqrt(2))

        inner_product = bra*ket
        @test inner_product isa Complex{Float64}
        @test inner_product ≈ 1/sqrt(2)

        outer_product = ket*bra
        @test outer_product isa Matrix{Complex{Float64}}
        @test outer_product ≈ [1 -im;0 0]/sqrt(2)
    end
end

@testset "bra and ket conversion via adjoint" begin
    ket = Ket([im,0])
    bra = Bra([im,0]')

    ψ = ket'
    @test ψ isa Bra{Complex{Int64}}
    @test ψ == [-im 0]

    ψ = adjoint(ket)
    @test ψ isa Bra{Complex{Int64}}
    @test ψ == [-im 0]

    ψ = bra'
    @test ψ isa Ket{Complex{Int64}}
    @test ψ == [im,0]

    ψ = adjoint(bra)
    @test ψ isa Ket{Complex{Int64}}
    @test ψ == [im,0]
end

@testset "kron() for bras and kets" begin
    @testset "ket" begin
        ket0 = Ket([1,0])
        ket1 = Ket([0,1])

        ket01 = kron(ket0,ket1)
        @test ket01 isa Ket{Int64}
        @test ket01 == [0,1,0,0]

        ket010 = kron(ket0,ket1,ket0)
        @test ket010 isa Ket{Int64}
        @test ket010 == [0,0,1,0,0,0,0,0]
    end

    @testset "ket using atol" begin
        ϵ = 1e-7
        ket0 = Ket([1,1], atol=sqrt(2)-1+ϵ)
        ket1 = Ket([1,-1], atol=sqrt(2)-1+ϵ)

        ket01 = kron(ket0,ket1, atol=2-1+ϵ)
        @test ket01 isa Ket{Int64}
        @test ket01 == [1,-1,1,-1]

        @test_throws DomainError kron(ket0,ket1, atol=2-1-ϵ)
        @test_throws DomainError kron(ket0,ket1)
    end

    @testset "bra" begin
        bra0 = Bra([1,0])
        bra1 = Bra([0,1])

        bra01 = kron(bra0,bra1)
        @test bra01 isa Bra{Int64}
        @test bra01 == [0,1,0,0]'

        bra010 = kron(bra0,bra1,bra0)
        @test bra010 isa Bra{Int64}
        @test bra010 == [0,0,1,0,0,0,0,0]'
    end

    @testset "bra using atol" begin
        ϵ = 1e-7
        bra0 = Bra([1,1], atol=sqrt(2)-1+ϵ)
        bra1 = Bra([1,-1], atol=sqrt(2)-1+ϵ)

        bra01 = kron(bra0,bra1, atol=2-1+ϵ)
        @test bra01 isa Bra{Int64}
        @test bra01 == [1 -1 1 -1]

        @test_throws DomainError kron(bra0,bra1, atol=2-1-ϵ)
        @test_throws DomainError kron(bra0,bra1)
    end
end

end


text:
using JuMP, Ipopt, Plots

function blockMoving(numCollocationPoints)
    model = Model(Ipopt.Optimizer)

    # Timestep
    hk = 30/(numCollocationPoints-1) # for now set as constant, later add variable timestep
    time = 0:hk:30

    # Variables and start values
    @variable(model, x[i=1:numCollocationPoints])
    set_start_value.(x,time)
    @variable(model, v[1:numCollocationPoints], start = 1)
    @variable(model, u[1:numCollocationPoints], start = 0)

    # Boundary constraints
    fix(x[1], 0)
    fix(x[end], 30)
    fix(v[1], 0)
    fix(v[end], 0)

    # System dynamics
    @NLexpression(model, xDot[i=1:numCollocationPoints], v[i])
    @NLexpression(model, vDot[i=1:numCollocationPoints], u[i])

    # Objective function
    @NLexpression(model, objective[i=1:numCollocationPoints], u[i]^2)
    for i in 1:numCollocationPoints - 1
        # Collocation points
        @NLconstraint(model, x[i+1] - x[i] - 0.5 * hk * (xDot[i+1] + xDot[i]) == 0)
        @NLconstraint(model, v[i+1] - v[i] - 0.5 * hk * (vDot[i+1] + vDot[i]) == 0)
    end

    # Objective function interpolation
    @NLexpression(model, objectiveInterp, sum(0.5 * hk * (objective[i] + objective[i+1]) for i in 1:numCollocationPoints - 1))
    @NLobjective(model, Min, objectiveInterp)
    optimize!(model)
    return value.(x), value.(v), value.(u), time
end

x,v,u,times = blockMoving(3000)

# Plotting
plotly()
plot(times, x, xlabel="time (s)", ylabel="x (m)", title="Location plot")
plot(times, v, xlabel="time (s)", ylabel="velocity (m/s)", title="Velocity plot")
plot(times, u, xlabel="time (s)", ylabel="Force (N)", title="Force plot")

max_stars_repo_name: sprockmonty/JuliaFYP

max_issues_repo_name: sprockmonty/JuliaFYP

max_forks_repo_name: sprockmonty/JuliaFYP


text:
struct Sampled{T,S} <: Evolution
    P::T
    rng::S
    Sampled(P::T, rng=Random.GLOBAL_RNG) where {T} = new{T,typeof(rng)}(P, rng)
end
eltype(::Type{Sampled{T}}) where {T} = eltype(T)

function evolve(X::Sampled, (t,x)::Pair, args...)
    u = evolve(X.P, (t=>x), args...)
    u === nothing && return nothing
    t, D = u
    x = rand(X.rng, D)
    t => x
end

struct Randn{T}
    x::T
end
eltype(::Type{Randn{T}}) where {T} = T
Randn(D::Randn) = D
#Randn(x::T) where {T} = Randn{T}(x)
rand(D::Randn) = randn(typeof(D.x))
rand(D::Randn{Array}) = randn(T, size(D.x))
rand(rng::AbstractRNG, D::Randn{T}) where {T} = randn(rng, typeof(D.x))
rand(rng::AbstractRNG, D::Randn{Array{T}}) where {T} = randn(rng, T, size(D.x))

struct WhiteNoise <: Evolution
end
evolve(::WhiteNoise, (t,x)::Pair{Int}) = (t+1) => Randn(x)
evolve(::WhiteNoise, x::Union{Randn, AbstractArray, Number}) where {T} = Randn(x)

#evolve(E::Evolution, (x,rng)::Sample{<:Start}) = evolve(E, Sample(x.value, rng))
#evolve(E::Evolution, (x,rng)::Sample) = rand(rng, evolve(E, x))
#evolve(E::Evolution, ((t, x), rng)::Sample{<:Pair{Int}}) = (t+1) => rand(rng, evolve(E, x))
dyniterate(E::Evolution, (x,rng)::Sample{<:Start}) = dyniterate(E, Sample(x.value, rng))
dyniterate(E::Evolution, u::Sample) = rand(u.rng, evolve(E, u.x)), u
function dyniterate(E::Evolution, ((t, x), rng)::Sample{<:Pair{Int}})
    y = rand(rng, evolve(E, x))
    (t+1) => y, Sample(t + 1 => y, rng)
end

#abstract type MarkovIterator
#end

struct InhomogeneousPoisson{T,S,R} <: Evolution
    λ::S
    λmax::T
    rng::R
end
InhomogeneousPoisson(λ, λmax) = InhomogeneousPoisson(λ, λmax, Random.GLOBAL_RNG)

function evolve(P::InhomogeneousPoisson, (t, i)::Pair)
    while true
        t = t - log(rand(P.rng))/P.λmax
        if rand(P.rng) ≤ P.λ(t)/P.λmax
            return t => i + 1
        end
    end
end

struct MetropolisHastings{T,F,G,RNG} <: Evolution
    P::T
    proposal::F
    logpdf::G
    rng::RNG
end
MetropolisHastings(P, proposal, logpdf=logpdf) = MetropolisHastings(P, proposal, logpdf, Random.GLOBAL_RNG)

function evolve(MH::MetropolisHastings, x)
    P = MH.P
    Q = MH.proposal(x)
    xᵒ = rand(Q)
    Qᵒ = MH.proposal(xᵒ)
    if log(rand(MH.rng)) < MH.logpdf(P, xᵒ) - MH.logpdf(P, x) + MH.logpdf(Qᵒ, x) - MH.logpdf(Q, xᵒ)
        x = xᵒ
    end
    x
end

#=
struct Link{S, T} <: DynamicIterator
    X::S
    Y::T
end
function dyniterate(P::Link, (p, q))
    u, p = @returnnothing dyniterate(P.X, p)
    v, q = @returnnothing dyniterate(P.Y, q, (control = u),)
    v, (p, q)
end

function dyniterate(S::Sampled{Link}, (p, q))
    u, p = @returnnothing dyniterate(Sampled(S.P.X, S.rng), p)
    v, q = @returnnothing dyniterate(Sampled(S.P.Y, S.rng), q, (control = u,))
    v, (p, q)
end
=#


text:
# --------------------------------------------------
# Cost Functions
# --------------------------------------------------
function cross_entropy_cost(l::Layer, D::Signal, Y::Signal)  # D = 教師
    cost = 0.0
    for n in 1:size(D,2)
        for k in 1:size(Y,1)
            cost += D[k,n] * log(Y[k,n])
        end
    end
    -(cost / size(D,2))
end

function cross_entropy_cost(diff::Symbol, l::Layer, D::Signal, Y::Signal)
    copy!(l.ΔE, (Y - D))
end

# --------------------------------------------------
function quadratic_cost(l::Layer, D::Signal, Y::Signal)      # D = 教師
    norm(Y - D)^2 / 2
end

function quadratic_cost(diff::Symbol, l::Layer, D::Signal, Y::Signal)
    copy!(l.ΔE, (Y - D))
end

# --------------------------------------------------
function rnn_quadratic_cost(l::Layer, D::Signal, Y::Signal)      # D = 教師
    norm(Y[:,end] - D[:,end])^2 / 2
end

function rnn_quadratic_cost(diff::Symbol, l::Layer, D::Signal, Y::Signal)
    copy!(l.ΔE, zeros(Float32, size(Y)))
    l.ΔE[:,end] = Y[:,end] - D[:,end]
    l.ΔE
end

# --------------------------------------------------
function rnn_cross_entropy_cost(l::Layer, D::Signal, Y::Signal)  # D = 教師
    cost = 0.0
    for k in 1:size(Y,1)
        cost += D[k,end] * log(Y[k,end])
    end
    -(cost)
end

function rnn_cross_entropy_cost(diff::Symbol, l::Layer, D::Signal, Y::Signal)
    copy!(l.ΔE, zeros(Float32, size(Y)))
    l.ΔE[:,end] = Y[:,end] - D[:,end]
    l.ΔE
end


text:
using Aqua
using Test
using LegendrePolynomials

@testset "Project quality assurance" begin
    Aqua.test_all(LegendrePolynomials, ambiguities = false)
end

using PerformanceTestTools

PerformanceTestTools.@include_foreach(
    "test_threaded.jl",
    [nothing,
    ["JULIA_NUM_THREADS" => Threads.nthreads() > 1 ? "1" : "2"],
    ],
)


text:
# Playing Cards: Standard Deck
import Base: +, -
import Random: shuffle!

struct Deck
    cards::Array{Card, 1}
    function Deck()
        cards = Card[]
        for rank in 1:13
            for suit in 1:4
                push!(cards, Card(rank, suit))
            end
        end
        return new(cards);
    end
end

@inline Base.length(deck::Deck) = length(deck.cards);
@inline Base.firstindex(deck::Deck) = 1;
@inline Base.lastindex(deck::Deck) = length(deck);
@inline Base.getindex(deck::Deck, i::Int64) = deck.cards[i];
@inline Base.getindex(deck::Deck, I::UnitRange{Int64}) = [deck[i] for i in I];
@inline Base.getindex(deck::Deck, I::Array{Int64, 1}) = [deck[i] for i in I];
@inline Base.iterate(deck::Deck, i::Int64 = 1) = (i > length(deck)) ? nothing : (deck[i], i + 1);

function Base.show(io::IO, deck::Deck)
    for card in deck
        print(io, card, " ")
    end
    return nothing;
end

@inline +(deck::Deck, card::Card) = push!(deck.cards, card);
@inline +(deck::Deck, cards::Array{Card, 1}) = append!(deck.cards, cards);
@inline -(deck::Deck, card::Card) = try deleteat!(deck.cards, findall(x -> x === card, deck.cards)) catch end;
function -(deck::Deck, cards::Array{Card, 1})
    for card in cards
        deck - card
    end
    return deck;
end

function shuffle(deck::Deck)
    shuffle!(deck.cards)
    return deck;
end

@inline deal(deck::Deck) = popfirst!(deck.cards);
function deal(deck::Deck, n::Int64)
    @assert(1 ≤ n ≤ length(deck.cards), "Invalid number of cards!")
    cards = Card[]
    for _ in 1:n
        push!(cards, deal(deck))
    end
    return cards;
end

max_stars_repo_name: Joshua-Robison/PlayingCards

max_issues_repo_name: Joshua-Robison/PlayingCards

max_forks_repo_name: Joshua-Robison/PlayingCards