text:
# # Forward model
#
# Here, we go through an example of simulating the full forward model,
# from the prior definition to the expected number of events in different
# bins of the detector response.

using QCDNUM, PartonDensity
using Plots, Printf, NaNMath, Parameters

# ## Define input PDFs
# We can use `DirichletPDFParams` or `ValencePDFParams`, as long
# as we do so according to the *PDF parametrisation and priors* docs.

weights = [3., 1., 5., 5., 1., 1., 1., 0.5, 0.5]
pdf_params = DirichletPDFParams(K_u=4.0, K_d=6.0, λ_g1=0.7, λ_g2=-0.4,
                                K_g=6.0, λ_q=-0.5, seed=5, weights=weights);

# Plot the input PDFs

plot_input_pdfs(pdf_params)

# Sanity check that sum = 1

int_xtotx(pdf_params) ≈ 1

# ## Define QCDNUM grids, weights and settings

grid = QCDNUMGrid(x_min=[1.0e-3], x_weights=[1], nx=100,
                  qq_bounds=[1.0e2, 3.0e4], qq_weights=[1.0, 1.0], nq=50, spline_interp=3)

qcdnum_params = QCDNUMParameters(order=2, α_S=0.118, q0=100.0, grid=grid,
                                 n_fixed_flav=5, iqc=1, iqb=1, iqt=1, weight_type=1);

# Initialise and set key parameters

QCDNUM.qcinit(-6, "")
QCDNUM.setord(qcdnum_params.order)
QCDNUM.setalf(qcdnum_params.α_S, qcdnum_params.q0)

# Build grids

g = qcdnum_params.grid
QCDNUM.gxmake(g.x_min, g.x_weights, g.x_num_bounds, g.nx,
              g.spline_interp);
QCDNUM.gqmake(g.qq_bounds, g.qq_weights, g.qq_num_bounds, g.nq);

# Define FFNS/VFNS

QCDNUM.setcbt(qcdnum_params.n_fixed_flav, qcdnum_params.iqc,
              qcdnum_params.iqb, qcdnum_params.iqt);

# Build weight tables

nw = QCDNUM.fillwt(qcdnum_params.weight_type)
nw = QCDNUM.zmfillw()

# ## Evolve the PDFs using QCDNUM
#
# Define input PDF function
# * See https://www.nikhef.nl/~h24/qcdnum-files/doc/qcdnum170115.pdf under `evolfg`
#
# There are functions available to generate the necessary input PDF function in
# the correct format for `QCDNUM.jl` (see `get_input_pdf_func()`), along with
# the mapping between this input function and quark species (see `input_pdf_map`).

# Get function and wrap with c-style pointer

my_func = get_input_pdf_func(pdf_params)
input_pdf = @cfunction(my_func, Float64, (Ref{Int32}, Ref{Float64}))

# Find index of starting scale and evolve

iq0 = QCDNUM.iqfrmq(qcdnum_params.q0)
pdf_loc = 1
eps = QCDNUM.evolfg(pdf_loc, input_pdf, input_pdf_map, iq0)

# ## Define necessary splines for cross section calculation

# For splines

splint_params = SPLINTParameters();
quark_coeffs = QuarkCoefficients();

# Define initial spline

QCDNUM.ssp_spinit(splint_params.nuser);
ia = QCDNUM.isp_s2make(splint_params.nsteps_x, splint_params.nsteps_q);
xnd = QCDNUM.ssp_unodes(ia, splint_params.nnodes_x, 0);
qnd = QCDNUM.ssp_vnodes(ia, splint_params.nnodes_q, 0);

# Check nodes and erase

QCDNUM.ssp_nprint(ia);
QCDNUM.ssp_erase(ia);

# Set nodes and fill spline with structure function

iaFLup = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaFLup, pdf_loc, quark_coeffs.proup, 1, 0.0);

iaF2up = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF2up, pdf_loc, quark_coeffs.proup, 2, 0.0);

iaF3up = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF3up, pdf_loc, quark_coeffs.valup, 3, 0.0);

iaFLdn = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaFLdn, pdf_loc, quark_coeffs.prodn, 1, 0.0);

iaF2dn = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF2dn, pdf_loc, quark_coeffs.prodn, 2, 0.0);

iaF3dn = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF3dn, 1, quark_coeffs.valdn, 3, 0.0);

# store spline addresses
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F2up, Float64(iaF2up));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F2dn, Float64(iaF2dn));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F3up, Float64(iaF3up));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F3dn, Float64(iaF3dn));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.FLup, Float64(iaFLup));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.FLdn, Float64(iaFLdn));

my_func = get_input_xsec_func()
input_xsec = @cfunction(my_func, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))

# plot

xsec_on_grid = zeros(g.nx, g.nq);

for ix = 1:g.nx
    for iq = 1:g.nq
        xsec_on_grid[ix, iq] = _fun_xsec_i(ix, iq)
    end
end

qcdnum_x_grid = QCDNUM.gxcopy(g.nx)
qcdnum_qq_grid =  QCDNUM.gqcopy(g.nq)
p1 = heatmap(qcdnum_x_grid, qcdnum_qq_grid, NaNMath.log10.(xsec_on_grid[:, :]'))
plot(p1, xlabel="x", ylabel="q2",
     xaxis=:log, yaxis=:log)

#

plot(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 4]),
    label="Q2=141 (input scale)", lw=3)
plot!(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 22]), label="Q2=1152", lw=3)
plot!(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 35]), label="Q2=5233", lw=3)
plot!(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 41]), label="Q2=10523", lw=3)
plot!(xaxis=:log, legend=:bottomleft, xlabel="x",
    ylabel="log10(cross section spline input)", ylims=(-7, 5))

#

set_lepcharge(1)
iaF_eP = QCDNUM.isp_s2make(1, 2);
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F_eP, Float64(iaF_eP));
QCDNUM.ssp_s2fill(iaF_eP, input_xsec, splint_params.rscut);

set_lepcharge(-1)
iaF_eM = QCDNUM.isp_s2make(1, 2);
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F_eM, Float64(iaF_eM));
QCDNUM.ssp_s2fill(iaF_eM, input_xsec, splint_params.rscut);

# plot spline

spline = zeros(g.nx, g.nq);

for ix = 1:g.nx
    for iq = 1:g.nq
        spline[ix, iq] = QCDNUM.dsp_funs2(iaF_eP, qcdnum_x_grid[ix],
            qcdnum_qq_grid[iq], 1)
    end
end

p1 = heatmap(qcdnum_x_grid, qcdnum_qq_grid, NaNMath.log10.(spline[:, :]'))
plot(p1, xlabel="x", ylabel="q2",
    xaxis=:log, yaxis=:log)

# ## Integrate over the cross section spline and find expected events numbers

nbins = size(xbins_M_begin)[1]
IntXsec_eP = zeros(nbins);
IntXsec_eM = zeros(nbins);
for i in 1:nbins
    IntXsec_eP[i] = QCDNUM.dsp_ints2(iaF_eP, xbins_M_begin[i], xbins_M_end[i],
        q2bins_M_begin[i], q2bins_M_end[i], 318., 4);
    IntXsec_eM[i] = QCDNUM.dsp_ints2(iaF_eM, xbins_M_begin[i], xbins_M_end[i],
        q2bins_M_begin[i], q2bins_M_end[i], 318., 4);
end

# 1 for e-p and 0 for e+p

ePp = 0;
eMp = 1;

TM_eP = get_TM_elements(ePp);
TM_eM = get_TM_elements(eMp);

K_eP = get_K_elements(ePp);
K_eM = get_K_elements(eMp);

nbins_out = size(TM_eP)[2];

xsec_pred_eP = zeros(nbins_out);
xsec_pred_eM = zeros(nbins_out);

for j in 1:nbins_out

    for i in 1:nbins

        xsec_pred_eP[j] += TM_eP[i, j] * (1.0/K_eP[i]) * IntXsec_eP[i];
        xsec_pred_eM[j] += TM_eM[i, j] * (1.0/K_eM[i]) * IntXsec_eM[i];

    end

end

xsec_pred_eM


text:
using Documenter, DocumenterCitations, ProximalAlgorithms
using Literate

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"))

src_path = joinpath(@__DIR__, "src/")

Literate.markdown(joinpath(src_path, "getting_started.jl"), src_path, documenter=true)
Literate.markdown(joinpath(src_path, "custom_objectives.jl"), src_path, documenter=true)

makedocs(
    bib,
    modules=[ProximalAlgorithms],
    sitename="ProximalAlgorithms.jl",
    pages=[
        "Home" => "index.md",
        "User guide" => [
            "getting_started.md",
            "implemented_algorithms.md",
            "custom_objectives.md",
        ],
        "Bibliography" => "bibliography.md",
    ],
)

deploydocs(
    repo="github.com/JuliaFirstOrder/ProximalAlgorithms.jl.git",
)

ext: jl


text:
using Plots, LaTeXStrings, ProgressBars
function LogisticMap(x::Matrix{<:Number}, r::Vector{<:Number})
    return x = 4 .* r .* x .* (1 .- x)
end
rcount= 2001
r_s = collect(range(0.0, 1.0, length=rcount))
x_s = rand(rcount, 10)
anim = @animate for n in ProgressBar(1:70)
    global x_s = LogisticMap(x_s, r_s)
    scatter(r_s, x_s, ms=0.5, msw=0, xlims=(0.0, 1.0), ylims = (0, 1), marker=:rect,
        legend=false, color=:black, alpha=0.25, xlabel=L"r", ylabel=L"x_n",
        title="\$\\textrm{Logistic\\ Map},\\quad n=$n\$",
        dpi=120)
end
gif(anim, "bifurcation.gif", fps=10)


text:
using LightGraphs
using MatrixMarket
using SimpleWeightedGraphs

name = splitext(splitdir(ARGS[1])[2])[1]
A = mmread(ARGS[1])
A = spdiagm(diag(A)) - A
g = SimpleWeightedGraph(A)

tris = Int(sum(triangles(g))//3)
clustering = global_clustering_coefficient(g)
lc = local_clustering_coefficient(g)
avg_clustering = mean(lc)
#= diam = parallel_diameter(g) =#

println("graph,triangles,clustering,avg_clustering")
@printf "%s,%s,%s,%s" name tris clustering avg_clustering

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
using SafeTestsets

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = ( Sys.iswindows() && haskey(ENV,"APPVEYOR") )
const is_TRAVIS = haskey(ENV,"TRAVIS")

if GROUP == "All"
    @time @safetestset "Exact coloring via contraction" begin include("test_contraction.jl") end
    @time @safetestset "Greedy distance-1 coloring" begin include("test_greedy_d1.jl") end
    @time @safetestset "Greedy star coloring" begin include("test_greedy_star.jl") end
    @time @safetestset "Acyclic coloring" begin include("test_acyclic.jl") end
    @time @safetestset "Matrix to graph conversion" begin include("test_matrix2graph.jl") end
    @time @safetestset "AD using colorvec vector" begin include("test_ad.jl") end
    @time @safetestset "Integration test" begin include("test_integration.jl") end
    @time @safetestset "Special matrices" begin include("test_specialmatrices.jl") end
    @time @safetestset "Jac Vecs and Hes Vecs" begin include("test_jaches_products.jl") end
end

if GROUP == "GPU"
    @time @safetestset "GPU AD" begin include("test_gpu_ad.jl") end
end


text:
export LowRankIntegralOperator

immutable LowRankIntegralOperator{S<:Space,M<:Space,T} <: AbstractLowRankOperator{T}
    U::Vector{Fun{S,T}}
    V::Vector{Fun{M,T}}

    function LowRankIntegralOperator(U::Vector{Fun{S,T}},V::Vector{Fun{M,T}})
        @assert length(U) == length(V)
        @assert length(U) > 0
        ds=space(first(V))
        for k=2:length(V)
            @assert space(V[k])==ds
        end
        rs=space(first(U))
        for k=2:length(U)
            @assert space(U[k])==rs
        end
        new(U,V)
    end
end

LowRankIntegralOperator{S,M,T}(U::Vector{Fun{S,T}},V::Vector{Fun{M,T}})=LowRankIntegralOperator{S,M,T}(U,V)

LowRankIntegralOperator{S,M,T1,T2}(U::Vector{Fun{S,T1}},V::Vector{Fun{M,T2}})=LowRankIntegralOperator(convert(Vector{Fun{S,promote_type(T1,T2)}},U),convert(Vector{Fun{M,promote_type(T1,T2)}},V))

LowRankIntegralOperator(A::Fun,B::Fun)=LowRankIntegralOperator([A],[B])

Base.promote_rule{S,M,T,S1,M1,T1}(::Type{LowRankIntegralOperator{S,M,T}},::Type{LowRankIntegralOperator{S1,M1,T1}})=LowRankIntegralOperator{promote_type(S,S1),promote_type(M,M1),promote_type(T,T1)}
Base.rank(L::AbstractLowRankOperator)=length(L.U)

domainspace(L::LowRankIntegralOperator)=space(first(L.V))
rangespace(L::LowRankIntegralOperator)=space(first(L.U))
promoterangespace(L::LowRankIntegralOperator,sp::Space)=LowRankIntegralOperator(map(u->Fun(u,sp),L.U),L.V)

+(L::LowRankIntegralOperator)=LowRankIntegralOperator(+L.U,L.V)
-(L::LowRankIntegralOperator)=LowRankIntegralOperator(-L.U,L.V)

*(L::LowRankIntegralOperator,f::Fun)=sum(map((u,v)->u*dot(v,f),L.U,L.V))

+(A::LowRankIntegralOperator,B::LowRankIntegralOperator)=LowRankIntegralOperator([A.U;B.U],[A.V;B.V])
-(A::LowRankIntegralOperator,B::LowRankIntegralOperator)=LowRankIntegralOperator([A.U;-B.U],[A.V;B.V])


text:
@testset "approximate" begin
    # Only test two-dimensional approximation
    # One-dimensional approximation gets tested in bases.jl
    @test coeffs(approximate((x1,x2)->x1*x2, Monomials(2))) ≈ [0 0; 0 1]
end


text:
module LibCSFML

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("CSFML was not built properly. Please run Pkg.build(\"CSFML\").")
end
include(depsjl_path)
# Module initialization function
function __init__()
    check_deps()
end

using CEnum

const Ctm = Base.Libc.TmStruct
const Ctime_t = UInt
const Cclock_t = UInt

export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "config.jl"))

include(joinpath(@__DIR__, "System", "system_common.jl"))
include(joinpath(@__DIR__, "System", "system_api.jl"))

include(joinpath(@__DIR__, "Audio", "audio_common.jl"))
include(joinpath(@__DIR__, "Audio", "audio_api.jl"))

include(joinpath(@__DIR__, "Window", "window_common.jl"))
include(joinpath(@__DIR__, "Window", "window_api.jl"))

include(joinpath(@__DIR__, "Graphics", "graphics_common.jl"))
include(joinpath(@__DIR__, "Graphics", "graphics_api.jl"))

include(joinpath(@__DIR__, "Network", "network_common.jl"))
include(joinpath(@__DIR__, "Network", "network_api.jl"))

# export everything
foreach(names(@__MODULE__, all=true)) do s
   if startswith(string(s), "sf")
       @eval export $s
   end
end

end # module


text:
using Test

using CompositeTypes
using CompositeTypes.Indexing
using CompositeTypes.Display

include("test_types.jl")
include("test_components.jl")
include("test_display.jl")


text:
using Clang

include("../src/find_ros2.jl")

isheader(header) = occursin(r".+\.(h|hpp)$", header)

function wrap_roslib(lib::AbstractString)
    lib_include = joinpath(ROS_INCLUDE_PATH, lib)
    lib_headers = String[]
    for (root, dirs, files) in walkdir(lib_include)
        for file in files
            if isheader(file)
                push!(lib_headers, joinpath(lib_include, root, file))
            end
        end
    end

    wc = init(; headers = lib_headers,
                output_file = joinpath(@__DIR__, "lib$(lib)_api.jl"),
                common_file = joinpath(@__DIR__, "lib$(lib)_common.jl"),
                clang_includes = vcat(lib_include, CLANG_INCLUDE),
                clang_args = ["-I", ROS_INCLUDE_PATH],
                header_wrapped = (root, current)->root == current,
                header_library = x->"lib$lib",
                clang_diagnostics = true,
                )

    run(wc)
end

wrap_roslib.(WRAPPED_ROS2_LIBS)


text:
using StatsBase
import POMDPs: MDP, Policy

#
# Define a policy that will act as an crossentropy policy
#
type CrossentropyPolicy <: Policy
    action_map::Array{Float64,2}
    action_space::Vector{Int64}
end

#
# Default uniform intializer based on a games input
#
function CrossentropyPolicy(pomdp::MDP)
    action_count = n_actions(pomdp)
    default_map = ones(action_count, n_states(pomdp)) * 1.0 / action_count
    return CrossentropyPolicy(default_map, 1:action_count)
end

#
# Initialize policy with a predefined values
#
function CrossentropyPolicy(pomdp::MDP, values::Array{Float64,2})
    action_count = n_actions(pomdp)
    return CrossentropyPolicy(values, 1:action_count)
end

#
# Select random action from a policy based on a state and action weight
#
action(policy::CrossentropyPolicy, s::Int64) = sample(policy.action_space, Weights(policy.action_map[:, s]))

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestCLI
using Test
import JSON

import SemanticFlowGraphs: CLI

const data_dir = joinpath(@__DIR__, "data")
const test_py = false
const test_r = false

mktempdir() do dir
  # Visualize raw flow graph.
  inpath = joinpath(data_dir, "clustering_kmeans.raw.graphml")
  outpath = joinpath(dir, "clustering_kmeans.raw.dot")
  CLI.main(["visualize", inpath, "--out", outpath])
  @test isfile(outpath)

  # Convert raw flow graph to semantic flow graph.
  outpath = joinpath(dir, "clustering_kmeans.semantic.graphml")
  CLI.main(["enrich", inpath, "--out", outpath])
  @test isfile(outpath)

  # Visualize semantic flow graph.
  inpath = outpath
  outpath = joinpath(dir, "clustering_kmeans.semantic.dot")
  CLI.main(["visualize", inpath, "--out", outpath])
  @test isfile(outpath)
end

mktempdir() do dir
  # Export concepts as JSON.
  outpath = joinpath(dir, "concepts.json")
  CLI.main(["ontology", "json", "--no-annotations", "--out", outpath])
  @test isfile(outpath)
  @test JSON.parsefile(outpath) isa AbstractVector
end

mktempdir() do dir
  # Export concepts as RDF/OWL.
  outpath = joinpath(dir, "concepts.ttl")
  CLI.main(["ontology", "rdf", "--no-annotations", "--out", outpath])
  @test isfile(outpath)
end

if test_py
  import PyCall

  mktempdir() do dir
    # Record Python raw flow graph.
    inpath = joinpath(data_dir, "sklearn_clustering_kmeans.py")
    outpath = joinpath(dir, "sklearn_clustering_kmeans.raw.graphml")
    CLI.main(["record", inpath, "--out", outpath])
    @test isfile(outpath)
  end
end

if test_r
  import RCall

  mktempdir() do dir
    # Record R raw flow graph.
    inpath = joinpath(data_dir, "clustering_kmeans.R")
    outpath = joinpath(dir, "clustering_kmeans.raw.graphml")
    CLI.main(["record", inpath, "--out", outpath])
    @test isfile(outpath)
  end
end

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
# Autogenerated wrapper script for libcxxwrap_julia_jll for aarch64-apple-darwin-julia_version+1.8.0
export libcxxwrap_julia, libcxxwrap_julia_stl

JLLWrappers.@generate_wrapper_header("libcxxwrap_julia")
JLLWrappers.@declare_library_product(libcxxwrap_julia, "@rpath/libcxxwrap_julia.0.dylib")
JLLWrappers.@declare_library_product(libcxxwrap_julia_stl, "@rpath/libcxxwrap_julia_stl.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcxxwrap_julia,
        "lib/libcxxwrap_julia.0.9.0.dylib",
        RTLD_GLOBAL,
    )

    JLLWrappers.@init_library_product(
        libcxxwrap_julia_stl,
        "lib/libcxxwrap_julia_stl.dylib",
        RTLD_GLOBAL,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LinearAlgebra

X = randn(40_000,40_000);
XX = X'X;

XX1 = copy(XX)
@time F1 = eigen(XX1)
@show F1.values[end]

XX2 = copy(XX)
@time F2 = eigvals(XX2)
@show F2[end]

XX3 = copy(XX)
@time F3 = LAPACK.syev!('N', 'U', XX3) #!! will change XX3 inplace
@show F3[end]

XX4 = copy(XX)
using Arpack
eigs(XX4,nev=1)[1][1]

#1. no need to Symmetric(XX2), because Julia will detect whether it is Symmetric.
#2. eigmax: first do eigvals(), then return the maximum
#3. both eigen!() and eigvals!() use BLAS function LAPACK.sygvd! if input is symmetric.
#4. eigen!() is slower than eigvals!(). Why?

#5. eigs() in Arpack.jl is faster than eigvals() because eigs() use Arnoldi iterative methods for top eigen values, but eigvals() use LAPACK (QR algorithm). QR algotirhm will return all the eigen values.
#http://hua-zhou.github.io/teaching/biostatm280-2017spring/slides/16-eigsvd/eigsvd.html#Lanczos/Arnoldi-iterative-method-for-top-eigen-pairs

max_stars_repo_name: zhaotianjing/JuliaTips

max_issues_repo_name: zhaotianjing/JuliaTips

max_forks_repo_name: zhaotianjing/JuliaTips


text:
"""
    Krylov solver
"""
mutable struct PetscKSP
    ksp::Ref{KSP}

    function PetscKSP()
        ksp = Ref{KSP}()
        ccall((:KSPCreate, library), PetscErrorCode, (comm_type, Ref{KSP}), MPI.COMM_WORLD, ksp)
        ccall((:KSPSetFromOptions, library), PetscErrorCode, (KSP,), ksp[])
        new(ksp)
    end
end

"""
    Set the matrix to use (A)
"""
function setOperators(ksp::PetscKSP, mat::PetscMat)
    ccall((:KSPSetOperators, library), PetscErrorCode, (KSP, Mat, Mat), ksp.ksp[], mat.mat[], mat.mat[])
end

"""
    Solve the linear system

    Ax=b
"""
function solve!(ksp::PetscKSP, b::PetscVec, x::PetscVec)
    ccall((:KSPSolve, library), PetscErrorCode, (KSP, Vec, Vec), ksp.ksp[], b.vec[], x.vec[])
end


text:
facts("Sparse (semi-)nonnegative Tucker decomposition") do

context("nonnegative decomposition") do
    # An example of nonnegative and semi-nonnegative Tucker decomposition
    tucker_orig = rand_tucker((4, 5, 6), (40, 50, 60), factors_nonneg=true, core_nonneg=true)
    tnsr_orig = compose(tucker_orig)

    tnsr_max = maximum(tnsr_orig)
    map!(x -> x / tnsr_max, tucker_orig.core)
    map!(x -> x / tnsr_max, tnsr_orig)

    tnsr = add_noise(tnsr_orig, 0.6, true)

    # Solve the problem
    @time tucker_spnn = spnntucker(tnsr, size(tucker_orig.core), tol=1E-4, ini_decomp=:hosvd,
                             core_nonneg=true,
                             max_iter=1000, verbose=true, lambdas=fill(0.1, 4))

    # Reporting
    @fact rel_residue(tucker_spnn) --> less_than(0.05)
    info("Relative error of decomposition : $(rel_residue(tucker_spnn))")
end

context("semi-nonnegative decomposition") do
    tucker_orig = rand_tucker((4, 5, 6), (40, 50, 60), factors_nonneg=true, core_nonneg=false)
    tnsr_orig = compose(tucker_orig)

    tnsr_max = maximum(tnsr_orig)
    map!(x -> x / tnsr_max, tucker_orig.core)
    map!(x -> x / tnsr_max, tnsr_orig)

    tnsr = add_noise(tnsr_orig, 0.6, false)

    # Solve the problem
    @time tucker_spnn = spnntucker(tnsr, size(tucker_orig.core), tol=1E-4, ini_decomp=:hosvd,
                             core_nonneg=false,
                             max_iter=1000, verbose=true, lambdas=fill(0.1, 4))

    # Reporting
    @fact rel_residue(tucker_spnn) --> less_than(0.05)
    info("Relative error of decomposition : $(rel_residue(tucker_spnn))")
end

end


text:
import Base.convert

convert(::Type{Vector{Bool}}, string::ASCIIString) = [(ch != '0') for ch in string]

function readbitstream(filename::ASCIIString)
    open(filename, "r") do file
        convert(Vector{Bool}, readall(file))
    end
end

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
using WeightedArrays
using Test

R = Weighted(rand(3,10), rand(10))

@test size(R) == (3,10)
@test lastlength(R) == 10

@test sum(R.weights) ≈ 1 ## normalised in construction

R.weights .+= 1
normalise!(R)
@test sum(R.weights) ≈ 1 ## restored

S = sobol(2,9)

@test size(unique(S + S)) == size(S) ## hcat & unique

# @test size(unique!(clamp!(S .+ 2))) == (2,1) ## broadcast & clamp!
# @test size(unique!(clamp!(wrandn(2,9) .+ 2))) == (2,9)

X = xgrid(2,0:0.1:1)
X = near(X, S, 0.2)
@test size(X) == (2,78) ## near

# T = unique(symm(wgrid(1, 0:0.3:1)) )
# @test weights(T)[1] ≈ 0.25


text:
"""
encode(input)

Program to implement atbash cipher for the given sentence.A full description of the algorithm can be found on [wikipedia](https://en.wikipedia.org/wiki/Atbash)

# Arguments:
- `input` : The sentence needed to rotate

# Examples/Tests
```julia
julia> encode("test")
gvhg

julia> encode("abcdefghijklmnopqrstuvwxyz")
zyxwvutsrqponmlkjihgfedcba

julia> encode("hello")
svool

```

# Algorithm:

```julia

for r in input
    part *= xform(r)
    if length(part) >= 5
      push!(parts, part)
      part = ""
    end
  end
  if part != ""
    push!(parts, part)
  end
  return join(parts, " ")

```

# References:
https://en.wikipedia.org/wiki/Atbash

```

# Contributed by:- [Ihjass Thasbekha](https://github.com/Ihjass)
"""
function encode(input)
  input = replace(lowercase(input), reject_re => "")
  parts = []
  part = ""
  for r in input
    part *= xform(r)
    if length(part) >= 5
      push!(parts, part)
      part = ""
    end
  end
  if !isempty(part)

    push!(parts, part)
  end
  return join(parts, " ")
end
reject_re = r"[^a-z\d]+"
xform(r) = (r >= 'a' && r <= 'z') ? r = 25 - (r - 'a') + 'a' : r

max_stars_repo_name: standardgalactic/Julia-algorithms


text:
export GetHelmholtzOperator,GetHelmholtzOperator,GetHelmholtzShiftOP,getABL,getSommerfeldBC,getHelmholtzFun,getMaximalFrequency

"""
function  GetHelmholtzOperator

Constructor for the acoustic Helmholtz equation

Special Fields:
omega can be real or complex for artificial attenuation.
Artificial attenuation through a complex omega keeps the phase fixed.
gamma - True attenuation parameter, corresponding to the time domain equation:  u_tt + gamma*ut + c^2*L*u = q
"""
function GetHelmholtzOperator(Msh::RegularMesh, mNodal::Array{Float64}, omega::Union{Float64,ComplexF64}, gamma::Array,
									NeumannAtFirstDim::Bool,ABLpad::Array{Int64},ABLamp::Float64,Sommerfeld::Bool)
if gamma == []
	gamma = getABL(Msh.n.+1,NeumannAtFirstDim,ABLpad,ABLamp);
else
	gamma += getABL(Msh.n.+1,NeumannAtFirstDim,ABLpad,ABLamp);
end
H = GetHelmholtzOperator(Msh, mNodal, omega, gamma, NeumannAtFirstDim,Sommerfeld);
return H,gamma
end

function GetHelmholtzOperator(Msh::RegularMesh, mNodal::Array{Float64}, omega::Union{Float64,ComplexF64}, gamma::Array{Float64},
							  NeumannAtFirstDim::Bool,Sommerfeld::Bool)
Lap   = getNodalLaplacianMatrix(Msh);
# this code computes a laplacian the long way, AND stores the gradient on Msh... So we avoid using it.
# Grad  = getNodalGradientMatrix(Msh)
# Lap   = Grad'*Grad

# mass = -((omega.^2).*(mNodal[:]).*(1.0.+1im*gamma[:]));
mass = -(omega.^2).*(mNodal[:]).*(1.0.-1im*gamma[:]./real(omega));

if Sommerfeld
	# println("Adding Sommerfeld");
	somm = getSommerfeldBC(Msh,mNodal,real(omega),NeumannAtFirstDim);
	mass -= somm[:];
end
H = Lap .+ sparse(Diagonal(mass));
return H;
end

function getMaximalFrequency(m::Union{Array{Float64},Array{Float32},Float64},M::RegularMesh)
## m here is in slowness squared
omegamax = (0.1*2*pi)./(maximum(M.h)*sqrt(maximum(m)));
return omegamax;
end

function GetHelmholtzShiftOP(mNodal::Array{Float64}, omega::Float64,shift::Float64)
return sparse(Diagonal(mNodal[:].*(1im*shift*omega^2)));
end

function getHelmholtzFun(ShiftedHelmholtT::SparseMatrixCSC,ShiftMat::SparseMatrixCSC,y::ArrayTypes,numCores::Int64)
function Hfun(x)
		# # here we avoid the storage of the Helmholtz matrix by using the shifted matrix plus the shift.
		Zero = zero(ComplexF64);
		One = one(ComplexF64);
		SpMatMul(One,ShiftedHelmholtT,x,Zero,y,numCores);
		SpMatMul(One,ShiftMat,x,One,y,numCores);
		return y;
end
return Hfun;
end

function getABL(n::Array{Int64},NeumannAtFirstDim::Bool,ABLpad::Array{Int64},ABLamp::Float64)
  pad = ABLpad;
  ntup = tuple(n...);

  if length(n)==2
	# x1 = linspace(-1,1,n[1]);
	# x2 = linspace(0,1,n[2]);
	# X1,X2 = ndgrid(x1,x2);
	# padx1 = ABLpad[1];
	# padx2 = ABLpad[2];
	# gammaxL = (X1 - x1[padx1]).^2;
	# gammaxL[padx1+1:end,:] = 0
	# gammaxR = (X1 - x1[end-padx1+1]).^2
	# gammaxR[1:end-padx1,:] = 0

	# gammax = gammaxL + gammaxR
	# gammax = gammax/maximum(gammax);

	# gammaz = (X2 - x2[end-padx2+1]).^2
	# gammaz[:,1:end-padx2] = 0
	# gammaz = gammaz/maximum(gammaz);

	# gamma = gammax + gammaz
	# gamma *= ABLamp;
	# gamma[gamma.>=ABLamp] = ABLamp;

	gamma = zeros(ntup);
	b_bwd1 = ((pad[1]:-1:1).^2)./pad[1]^2;
	b_bwd2 = ((pad[2]:-1:1).^2)./pad[2]^2;

	b_fwd1 = ((1:pad[1]).^2)./pad[1]^2;
	b_fwd2 = ((1:pad[2]).^2)./pad[2]^2;
	I1 = (n[1] - pad[1] + 1):n[1];
	I2 = (n[2] - pad[2] + 1):n[2];

	if NeumannAtFirstDim==false
		gamma[:,1:pad[2]] += ones(n[1],1)*b_bwd2';
		gamma[1:pad[1],1:pad[2]] -= b_bwd1*b_bwd2';
		gamma[I1,1:pad[2]] -= b_fwd1*b_bwd2';
	end

	gamma[:,I2] +=  ones(n[1],1)*b_fwd2';
	gamma[1:pad[1],:] += b_bwd1*ones(1,n[2]);
	gamma[I1,:] += b_fwd1*ones(1,n[2]);
	gamma[1:pad[1],I2] -= b_bwd1*b_fwd2';
	gamma[I1,I2] -= b_fwd1*b_fwd2';
	gamma *= ABLamp;
	# figure()
	# imshow(gamma'); colorbar()

  else
	x1 = linspace(-1,1,n[1]);
	x2 = linspace(-1,1,n[2]);
	x3 = linspace( 0,1,n[3]);
	X1,X2,X3 = ndgrid(x1,x2,x3);
	padx1 = ABLpad[1];
	padx2 = ABLpad[2];
	padx3 = ABLpad[3];
	gammaL = (X1 - x1[padx1]).^2;
	gammaL[padx1+1:end,:,:] = 0.0
	gammaR = (X1 - x1[end-padx1+1]).^2
	gammaR[1:end-padx1,:,:] = 0.0

	gammat = gammaL + gammaR;
	gammat = gammat/maximum(gammat);
	gamma = copy(gammat);
	gammat[:] = 0.0;

	gammaL = (X2 - x2[padx2]).^2;
	gammaL[:,padx2+1:end,:] = 0.0
	gammaR = (X2 - x2[end-padx2+1]).^2
	gammaR[:,1:end-padx2,:] = 0.0

	gammat = gammaL + gammaR
	gammat = gammat/maximum(gammat);
	gamma += gammat;

	gammat = (X3 - x3[end-padx3+1]).^2
	gammat[:,:,1:end-padx3] = 0.0
	gammat = gammat/maximum(gammat);
	gamma += gammat;
	gamma *= ABLamp;
	gamma[gamma.>=ABLamp] = ABLamp;
  end
  return gamma;
end

function getSommerfeldBC(Msh::RegularMesh,mNodal::Array{Float64}, omega::Float64,NeumannOnTop::Bool)
ntup = tuple((Msh.n .+ 1)...);
Somm = zeros(ComplexF64,ntup);
mNodal = reshape(mNodal,ntup);
h = Msh.h;

if Msh.dim==2
	if !NeumannOnTop
		Somm[:,1] += -1im*omega*(BC()/h[2]).*sqrt.(mNodal[:,1]);
	end
	Somm[:,end] += (-1im*omega*(BC()/h[2])).*sqrt.(mNodal[1:end,end]);
	Somm[end,:] += (-1im*omega*(BC()/h[1])).*sqrt.(mNodal[end,:]);
	Somm[1,:] += (-1im*omega*(BC()/h[1])).*sqrt.(mNodal[1,:]);
else
	if !NeumannOnTop
		Somm[:,:,1] += -1im*omega*(BC()/h[3]).*sqrt.(mNodal[:,:,1]);
	end
	Somm[:,:,end] += -1im*omega*(BC()/h[3]).*sqrt.(mNodal[:,:,end]);
	Somm[:,1,:] += -1im*omega*(BC()/h[2]).*sqrt.(mNodal[:,1,:]);
	Somm[:,end,:] += -1im*omega*(BC()/h[2]).*sqrt.(mNodal[:,end,:]);
	Somm[1,:,:] += -1im*omega*(BC()/h[1]).*sqrt.(mNodal[1,:,:]);
	Somm[end,:,:] += -1im*omega*(BC()/h[1]).*sqrt.(mNodal[end,:,:]);
end
return Somm;
end


text:
@define_scenario DayAheadScenario = begin
    ρ::PriceCurve{24,Float64}
    Q̃::Inflows{typeof(Skellefteälven),Float64}

    @expectation begin
        ρ = mean([s.ρ for s in scenarios])
        Q̃ = mean([s.Q̃ for s in scenarios])
        return DayAheadScenario(ρ, Q̃, probability = 1.0)
    end
end

function StochasticPrograms.expected(ξ₁::DayAheadScenario, ξ₂::DayAheadScenario)
    ρ = mean([ξ₁.ρ, ξ₂.ρ])
    Q̃ = mean([ξ₁.Q̃, ξ₂.Q̃])
    expected = DayAheadScenario(ρ, Q̃)
    return ExpectedScenario(expected)
end

function penalty(scenario::DayAheadScenario, t)
    if 8 <= t <= 12 || 17 <= t <= 24
        return 1.15*scenario.ρ[t]
    else
        return 1.1*scenario.ρ[t]
    end
end
function penalty(scenario::DayAheadScenario, α, t)
    α == 0.0 && return penalty(scenario, t)
    return α*scenario.ρ[t]
end
function reward(scenario::DayAheadScenario, t)
    if 8 <= t <= 12 || 17 <= t <= 24
        return 0.85*scenario.ρ[t]
    else
        return 0.9*scenario.ρ[t]
    end
end
function reward(scenario::DayAheadScenario, β, t)
    β == 0.0 && return reward(scenario, t)
    return β*scenario.ρ[t]
end

@sampler RecurrentDayAheadSampler = begin
    date::Date
    plants::PlantCollection
    price_forecaster::Forecaster
    flow_forecaster::Forecaster

    @sample DayAheadScenario begin
        price_curve = forecast(sampler.price_forecaster, month(sampler.date))
        while !all(price_curve .>= 0)
            price_curve = forecast(sampler.price_forecaster, month(sampler.date))
        end
        flows = forecast(sampler.flow_forecaster, week(sampler.date))
        while !all(flows[:, 2] .>= 0)
            flows = forecast(sampler.flow_forecaster, week(sampler.date))
        end
        return DayAheadScenario(PriceCurve(price_curve), Inflows(sampler.plants, flows[:, 2]))
    end
end

function generate_bidlevels(sampler::AbstractSampler{DayAheadScenario})
    price_curves = [sampler().ρ for i = 1:1000]
    bidlevels = zeros(5,24)
    for i = 1:24
        hourly_prices = [p[i] for p in price_curves]
        μ = mean(hourly_prices)
        σ = std(hourly_prices)
        bidlevels[1,i] = μ - 2*σ
        bidlevels[2,i] = μ - σ
        bidlevels[3,i] = μ
        bidlevels[4,i] = μ + σ
        bidlevels[5,i] = μ + 2*σ
    end
    return bidlevels
end


text:
using Sockets # for version 1.0
println("Echo server on port 12321")
try
    server = listen(12321)
    instance = 0
    while true
        sock = accept(server)
        instance += 1
        socklabel = "$(getsockname(sock)) number $instance"
        @async begin
            println("Server connected to socket $socklabel")
            write(sock, "Connected to echo server.\r\n")
            while isopen(sock)
                str = readline(sock)
                write(sock,"$str\r\n")
                println("Echoed $str to socket $socklabel")
            end
            println("Closed socket $socklabel")
        end
    end
catch y
    println("Caught exception: $y")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text: ## This file will be automatically completed during the package build process


text:
module ShapeReconstructionPaLS
using jInv
using ParamLevelSet
include("Utils/Utils.jl")
include("DipTransform/DipTransform.jl")
include("ShapeFromSilhouette/ShapeFromSilhouette.jl")
include("Direct/Direct.jl")
include("PointCloud/PointCloud.jl")
end # module

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
using Hexagons
using Test
using Random: seed!
using Statistics: mean

# Test a few identities for the hexagon containing the point x, y
function run_point_test(x, y, hextype)
    hex = hex_containing(hextype, x, y)
    verts = vertices(hex)
    # x,y should be in the bounding box of the hexagon vertices
    @test minimum(v[1] for v in verts) <= x <= maximum(v[1] for v in verts)
    @test minimum(v[2] for v in verts) <= y <= maximum(v[2] for v in verts)
    # the center of the hexagon should be near its vertex mean
    mean_vert_x = mean(v[1] for v in verts)
    mean_vert_y = mean(v[2] for v in verts)
    hex_center = center(hex)
    @test isapprox(hex_center[1], mean_vert_x; atol = 1e-6)
    @test isapprox(hex_center[2], mean_vert_y; atol = 1e-6)
    # a string of type conversions should recover the original hex
    hex_axial = to_coord_form(CoordAxial, hex)
    hex_cubic = to_coord_form(CoordCubic, hex_axial)
    @test hex == hex_cubic
end

test_points = [
    (0, 0),
    (0, 1),
    (1, 0),
    (1, 1),
    (1, -1),
    (0, 47),
    (-4.7, 0),
    (1.234, 5.678),
    (1e6, 2e6),
]
# bunch of random test points
seed!(1234)
for _ in 1:1000
    push!(test_points, (rand() * 100, rand() * 100))
end
for point in test_points
    run_point_test(point..., HexPointyTop)
    run_point_test(point..., HexFlatTop)
end

# Test some hash and equality identities
ax0_int = CoordAxial(0, 0)
ax0_float = CoordAxial(0.0, 0.0)
ax1_int = CoordAxial(1, 1)
@test ax0_float == ax0_int
@test hash(ax0_float) == hash(ax0_int)
@test ax1_int != ax0_int
@test hash(ax1_int) != hash(ax0_int)
for T in (HexFlatTop, HexPointyTop)
    @test T(ax0_float) == T(ax0_int)
    @test hash(T(ax0_float)) == hash(T(ax0_int))
end
@test HexFlatTop(ax0_int) != HexPointyTop(ax0_int)


text:
"""
    struct Name
        particle::String
        last::String
        junior::String
        first::String
        middle::String
    end
Store a name with no ambiguities. For additional parts to define a name, please fill an issue or make a pull request.
"""
# Name type
struct Name
    particle::String
    last::String
    junior::String
    first::String
    middle::String
end
Names = Vector{Name}

"""
    Name(str::String)
Decompose without ambiguities a name as `particle` (optional) `last`, `junior` (optional), `first` `middle` (optional) based on BibTeX possible input. As for BibTeX, the decomposition of a name in the form of `first` `last` is also possible, but ambiguities can occur.
"""
function Name(str::String)
    subnames = split(str, r"[\n\r ]*,[\n\r ]*")

    # subnames containers
    first = ""
    middle = ""
    particle = ""
    last = ""
    junior = ""

    # mark for string parsing
    mark_in = 1
    mark_out = 0

    # BibTeX form 1: First Second von Last
    if length(subnames) == 1
        aux = split(subnames[1], r"[\n\r ]+")
        mark_out = length(aux) - 1
        last = aux[end]
        if length(aux) > 1 && isuppercase(aux[1][1])
            first = aux[1]
            mark_in += 1
            for s in aux[2:end-1]
                mark_in += 1
                if islowercase(s[1])
                    break;
                end
                middle *= " $s"
            end
            for s in reverse(aux[mark_in:mark_out])
                if islowercase(s[1])
                    break;
                end
                mark_out -= 1
                last = "$s " * last
            end
            for s in aux[mark_in:mark_out]
                particle *= " $s"
            end
        end
    end
    # BibTeX form 2: von Last, First Second
    if length(subnames) == 2
        aux = split(subnames[1], r"[\n\r ]+") # von Last
        mark_out = length(aux) - 1
        last = string(aux[end])
        for s in reverse(aux[1:mark_out])
            if islowercase(s[1])
                break;
            end
            mark_out -= 1
            last = "$s " * last
        end
        for s in aux[1:mark_out]
            particle *= " $s"
        end
        aux = split(subnames[2], r"[\n\r ]+")
        first = aux[1]
        if length(aux) > 1
            for s in aux[2:end]
                middle *= " $s"
            end
        end
    end
    if length(subnames) == 3
        aux = split(subnames[1], r"[\n\r ]+") # von Last
        mark_out = length(aux) - 1
        last = aux[end]
        for s in reverse(aux[1:mark_out])
            if islowercase(s[1])
                break;
            end
            mark_out -= 1
            last = "$s " * last
        end
        for s in aux[1:mark_out]
            particle *= " $s"
        end
        junior = subnames[2]
        aux =split(subnames[3], r"[\n\r ]+")
        first = aux[1]
        if length(aux) > 1
            for s in aux[2:end]
                middle *= " $s"
            end
        end
    end
    return Name(particle, last, junior, first, middle)
end

"""
    names(str::String)
Decompose into parts a list of names in BibTeX compatible format. That is names sparated by `and`.
"""
function names(str::String)
    aux = split(str, r"[\n\r ]and[\n\r ]")
    return map(x -> Name(String(x)), aux)
end

"""
    struct Access
        doi::String
        howpublished::String
        url::String
    end
Store the online access of an entry as a String. Handles the fields `doi` and `url` and the `arXiV` entries. For additional fields or entries, please fill an issue or make a pull request.
"""
struct Access
    doi::String
    howpublished::String
    url::String
end

"""
    Access(fields::Fields)
Construct the online access information based on the entry fields.
"""
function Access(fields::Fields)
    doi = get_delete!(fields, "doi")
    howpublished = get_delete!(fields, "howpublished")
    url = fields["type"] == "eprint" ? arxive_url(fields) : get_delete!(fields, "url")
    Access(doi, howpublished, url)
end

"""
    struct Date
        day::String
        month::String
        year::String
    end
Store the date information as `day`, `month`, and `year`.
"""
struct Date
    day::String
    month::String
    year::String
end

"""
    Date(fields::Fields)
Construct the date information based on the entry fields.
"""
function Date(fields::Fields)
    day = get_delete!(fields, "day")
    month = get_delete!(fields, "month")
    year = get_delete!(fields, "year")
    Date(day, month, year)
end

"""
    struct Eprint
        archive_prefix::String
        eprint::String
        primary_class::String
    end
Store the information related to arXiv eprint format.
"""
struct Eprint
    archive_prefix::String
    eprint::String
    primary_class::String
end

"""
    Eprint(fields::Fields)
Construct the eprint arXiv information based on the entry fields. Handle old and current arXiv format.
"""
function Eprint(fields::Fields)
    archive_prefix = get_delete!(fields, "archivePrefix")
    eprint = get_delete!(fields, "eprint")
    primary_class = get_delete!(fields, "primaryClass")
    Eprint(archive_prefix, eprint, primary_class)
end

"""
    struct In
        address::String
        chapter::String
        edition::String
        institution::String
        journal::String
        number::String
        organization::String
        pages::String
        publisher::String
        school::String
        series::String
        volume::String
    end
Store all the information related to how an entry was published.
"""
struct In
    address::String
    chapter::String
    edition::String
    institution::String
    journal::String
    number::String
    organization::String
    pages::String
    publisher::String
    school::String
    series::String
    volume::String
end

"""
    In(fields::Fields)
Construct the information of how an entry was published based on its fields
"""
function In(fields::Fields)
    address = get_delete!(fields, "address")
    chapter = get_delete!(fields, "chapter")
    edition = get_delete!(fields, "edition")
    institution = get_delete!(fields, "institution")
    journal = get_delete!(fields, "journal")
    number = get_delete!(fields, "number")
    organization = get_delete!(fields, "organization")
    pages = get_delete!(fields, "pages")
    publisher = get_delete!(fields, "publisher")
    school = get_delete!(fields, "school")
    series = get_delete!(fields, "series")
    volume = get_delete!(fields, "volume")
    In(address, chapter, edition, institution, journal, number, organization, pages, publisher, school, series, volume)
end

"""
    struct Entry <: AbstractEntry
        access::Access
        authors::Names
        booktitle::String
        date::Date
        editors::Names
        eprint::Eprint
        id::String
        in::In
        fields::Dict{String,String}
        title::String
        type::String
    end
Generic Entry type. If some construction rules are required, it should be done beforehand. Check `bibtex.jl` as the example of rules implementation for BibTeX format.
"""
struct Entry <: AbstractEntry
    access::Access
    authors::Names
    booktitle::String
    date::Date
    editors::Names
    eprint::Eprint
    id::String
    in::In
    fields::Dict{String,String}
    title::String
    type::String
end

"""
    Entry(id::String, fields::Fields)
Construct an entry with a unique id and a list of `Fields`.
"""
function Entry(id::String, fields::Fields)
    # foreach((name, field) -> fields[name] = erase_spaces(field), fields)
    access = Access(fields)
    authors = names(get_delete!(fields, "author"))
    booktitle = get_delete!(fields, "booktitle")
    date = Date(fields)
    editors = names(get_delete!(fields, "editor"))
    eprint = Eprint(fields)
    in_ = In(fields)
    title = get_delete!(fields, "title")
    type = get_delete!(fields, "type")
    return Entry(access, authors, booktitle, date, editors, eprint, id, in_, fields, title, type)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.2'
#       jupytext_version: 1.2.1
#   kernelspec:
#     display_name: Julia 1.6.0-beta1 depwarn -O3
#     language: julia
#     name: julia-1.6-depwarn-o3
# ---

# %% [markdown]
# # Ising2D
#
# * Copyright (c) 2020 Gen Kuroki
# * License: MIT
# * Repository: https://github.com/genkuroki/Ising2D.jl

# %% [markdown] {"toc": true}
#
Table of Contents

#

1  Benchmarks

2  Examples

3  Documents

# %%
if isfile("Project.toml")
    using Pkg
    Pkg.activate(".")
    using Revise
end

# %%
using Ising2D, Plots, Random

# %% [markdown]
# ## Benchmarks

# %%
using Ising2D, Random
Random.seed!(464937337564)
print("IfElse() on v$VERSION:")
@time s = ising2d!(rand_ising2d(), β_ising2d, 10^5; algorithm=IfElse())
plot_ising2d(s)

# %%
using Ising2D, Random
Random.seed!(464937337564)
print("MultiFor() on v$VERSION:")
@time s = ising2d!(rand_ising2d(), β_ising2d, 10^5; algorithm=MultiFor())
plot_ising2d(s)

# %%
using Ising2D, BenchmarkTools, Random
Random.seed!(464937337564)
s = rand_ising2d()
print("IfElse() on v$VERSION:")
@btime ising2d!($s; algorithm=IfElse());

# %%
using Ising2D, BenchmarkTools, Random
Random.seed!(464937337564)
s = rand_ising2d()
print("MultiFor() on v$VERSION:")
@btime ising2d!($s; algorithm=MultiFor());

# %% [markdown]
# ## Examples

# %%
Random.seed!(464937337564);

# %%
s = rand_ising2d(200)
P0 = plot_ising2d(s)

# %%
ising2d!(s, β_ising2d, 500)
P1 = plot_ising2d(s)

# %%
png(P0, "images/s0.png")
png(P1, "images/s1.png")

# %%
s = rand_ising2d(200)
gif_ising2d(s=s, k=1.0, nwarmups=0, nskips=1, niters=500, fps=15,
    gifname="images/ising2d.gif")

# %%
s₀ = rand_ising2d()
@show energy_density_ising2d(s₀)
@show magnetization_ising2d(s₀)
s = copy(s₀)
S = mcmc_ising2d!(s=s, k=1.0, nwarmups=1000, nskips=100, niters=5000)
E = energy_density_ising2d.(S)
M = magnetization_ising2d.(S)
@show mean(E)
@show mean(M)

histogram_mcmc_ising2d(S)

# %%
plot_mcmc_ising2d(S)

# %%
plot_mcmc_ising2d(S; t=4000)

# %%
gif_mcmc_ising2d(S[10:10:end]; gifname="images/ising2d_mcmc.gif")

# %% [markdown]
# ## Documents

# %%
?Ising2D

# %%
?Ising2D.β_ising2d

# %%
?Ising2D.rand_ising2d

# %%
?Ising2D.ones_ising2d

# %%
?Ising2D.ising2d!

# %%
?Ising2D.plot_ising2d

# %%
?Ising2D.gif_ising2d

# %%
?Ising2D.mcmc_ising2d!

# %%
?Ising2D.energy_density_ising2d

# %%
?Ising2D.magnetization_ising2d

# %%
?Ising2D.plot_mcmc_ising2d

# %%
?Ising2D.histogram_mcmc_ising2d

# %%
?Ising2D.gif_mcmc_ising2d

# %%


text:
#=
dev.jl is a draft file that I use during the development of a new functionality if I need to call some function, etc. and see how it works.
I.e. at the beginning of the development of a new functionality, I disable all tests except dev.jl,
after the overall architecture emerges, I disable dev.jl and start writing full-fledged tests
=#
using Test
using Dash
using JSON3
using HTTP
@testset "dev" begin
    r = Dash.load_meta("dash_core_components")
    c = Dash.generate_embeded_components()
    println(c)
end


text:
"""
    struct Sub

Represents a subreddit on Reddit.
"""
struct Subreddit
    name::String
end

"""
    about(sub::Subreddit)

Get information about a specified subreddit using default credentials.
"""
function about(sub::Subreddit)
    about(sub, default())
end

"""
    about(sub::Subreddit, a::Authorized)

Get information about a specified subreddit.
"""
function about(sub::Subreddit, a::Authorized)
    JSON.parse(get("/r/$(sub.name)/about", a))["data"]
end

"""
    subscribers(sub::AbstractString)

Get total number of subscribers for a subreddit by name, using default
credentials.
"""
function subscribers(sub::AbstractString)
    subscribers(sub, default())
end

"""
    subscribers(sub::Subreddit)

Get total number of subscribers for a subreddit with Subreddit type.
"""
function subscribers(sub::Subreddit)
    subscribers(sub.name, default())
end

"""
    subscribers(sub::Subreddit, a::Authorized)

Get total number of subscribers for a subreddit with Subreddit type.
"""
function subscribers(sub::Subreddit, a::Authorized)
    subscribers(sub.name, a)
end

"""
    subscribers(sub::AbstractString, a::Authorized)

Get total number of subscribers for a subreddit by name.
"""
function subscribers(sub::AbstractString, a::Authorized)
    about(Subreddit(sub), a)["subscribers"]
end


text:
module JuliaBerry
__precompile__(false)
using PiGPIO

export Pin, InputPin, OuputPin,
	LED, Motor, Buzzer,
	forward, backward, speed, stop, on, off, is_on, is_off, frequency

# package code goes here

abstract type Thing
end

struct Pin <: Thing
    pin::Int

    function Pin(pin, mode)
        PiGPIO.set_mode(_pi[], pin, mode)
        new(pin)
    end
end

InputPin(pin) = Pin(pin, PiGPIO.INPUT)
OutputPin(pin) = Pin(pin, PiGPIO.OUTPUT)

struct LED <: Thing
    pin::Int

    function LED(pin)
        PiGPIO.set_mode(_pi[], pin, PiGPIO.OUTPUT)
        new(pin)
    end
end

struct Buzzer <: Thing
	pin::Int

	function Buzzer(pin, freq)
		PiGPIO.set_mode(_pi[], pin, PiGPIO.OUTPUT)
		PiGPIO.set_PWM_frequency(_pi[], pin, freq)
		new(pin)
	end
end

Buzzer(pin) = Buzzer(pin, 150)
function frequency(b::Buzzer, freq::Int)
	PiGPIO.set_PWM_frequency(_pi[], b.pin, freq)
end

on(b::Buzzer) = PiGPIO.set_PWM_dutycycle(_pi[], b.pin, 100)
off(b::Buzzer) = PiGPIO.set_PWM_dutycycle(_pi[], b.pin, 0)

function on(x::Thing)
    PiGPIO.write(_pi[], x.pin, 1)
end

function off(x::Thing)
    PiGPIO.write(_pi[], x.pin, 0)
end

on(x::Array{T}) where T<:Thing = on.(x)
off(x::Array{T}) where T<:Thing = off.(x)

is_off(x::Thing) = PiGPIO.read(_pi[], x.pin) == 0
is_on(x::Thing) = PiGPIO.read(_pi[], x.pin) == 1

struct Motor <: Thing
    fw_pin::Int
    bw_pin::Int

    function Motor(fw_pin, bw_pin)
        PiGPIO.set_mode(_pi[], fw_pin, PiGPIO.OUTPUT)
        PiGPIO.set_mode(_pi[], bw_pin, PiGPIO.OUTPUT)
        PiGPIO.set_PWM_range(_pi[], fw_pin, 100)
        PiGPIO.set_PWM_range(_pi[], bw_pin, 100)
        PiGPIO.set_PWM_dutycycle(_pi[], fw_pin, 0)
        PiGPIO.set_PWM_dutycycle(_pi[], bw_pin, 0)
	new(fw_pin, bw_pin)
    end
end

function speed(x::Motor, speed::Int=100)
    if speed>0
        forward(x, speed)
    elseif speed<0
        backward(x, abs(speed))
    end
end

function stop(x::Motor)
    PiGPIO.set_PWM_dutycycle(_pi[], x.fw_pin, 0)
    PiGPIO.set_PWM_dutycycle(_pi[], x.bw_pin, 0)
end

function forward(x::Motor, speed::Int)
    PiGPIO.set_PWM_dutycycle(_pi[], x.fw_pin, speed)
    PiGPIO.set_PWM_dutycycle(_pi[], x.bw_pin, 0)
end

function backward(x::Motor, speed::Int)
    PiGPIO.set_PWM_dutycycle(_pi[], x.fw_pin, 0)
    PiGPIO.set_PWM_dutycycle(_pi[], x.bw_pin, speed)
end

global const _pi = Ref{Pi}()
function __init__()
    _pi[] = Pi()
end
#TODO remove
__init__()

setpi(p::Pi) = _pi[]=p

include("explorerhat.jl")
end # module


text:
abstract type CPAAlphaModel <: AlphaModel end

struct CPAAlphaParam <: EoSParam
    c1::SingleParam{Float64}
end

@newmodelsimple CPAAlpha CPAAlphaModel CPAAlphaParam

export CPAAlpha
function CPAAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["SAFT/CPA/CPA_like.csv"]; userlocations=userlocations, ignore_missing_singleparams=["Mw"], verbose=verbose)
    c1 = params["c1"]
    packagedparams = CPAAlphaParam(c1)
    model = CPAAlpha(packagedparams, verbose=verbose)
    return model
end

function α_function(model::CubicModel,V,T,z,alpha_model::CPAAlphaModel)
    Tc = model.params.Tc.values
    Tr = @. T/Tc
    c1  = alpha_model.params.c1.values
    α  = @. (1+c1*(1-√(Tr)))^2
    return α
end


text:
@testset "4.1.1.1 (a+b sin)^n" begin
    (a, b, c, d, n, x, ) = @variables a b c d n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(a+b*sin(c+d*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+a*sin(c+d*x))^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(a+a*sin(c+d*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+a*sin(c+d*x))^(n/2)=#

    @test_int [(a + a*sin(c + d*x))^(7/2), x, 4, -((256*a^4*cos(c + d*x))/(35*d*sqrt(a + a*sin(c + d*x)))) - (64*a^3*cos(c + d*x)*sqrt(a + a*sin(c + d*x)))/(35*d) - (24*a^2*cos(c + d*x)*(a + a*sin(c + d*x))^(3/2))/(35*d) - (2*a*cos(c + d*x)*(a + a*sin(c + d*x))^(5/2))/(7*d)]
    @test_int [(a + a*sin(c + d*x))^(5/2), x, 3, -((64*a^3*cos(c + d*x))/(15*d*sqrt(a + a*sin(c + d*x)))) - (16*a^2*cos(c + d*x)*sqrt(a + a*sin(c + d*x)))/(15*d) - (2*a*cos(c + d*x)*(a + a*sin(c + d*x))^(3/2))/(5*d)]
    @test_int [(a + a*sin(c + d*x))^(3/2), x, 2, -((8*a^2*cos(c + d*x))/(3*d*sqrt(a + a*sin(c + d*x)))) - (2*a*cos(c + d*x)*sqrt(a + a*sin(c + d*x)))/(3*d)]
    @test_int [(a + a*sin(c + d*x))^(1/2), x, 1, -((2*a*cos(c + d*x))/(d*sqrt(a + a*sin(c + d*x))))]
    @test_int [1/(a + a*sin(c + d*x))^(1/2), x, 2, -((sqrt(2)*atanh((sqrt(a)*cos(c + d*x))/(sqrt(2)*sqrt(a + a*sin(c + d*x)))))/(sqrt(a)*d))]
    @test_int [1/(a + a*sin(c + d*x))^(3/2), x, 3, -(atanh((sqrt(a)*cos(c + d*x))/(sqrt(2)*sqrt(a + a*sin(c + d*x))))/(2*sqrt(2)*a^(3/2)*d)) - cos(c + d*x)/(2*d*(a + a*sin(c + d*x))^(3/2))]
    @test_int [1/(a + a*sin(c + d*x))^(5/2), x, 4, -((3*atanh((sqrt(a)*cos(c + d*x))/(sqrt(2)*sqrt(a + a*sin(c + d*x)))))/(16*sqrt(2)*a^(5/2)*d)) - cos(c + d*x)/(4*d*(a + a*sin(c + d*x))^(5/2)) - (3*cos(c + d*x))/(16*a*d*(a + a*sin(c + d*x))^(3/2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+a*sin(c+d*x))^(n/3)=#

    @test_int [(a + a*sin(c + d*x))^(4/3), x, 2, -((2*2^(5/6)*a*cos(c + d*x)*HypergeometricFunctions._₂F₁(-(5/6), 1/2, 3/2, (1/2)*(1 - sin(c + d*x)))*(a + a*sin(c + d*x))^(1/3))/(d*(1 + sin(c + d*x))^(5/6)))]
    @test_int [(a + a*sin(c + d*x))^(2/3), x, 2, -((2*2^(1/6)*cos(c + d*x)*HypergeometricFunctions._₂F₁(-(1/6), 1/2, 3/2, (1/2)*(1 - sin(c + d*x)))*(a + a*sin(c + d*x))^(2/3))/(d*(1 + sin(c + d*x))^(7/6)))]
    @test_int [(a + a*sin(c + d*x))^(1/3), x, 2, -((2^(5/6)*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/6, 1/2, 3/2, (1/2)*(1 - sin(c + d*x)))*(a + a*sin(c + d*x))^(1/3))/(d*(1 + sin(c + d*x))^(5/6)))]
    @test_int [1/(a + a*sin(c + d*x))^(1/3), x, 2, -((2^(1/6)*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 5/6, 3/2, (1/2)*(1 - sin(c + d*x))))/(d*(1 + sin(c + d*x))^(1/6)*(a + a*sin(c + d*x))^(1/3)))]
    @test_int [1/(a + a*sin(c + d*x))^(2/3), x, 2, -((cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 7/6, 3/2, (1/2)*(1 - sin(c + d*x)))*(1 + sin(c + d*x))^(1/6))/(2^(1/6)*d*(a + a*sin(c + d*x))^(2/3)))]
    @test_int [1/(a + a*sin(c + d*x))^(4/3), x, 2, -((cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 11/6, 3/2, (1/2)*(1 - sin(c + d*x))))/(2^(5/6)*a*d*(1 + sin(c + d*x))^(1/6)*(a + a*sin(c + d*x))^(1/3)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+a*sin(c+d*x))^n*with*n*symbolic=#

    @test_int [(a + a*sin(c + d*x))^n, x, 2, -((2^(1/2 + n)*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - n, 3/2, (1/2)*(1 - sin(c + d*x)))*(1 + sin(c + d*x))^(-(1/2) - n)*(a + a*sin(c + d*x))^n)/d)]
    @test_int [(a - a*sin(c + d*x))^n, x, 2, (2^(1/2 + n)*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - n, 3/2, (1/2)*(1 + sin(c + d*x)))*(1 - sin(c + d*x))^(-(1/2) - n)*(a - a*sin(c + d*x))^n)/d]

    @test_int [(2 + 2*sin(c + d*x))^n, x, 1, -((2^(1/2 + 2*n)*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - n, 3/2, (1/2)*(1 - sin(c + d*x))))/(d*sqrt(1 + sin(c + d*x))))]
    @test_int [(2 - 2*sin(c + d*x))^n, x, 1, (2^(1/2 + 2*n)*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 1/2 - n, 3/2, (1/2)*(1 + sin(c + d*x))))/(d*sqrt(1 - sin(c + d*x)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(c+d*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(c+d*x))^n=#

    @test_int [1/(5 + 3*sin(c + d*x)), x, 1, x/4 + atan(cos(c + d*x)/(3 + sin(c + d*x)))/(2*d)]
    @test_int [1/(5 + 3*sin(c + d*x))^2, x, 3, (5*x)/64 + (5*atan(cos(c + d*x)/(3 + sin(c + d*x))))/(32*d) + (3*cos(c + d*x))/(16*d*(5 + 3*sin(c + d*x)))]
    @test_int [1/(5 + 3*sin(c + d*x))^3, x, 4, (59*x)/2048 + (59*atan(cos(c + d*x)/(3 + sin(c + d*x))))/(1024*d) + (3*cos(c + d*x))/(32*d*(5 + 3*sin(c + d*x))^2) + (45*cos(c + d*x))/(512*d*(5 + 3*sin(c + d*x)))]
    @test_int [1/(5 + 3*sin(c + d*x))^4, x, 5, (385*x)/32768 + (385*atan(cos(c + d*x)/(3 + sin(c + d*x))))/(16384*d) + cos(c + d*x)/(16*d*(5 + 3*sin(c + d*x))^3) + (25*cos(c + d*x))/(512*d*(5 + 3*sin(c + d*x))^2) + (311*cos(c + d*x))/(8192*d*(5 + 3*sin(c + d*x)))]

    @test_int [1/(5 - 3*sin(c + d*x)), x, 1, x/4 - atan(cos(c + d*x)/(3 - sin(c + d*x)))/(2*d)]
    @test_int [1/(5 - 3*sin(c + d*x))^2, x, 3, (5*x)/64 - (5*atan(cos(c + d*x)/(3 - sin(c + d*x))))/(32*d) - (3*cos(c + d*x))/(16*d*(5 - 3*sin(c + d*x)))]
    @test_int [1/(5 - 3*sin(c + d*x))^3, x, 4, (59*x)/2048 - (59*atan(cos(c + d*x)/(3 - sin(c + d*x))))/(1024*d) - (3*cos(c + d*x))/(32*d*(5 - 3*sin(c + d*x))^2) - (45*cos(c + d*x))/(512*d*(5 - 3*sin(c + d*x)))]
    @test_int [1/(5 - 3*sin(c + d*x))^4, x, 5, (385*x)/32768 - (385*atan(cos(c + d*x)/(3 - sin(c + d*x))))/(16384*d) - cos(c + d*x)/(16*d*(5 - 3*sin(c + d*x))^3) - (25*cos(c + d*x))/(512*d*(5 - 3*sin(c + d*x))^2) - (311*cos(c + d*x))/(8192*d*(5 - 3*sin(c + d*x)))]

    @test_int [1/(-5 + 3*sin(c + d*x)), x, 1, -(x/4) + atan(cos(c + d*x)/(3 - sin(c + d*x)))/(2*d)]
    @test_int [1/(-5 + 3*sin(c + d*x))^2, x, 3, (5*x)/64 - (5*atan(cos(c + d*x)/(3 - sin(c + d*x))))/(32*d) - (3*cos(c + d*x))/(16*d*(5 - 3*sin(c + d*x)))]
    @test_int [1/(-5 + 3*sin(c + d*x))^3, x, 4, -((59*x)/2048) + (59*atan(cos(c + d*x)/(3 - sin(c + d*x))))/(1024*d) + (3*cos(c + d*x))/(32*d*(5 - 3*sin(c + d*x))^2) + (45*cos(c + d*x))/(512*d*(5 - 3*sin(c + d*x)))]
    @test_int [1/(-5 + 3*sin(c + d*x))^4, x, 5, (385*x)/32768 - (385*atan(cos(c + d*x)/(3 - sin(c + d*x))))/(16384*d) - cos(c + d*x)/(16*d*(5 - 3*sin(c + d*x))^3) - (25*cos(c + d*x))/(512*d*(5 - 3*sin(c + d*x))^2) - (311*cos(c + d*x))/(8192*d*(5 - 3*sin(c + d*x)))]

    @test_int [1/(-5 - 3*sin(c + d*x)), x, 1, -(x/4) - atan(cos(c + d*x)/(3 + sin(c + d*x)))/(2*d)]
    @test_int [1/(-5 - 3*sin(c + d*x))^2, x, 3, (5*x)/64 + (5*atan(cos(c + d*x)/(3 + sin(c + d*x))))/(32*d) + (3*cos(c + d*x))/(16*d*(5 + 3*sin(c + d*x)))]
    @test_int [1/(-5 - 3*sin(c + d*x))^3, x, 4, -((59*x)/2048) - (59*atan(cos(c + d*x)/(3 + sin(c + d*x))))/(1024*d) - (3*cos(c + d*x))/(32*d*(5 + 3*sin(c + d*x))^2) - (45*cos(c + d*x))/(512*d*(5 + 3*sin(c + d*x)))]
    @test_int [1/(-5 - 3*sin(c + d*x))^4, x, 5, (385*x)/32768 + (385*atan(cos(c + d*x)/(3 + sin(c + d*x))))/(16384*d) + cos(c + d*x)/(16*d*(5 + 3*sin(c + d*x))^3) + (25*cos(c + d*x))/(512*d*(5 + 3*sin(c + d*x))^2) + (311*cos(c + d*x))/(8192*d*(5 + 3*sin(c + d*x)))]

    @test_int [1/(3 + 5*sin(c + d*x)), x, 4, -(log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x)))/(4*d)) + log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(3 + 5*sin(c + d*x))^2, x, 6, (3*log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(64*d) - (3*log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x))))/(64*d) - (5*cos(c + d*x))/(16*d*(3 + 5*sin(c + d*x)))]
    @test_int [1/(3 + 5*sin(c + d*x))^3, x, 7, -((43*log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(2048*d)) + (43*log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x))))/(2048*d) - (5*cos(c + d*x))/(32*d*(3 + 5*sin(c + d*x))^2) + (45*cos(c + d*x))/(512*d*(3 + 5*sin(c + d*x)))]
    @test_int [1/(3 + 5*sin(c + d*x))^4, x, 8, (279*log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(32768*d) - (279*log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x))))/(32768*d) - (5*cos(c + d*x))/(48*d*(3 + 5*sin(c + d*x))^3) + (25*cos(c + d*x))/(512*d*(3 + 5*sin(c + d*x))^2) - (995*cos(c + d*x))/(24576*d*(3 + 5*sin(c + d*x)))]

    @test_int [1/(3 - 5*sin(c + d*x)), x, 4, -(log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x)))/(4*d)) + log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(3 - 5*sin(c + d*x))^2, x, 6, (3*log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x))))/(64*d) - (3*log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(64*d) + (5*cos(c + d*x))/(16*d*(3 - 5*sin(c + d*x)))]
    @test_int [1/(3 - 5*sin(c + d*x))^3, x, 7, -((43*log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x))))/(2048*d)) + (43*log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(2048*d) + (5*cos(c + d*x))/(32*d*(3 - 5*sin(c + d*x))^2) - (45*cos(c + d*x))/(512*d*(3 - 5*sin(c + d*x)))]
    @test_int [1/(3 - 5*sin(c + d*x))^4, x, 8, (279*log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x))))/(32768*d) - (279*log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(32768*d) + (5*cos(c + d*x))/(48*d*(3 - 5*sin(c + d*x))^3) - (25*cos(c + d*x))/(512*d*(3 - 5*sin(c + d*x))^2) + (995*cos(c + d*x))/(24576*d*(3 - 5*sin(c + d*x)))]

    @test_int [1/(-3 + 5*sin(c + d*x)), x, 4, log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x)))/(4*d) - log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(-3 + 5*sin(c + d*x))^2, x, 6, (3*log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x))))/(64*d) - (3*log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(64*d) + (5*cos(c + d*x))/(16*d*(3 - 5*sin(c + d*x)))]
    @test_int [1/(-3 + 5*sin(c + d*x))^3, x, 7, (43*log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x))))/(2048*d) - (43*log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(2048*d) - (5*cos(c + d*x))/(32*d*(3 - 5*sin(c + d*x))^2) + (45*cos(c + d*x))/(512*d*(3 - 5*sin(c + d*x)))]
    @test_int [1/(-3 + 5*sin(c + d*x))^4, x, 8, (279*log(cos((1/2)*(c + d*x)) - 3*sin((1/2)*(c + d*x))))/(32768*d) - (279*log(3*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(32768*d) + (5*cos(c + d*x))/(48*d*(3 - 5*sin(c + d*x))^3) - (25*cos(c + d*x))/(512*d*(3 - 5*sin(c + d*x))^2) + (995*cos(c + d*x))/(24576*d*(3 - 5*sin(c + d*x)))]

    @test_int [1/(-3 - 5*sin(c + d*x)), x, 4, log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x)))/(4*d) - log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(-3 - 5*sin(c + d*x))^2, x, 6, (3*log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(64*d) - (3*log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x))))/(64*d) - (5*cos(c + d*x))/(16*d*(3 + 5*sin(c + d*x)))]
    @test_int [1/(-3 - 5*sin(c + d*x))^3, x, 7, (43*log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(2048*d) - (43*log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x))))/(2048*d) + (5*cos(c + d*x))/(32*d*(3 + 5*sin(c + d*x))^2) - (45*cos(c + d*x))/(512*d*(3 + 5*sin(c + d*x)))]
    @test_int [1/(-3 - 5*sin(c + d*x))^4, x, 8, (279*log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(32768*d) - (279*log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x))))/(32768*d) - (5*cos(c + d*x))/(48*d*(3 + 5*sin(c + d*x))^3) + (25*cos(c + d*x))/(512*d*(3 + 5*sin(c + d*x))^2) - (995*cos(c + d*x))/(24576*d*(3 + 5*sin(c + d*x)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(c+d*x))^(n/2)=#

    @test_int [(a + b*sin(c + d*x))^(7/2), x, 8, -((2*b*(71*a^2 + 25*b^2)*cos(c + d*x)*sqrt(a + b*sin(c + d*x)))/(105*d)) - (24*a*b*cos(c + d*x)*(a + b*sin(c + d*x))^(3/2))/(35*d) - (2*b*cos(c + d*x)*(a + b*sin(c + d*x))^(5/2))/(7*d) + (32*a*(11*a^2 + 13*b^2)*Elliptic.E((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt(a + b*sin(c + d*x)))/(105*d*sqrt((a + b*sin(c + d*x))/(a + b))) - (2*(71*a^4 - 46*a^2*b^2 - 25*b^4)*Elliptic.F((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt((a + b*sin(c + d*x))/(a + b)))/(105*d*sqrt(a + b*sin(c + d*x)))]
    @test_int [(a + b*sin(c + d*x))^(5/2), x, 7, -((16*a*b*cos(c + d*x)*sqrt(a + b*sin(c + d*x)))/(15*d)) - (2*b*cos(c + d*x)*(a + b*sin(c + d*x))^(3/2))/(5*d) + (2*(23*a^2 + 9*b^2)*Elliptic.E((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt(a + b*sin(c + d*x)))/(15*d*sqrt((a + b*sin(c + d*x))/(a + b))) - (16*a*(a^2 - b^2)*Elliptic.F((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt((a + b*sin(c + d*x))/(a + b)))/(15*d*sqrt(a + b*sin(c + d*x)))]
    @test_int [(a + b*sin(c + d*x))^(3/2), x, 6, -((2*b*cos(c + d*x)*sqrt(a + b*sin(c + d*x)))/(3*d)) + (8*a*Elliptic.E((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt(a + b*sin(c + d*x)))/(3*d*sqrt((a + b*sin(c + d*x))/(a + b))) - (2*(a^2 - b^2)*Elliptic.F((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt((a + b*sin(c + d*x))/(a + b)))/(3*d*sqrt(a + b*sin(c + d*x)))]
    @test_int [(a + b*sin(c + d*x))^(1/2), x, 2, (2*Elliptic.E((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt(a + b*sin(c + d*x)))/(d*sqrt((a + b*sin(c + d*x))/(a + b)))]
    @test_int [1/(a + b*sin(c + d*x))^(1/2), x, 2, (2*Elliptic.F((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt((a + b*sin(c + d*x))/(a + b)))/(d*sqrt(a + b*sin(c + d*x)))]
    @test_int [1/(a + b*sin(c + d*x))^(3/2), x, 4, (2*b*cos(c + d*x))/((a^2 - b^2)*d*sqrt(a + b*sin(c + d*x))) + (2*Elliptic.E((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt(a + b*sin(c + d*x)))/((a^2 - b^2)*d*sqrt((a + b*sin(c + d*x))/(a + b)))]
    @test_int [1/(a + b*sin(c + d*x))^(5/2), x, 7, (2*b*cos(c + d*x))/(3*(a^2 - b^2)*d*(a + b*sin(c + d*x))^(3/2)) + (8*a*b*cos(c + d*x))/(3*(a^2 - b^2)^2*d*sqrt(a + b*sin(c + d*x))) + (8*a*Elliptic.E((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt(a + b*sin(c + d*x)))/(3*(a^2 - b^2)^2*d*sqrt((a + b*sin(c + d*x))/(a + b))) - (2*Elliptic.F((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt((a + b*sin(c + d*x))/(a + b)))/(3*(a^2 - b^2)*d*sqrt(a + b*sin(c + d*x)))]
    @test_int [1/(a + b*sin(c + d*x))^(7/2), x, 8, (2*b*cos(c + d*x))/(5*(a^2 - b^2)*d*(a + b*sin(c + d*x))^(5/2)) + (16*a*b*cos(c + d*x))/(15*(a^2 - b^2)^2*d*(a + b*sin(c + d*x))^(3/2)) + (2*b*(23*a^2 + 9*b^2)*cos(c + d*x))/(15*(a^2 - b^2)^3*d*sqrt(a + b*sin(c + d*x))) + (2*(23*a^2 + 9*b^2)*Elliptic.E((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt(a + b*sin(c + d*x)))/(15*(a^2 - b^2)^3*d*sqrt((a + b*sin(c + d*x))/(a + b))) - (16*a*Elliptic.F((1/2)*(c - pi/2 + d*x), (2*b)/(a + b))*sqrt((a + b*sin(c + d*x))/(a + b)))/(15*(a^2 - b^2)^2*d*sqrt(a + b*sin(c + d*x)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(c+d*x))^(n/3)=#

    @test_int [(a + b*sin(c + d*x))^(4/3), x, 3, -((sqrt(2)*(a + b)*AppellF1(1/2, 1/2, -(4/3), 3/2, (1/2)*(1 - sin(c + d*x)), (b*(1 - sin(c + d*x)))/(a + b))*cos(c + d*x)*(a + b*sin(c + d*x))^(1/3))/(d*sqrt(1 + sin(c + d*x))*((a + b*sin(c + d*x))/(a + b))^(1/3)))]
    @test_int [(a + b*sin(c + d*x))^(2/3), x, 3, -((sqrt(2)*AppellF1(1/2, 1/2, -(2/3), 3/2, (1/2)*(1 - sin(c + d*x)), (b*(1 - sin(c + d*x)))/(a + b))*cos(c + d*x)*(a + b*sin(c + d*x))^(2/3))/(d*sqrt(1 + sin(c + d*x))*((a + b*sin(c + d*x))/(a + b))^(2/3)))]
    @test_int [(a + b*sin(c + d*x))^(1/3), x, 3, -((sqrt(2)*AppellF1(1/2, 1/2, -(1/3), 3/2, (1/2)*(1 - sin(c + d*x)), (b*(1 - sin(c + d*x)))/(a + b))*cos(c + d*x)*(a + b*sin(c + d*x))^(1/3))/(d*sqrt(1 + sin(c + d*x))*((a + b*sin(c + d*x))/(a + b))^(1/3)))]
    @test_int [1/(a + b*sin(c + d*x))^(1/3), x, 3, -((sqrt(2)*AppellF1(1/2, 1/2, 1/3, 3/2, (1/2)*(1 - sin(c + d*x)), (b*(1 - sin(c + d*x)))/(a + b))*cos(c + d*x)*((a + b*sin(c + d*x))/(a + b))^(1/3))/(d*sqrt(1 + sin(c + d*x))*(a + b*sin(c + d*x))^(1/3)))]
    @test_int [1/(a + b*sin(c + d*x))^(2/3), x, 3, -((sqrt(2)*AppellF1(1/2, 1/2, 2/3, 3/2, (1/2)*(1 - sin(c + d*x)), (b*(1 - sin(c + d*x)))/(a + b))*cos(c + d*x)*((a + b*sin(c + d*x))/(a + b))^(2/3))/(d*sqrt(1 + sin(c + d*x))*(a + b*sin(c + d*x))^(2/3)))]
    @test_int [1/(a + b*sin(c + d*x))^(4/3), x, 3, -((sqrt(2)*AppellF1(1/2, 1/2, 4/3, 3/2, (1/2)*(1 - sin(c + d*x)), (b*(1 - sin(c + d*x)))/(a + b))*cos(c + d*x)*((a + b*sin(c + d*x))/(a + b))^(1/3))/((a + b)*d*sqrt(1 + sin(c + d*x))*(a + b*sin(c + d*x))^(1/3)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*sin(c+d*x))^n*with*n*symbolic=#

    @test_int [(a + b*sin(c + d*x))^n, x, 3, -((sqrt(2)*AppellF1(1/2, 1/2, -n, 3/2, (1/2)*(1 - sin(c + d*x)), (b*(1 - sin(c + d*x)))/(a + b))*cos(c + d*x)*(a + b*sin(c + d*x))^n)/(((a + b*sin(c + d*x))/(a + b))^n*(d*sqrt(1 + sin(c + d*x)))))]

    @test_int [(3 + 4*sin(c + d*x))^n, x, 2, -((sqrt(2)*7^n*AppellF1(1/2, 1/2, -n, 3/2, (1/2)*(1 - sin(c + d*x)), (4/7)*(1 - sin(c + d*x)))*cos(c + d*x))/(d*sqrt(1 + sin(c + d*x))))]
    @test_int [(3 - 4*sin(c + d*x))^n, x, 2, (sqrt(2)*7^n*AppellF1(1/2, -n, 1/2, 3/2, (4/7)*(1 + sin(c + d*x)), (1/2)*(1 + sin(c + d*x)))*cos(c + d*x))/(d*sqrt(1 - sin(c + d*x)))]

    @test_int [(4 + 3*sin(c + d*x))^n, x, 2, (sqrt(2)*AppellF1(1/2, 1/2, -n, 3/2, (1/2)*(1 + sin(c + d*x)), -3*(1 + sin(c + d*x)))*cos(c + d*x))/(d*sqrt(1 - sin(c + d*x)))]
    @test_int [(4 - 3*sin(c + d*x))^n, x, 2, (sqrt(2)*7^n*AppellF1(1/2, -n, 1/2, 3/2, (3/7)*(1 + sin(c + d*x)), (1/2)*(1 + sin(c + d*x)))*cos(c + d*x))/(d*sqrt(1 - sin(c + d*x)))]

    @test_int [(-3 + 4*sin(c + d*x))^n, x, 2, -((sqrt(2)*AppellF1(1/2, 1/2, -n, 3/2, (1/2)*(1 - sin(c + d*x)), 4*(1 - sin(c + d*x)))*cos(c + d*x))/(d*sqrt(1 + sin(c + d*x))))]
    @test_int [(-3 - 4*sin(c + d*x))^n, x, 2, (sqrt(2)*AppellF1(1/2, -n, 1/2, 3/2, 4*(1 + sin(c + d*x)), (1/2)*(1 + sin(c + d*x)))*cos(c + d*x))/(d*sqrt(1 - sin(c + d*x)))]

    @test_int [(-4 + 3*sin(c + d*x))^n, x, 3, (sqrt(2)*7^n*AppellF1(1/2, -n, 1/2, 3/2, (3/7)*(1 + sin(c + d*x)), (1/2)*(1 + sin(c + d*x)))*cos(c + d*x)*(-4 + 3*sin(c + d*x))^n)/((4 - 3*sin(c + d*x))^n*(d*sqrt(1 - sin(c + d*x))))]
    @test_int [(-4 - 3*sin(c + d*x))^n, x, 3, -((AppellF1(1 + n, 1/2, 1/2, 2 + n, 4 + 3*sin(c + d*x), (1/7)*(4 + 3*sin(c + d*x)))*cos(c + d*x)*(-4 - 3*sin(c + d*x))^(1 + n)*sqrt(-1 - sin(c + d*x)))/(sqrt(7)*d*(1 + n)*sqrt(1 - sin(c + d*x))*(1 + sin(c + d*x))))]
    end

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.1.1 (a+b sin)^n.jl


text:
## info command format:
## info time 203 nodes 11111111 score cp 11168 pv 5e9i+ ....

const PV = 1
const NonPV = 2
const Root = 3

function rememberPV(gs::GameStatus)
    # remember the last PV, and also the 5 previous ones because
    # they usually contain good moves to try
    i::Int = 0

    gs.lastPVLength = gs.triangularLength[1]
    #println("last length = ", gs.triangularLength[1])
    for i = 1:gs.triangularLength[1]
	gs.lastPV[i] = gs.triangularArray[1,i];
    end
end

# alpha-beta search
function AlphaBeta( gs::GameStatus, ply::Int, depth::Float64, alpha::Int, beta::Int)
    gs.triangularLength[ply+1] = ply
    ev_sign = ((ply & 0x01) == 0)?-1:1
    teban::Int = ((ply & 0x1) == 0)?gs.side:((gs.side==SENTE)?GOTE:SENTE)

    if depth <= 0.0
        return Eval( gs.board.nextMove, gs.board)
    end

    if gs.inodes > 100000
        gs.timedout = true
        return 0
    end

    gs.moveBufLen[ply+1+1] = generateMoves(gs.board, gs.moveBuf, teban, gs.moveBufLen[ply+1], gs)

    for i = gs.moveBufLen[ply+1]+1:gs.moveBufLen[ply+1+1] # 1 origin
        makeMove( gs.board, i, gs.moveBuf, teban)
        if in_check( teban, gs.board)
            #println("check!")
            takeBack( gs.board, i, gs.moveBuf, teban)
        else
            gs.inodes += 1
            val::Int = -AlphaBeta( gs, ply+1, depth-1.0, -beta, -alpha)
            takeBack( gs.board, i, gs.moveBuf, teban)
            #if val >= beta
            #    return beta
            #end
            if val > alpha
                alpha = val
		# both sides want to maximize from *their* perspective

                gs.triangularArray[ply+1,ply+1] = gs.moveBuf[i]
		# save this move
                for j = (ply +1+1):gs.triangularLength[ply+1+1]
                    gs.triangularArray[ply+1,j] = gs.triangularArray[ply+1+1,j]
                    # and append the latest best PV from deeper plies
                end
                gs.triangularLength[ply+1] = gs.triangularLength[ply + 1+1]
                rememberPV(gs)
            end
        end
    end
    alpha
end

# search driver
function think( sengo::Int, gs::GameStatus)
    move::Move = Move(0,0,0,0,0,0) # null move
    gs.lastPV = [Move(0,0,0,0,0,0) for x = 1:MaxPly]::Array{Move,1}
    gs.lastPVLength = 0
    gs.whiteHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.blackHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.nsStart = time_ns()
    gs.inodes = 0
    gs.timedout = false
    gs.side = sengo
    gs.board.nextMove = sengo

    for IDdepth = 1.0:1.0:15.0 #MaxPly
        gs.moveBufLen = [0 for x = 1:MaxPly]::Array{Int,1}
        gs.moveBuf = [Move(0,0,0,0,0,0) for x = 1:MaxMoves]
        gs.triangularLength = [0 for x = 1:MaxPly]::Array{Int,1}
        gs.triangularArray  = [Move(0,0,0,0,0,0) for x = 1:MaxPly, y = 1:MaxPly]::Array{Move,2}
        gs.followpv = true
        gs.allownull = true
        gs.depth = IDdepth
        # score::Int = AlphaBeta(gs, 0, IDdepth,-Infinity,Infinity) # ply=0
        score::Int = PVS(gs, 0, IDdepth,-Infinity,Infinity) # ply=0
        # println("Depth = ", IDdepth, ", Score = ", score)
        if gs.timedout
            return gs.lastPV[1]
        end
        rememberPV(gs)
        timeInMSecs::Int = int((time_ns() - gs.nsStart)/1000000)
        NPS::Int = int(gs.inodes / ((time_ns() - gs.nsStart)/1000000000))
        ## info time 203 nodes 11111111 score cp 11168 pv 5e9i+ ....
        print("info time ",timeInMSecs," depth ", int(IDdepth), " nodes ", gs.inodes, " score cp ", score, " nps ",NPS," pv")
        for i = 1:gs.triangularLength[1]
	    print(" ",move2USIString(gs.lastPV[i]))
        end
        println()
    end
    move = gs.lastPV[1]
    return move
end

function showPVforUSI(IDdepth::Float64, score::Int64, gs::GameStatus)
    i = 0::Int
    timeInMSecs = int((time_ns() - gs.nsStart)/1000000)::Int
    NPS = int(gs.inodes / ((time_ns() - gs.nsStart)/1000000000))::Int
    ## info time 203 nodes 11111111 score cp 11168 pv 5e9i+ ....
    print("info time ",timeInMSecs," depth ", int(IDdepth), " nodes ", gs.inodes, " score cp ", score, " nps ",NPS," pv")
    for i = 1:gs.triangularLength[1]
	print(" ",move2USIString(gs.lastPV[i]))
    end
    println()
end

# search driver
function thinkASP( sengo::Int, gs::GameStatus)
    move::Move = Move(0,0,0,0,0,0) # null move
    bestMove::Move = Move(0,0,0,0,0,0) # null move # added recently
    gs.lastPV = [Move(0,0,0,0,0,0) for x = 1:MaxPly]::Array{Move,1}
    gs.lastPVLength = 0
    gs.whiteHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.blackHeuristics = [0 for x = 1:0xff, y = 1:0xff]::Array{Int,2}
    gs.nsStart = time_ns()
    gs.inodes = 0
    gs.timedout = false
    gs.side = sengo
    gs.board.nextMove = sengo
    gs.tt = Dict{Uint64,TransP}()
    delta::Int = 48
    score::Int = 0
    middle::Int = 0
    smallAlpha::Int = 0
    smallBeta::Int  = 0

    ev_sign = gs.side == SENTE? 1: -1

    for IDdepth = 1.0:1.0:60.0 #MaxPly
        hitScore::Bool = false

        while hitScore == false
            gs.moveBufLen = [0 for x = 1:MaxPly]::Array{Int,1}
            gs.moveBuf = [Move(0,0,0,0,0,0) for x = 1:MaxMoves]
            gs.triangularLength = [0 for x = 1:MaxPly]::Array{Int,1}
            gs.triangularArray  = [Move(0,0,0,0,0,0) for x = 1:MaxPly, y = 1:MaxPly]::Array{Move,2}
            gs.followpv = true
            gs.allownull = true
            gs.MoveBeginIndex = 0

            # score::Int = AlphaBeta(gs, 0, IDdepth,-Infinity,Infinity) # ply=0
            if IDdepth == 1.0
                score = ev_sign * SimplePVS(gs, 0, IDdepth,-Infinity,Infinity,Root,false) # ply=0
                middle = score
                hitScore = true
                smallAlpha = middle - delta
                smallBeta  = middle + delta
                bestMove = gs.lastPV[1]
            else
                score = ev_sign * SimplePVS(gs, 0, IDdepth,smallAlpha,smallBeta,Root,false) # ply=0
                middle = score

                if smallAlpha == -Infinity
                    hitScore = true
                end
                if smallBeta  == Infinity
                    hitScore = true
                end

                if smallAlpha < score < smallBeta
                    hitScore = true
                elseif score <= smallAlpha
                    delta *= 2
                    smallAlpha = middle - delta
                    if delta >= 2000
                        smallAlpha = -Infinity
                    end
                    showPVforUSI(IDdepth,score,gs)
                elseif score >= smallBeta
                    delta *= 2
                    smallBeta = middle + delta
                    if delta >= 2000
                        smallBeta = Infinity
                    end
                    showPVforUSI(IDdepth,score,gs)
                end
            end
        end
        # println("Depth = ", IDdepth, ", Score = ", score)
        if gs.timedout
            return bestMove
        end
        rememberPV(gs)
        showPVforUSI(IDdepth,score,gs)
        bestMove = gs.lastPV[1]
    end
    move = bestMove
    return move
end

function ThinkTest(gs::GameStatus)
    m::Move = think(SENTE,gs)
end

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/JapaneseChess

max_forks_repo_name: kimrin/JapaneseChess


text:
    using Revise, Parameters, Logging, Test
    using CUDA
    includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
    using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
    using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils, ..Housdorff
    using ..ResultProcess

fpdummy = 100
fndummy = 100

resList,resListIndicies = allocateResultLists(fpdummy,fndummy)
reference= zeros(UInt8,100)
for i in 1:100
  resList[i,6]=i
  reference[i]=i
end
#only non 0 should be taken into acount

for i in 1:90
  resListIndicies[i]=i
end

function testGetAverage()
  getAverage(resList,resListIndicies,entriesPerBlock,totalLength,iterLoopResList ,globalSum )
return
end

@cuda threads=(16,16) blocks=1 testKernelForPaddingAnalysis(testArrIn,referenceArray,blockBeginingX,blockBeginingY,blockBeginingZ,  resArray,


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

using Markdown
using InteractiveUtils

# ╔═╡ ad4f98b4-4d34-40e6-a546-1013badd310a
begin
    import Pkg
    # activate the shared project environment
    Pkg.activate(Base.current_project())
    using Omega, Distributions, UnicodePlots
end

# ╔═╡ 798f9234-39f5-42dc-9fe7-48f7c77cf2d2
md"""
# Prelude: Thinking About Assembly Lines

Imagine a factory where the widget-maker is used to make widgets, but they are sometimes faulty. The tester tests them and lets through only the good ones. You don’t know what tolerance the widget tester is set to, and wish to infer it. We can represent this as:
"""

# ╔═╡ f9c9e3cf-cb7b-4591-94ce-a4aae8679c92
widget_machine_choice = Categorical([.05, .1, .2, .3, .2, .1, .05])

# ╔═╡ 9c0214f4-2c2d-4828-a77b-9b0d144f7bbc
widget_machine(i, ω) = [.2 , .3, .4, .5, .6, .7, .8 ][widget_machine_choice(i, ω)]

# ╔═╡ 4b961d20-f630-4b34-8809-884a4bdf26bc
actual_weights = [.6, .7, .8]

# ╔═╡ 01c5334f-ff4b-4087-96ac-57131e65886b
tolerance = @~ Uniform(0.3, 0.7)

# ╔═╡ dae8d29a-fa2f-4322-8cf1-23acabea44e3
function get_good_widget(i, ω)
	widget = (@uid, i) ~ widget_machine
	widget(ω) > tolerance(ω) ? widget(ω) : get_good_widget(i + 1, ω)
end

# ╔═╡ 68979be3-f3c4-4926-90d4-b095d42ccbb1
actual_widgets = [0.6, 0.7, 0.8]

# ╔═╡ 2c88252c-ba83-423e-a340-c963f8cae334
random_widgets = manynth(get_good_widget, 1:length(actual_widgets))

# ╔═╡ 772edb3b-9d81-4cf1-bc8c-7f8373334a6c
randsample(random_widgets)

# ╔═╡ 7734b963-5cd6-42e6-9c09-fa1f4691109f
tolerance_ = randsample(tolerance |ᶜ (random_widgets ==ₚ actual_widgets), 1000)

# ╔═╡ 7d411137-313e-473a-b9ef-4c82ccb3e2d8
histogram(tolerance_)

# ╔═╡ 38576ae3-8cee-4ba0-9eb4-d81c167c808b
md"But notice that the definition of getGoodWidget is exactly like the definition of rejection sampling! We can re-write this much more simply"

# ╔═╡ b90bf4ab-7a19-439d-bae8-5d6a834f0d9a
get_good_widget_simple = widget_machine |ᶜ widget > tolerance

# ╔═╡ 9ea0aab7-1228-43d5-9a68-47b4d117036c
md"# Social Cognition
How can we capture our intuitive theory of other people? Central to our understanding is the principle of rationality: an agent tends to choose actions that she expects to lead to outcomes that satisfy her goals. (This is a slight restatement of the principle as discussed in Baker et al. (2009), building on earlier work by Dennett (1989), among others.) We can represent this in WebPPL as an inference over actions—an agent reasons about actions that lead to their goal being satisfied:
"

# ╔═╡ 1d1b7607-e0f5-4234-b05b-40daba65bc36
action_prior = @~ Categorical([0.5, 0.5])

# ╔═╡ 96fbba9c-516e-47cc-a3bc-7eefcbfdfc17
function vending_machine(action)
	if action == 1
		"bagel"
	elseif action == 2
		"cookie"
	end
end

# ╔═╡ db8999c9-6293-4e3c-9a79-6b5c29770261
function choose_action(goal_state, transition)
	action_prior |ᶜ ω -> transition(action_prior(ω)) == goal_state
end

# ╔═╡ 94a623bd-0c48-478c-802a-21288ce6c9a5
randsample(choose_action("cookie", vending_machine))

# ╔═╡ 58ec1fc3-2ac8-4ed3-9028-e09dbc10cfc7
md"""We see, unsurprisingly, that if Sally wants a cookie, she will always press button b. In a world that is not quite so deterministic Sally’s actions will be more stochastic:"""

# ╔═╡ 1fab7ab9-258d-4526-ac80-52b107058fe2
pget(xs) = i -> xs[i]

# ╔═╡ deaa4d92-1863-4963-8a54-25b925240c5e
function vending_machine_stochastic(ω, action)
	choices = ["bagel", "cookie"]
	if action == 1
		choices[(@~ Categorical([0.9, 0.1]))(ω)]
	elseif action == 2
		choices[(@~ Categorical([0.1, 0.9]))(ω)]
	end
end

# ╔═╡ 64c941d8-900e-41a8-b5ab-3539b45dccfb
function choose_action_stochastic(goal_state, transition)
	action_prior |ᶜ ω -> transition(ω, action_prior(ω)) == goal_state
end

# ╔═╡ 041bd9be-7161-4e2c-8ea1-f1e93c63836e
action_samples = randsample(choose_action_stochastic("cookie", vending_machine_stochastic), 100)

# ╔═╡ 95cc8a2b-3280-4ed0-9d64-488062b09a10
histogram(action_samples)

# ╔═╡ 6f619867-9af2-47d9-9621-759b65593a84
md"## Inferring Goals
Now imagine that we don’t know Sally’s goal (which food she wants), but we observe her pressing button b. We can use Infer to infer her goal (this is sometimes called “inverse planning”, since the outer infer “inverts” the inference inside chooseAction)."

# ╔═╡ 344b85fc-d799-43a0-856c-297679084156
goal = pget(["bagel", "cookie"]) ∘ @~ Categorical([.5, .5])

# ╔═╡ 1f10beb3-9483-4619-83ba-49f743135898
action_dist(ω) = choose_action_stochastic(goal(ω), vending_machine_stochastic)(ω)

# ╔═╡ fcf06a26-35e4-47ac-878b-dae9b34182c8
goal_posterior = goal |ᶜ (action_dist ==ₚ 2)

# ╔═╡ 7d2f1234-84c3-4a63-87b9-f36b492fd3e0
randsample(goal_posterior, 100)

# ╔═╡ 402d913d-8788-4ade-a11e-d0f09853e5f5
md"Now let’s imagine a more ambiguous case: button b is “broken” and will (uniformly) randomly result in a food from the machine. If we see Sally press button b, what goal is she most likely to have?"

# ╔═╡ 0a4ac84e-f51a-41f3-88a2-65cb695b45cf

# ╔═╡ 3130e952-63ed-4236-bca6-f6629c38440c
md"## Inferring preferences

If we have some prior knowledge about Sally’s preferences (which goals she is likely to have) we can incorporate this immediately into the prior over goals (which above was uniform).

A more interesting situation is when we believe that Sally has some preferences, but we don’t know what they are. We capture this by adding a higher level prior (a uniform) over preferences. Using this we can learn about Sally’s preferences from her actions: after seeing Sally press button b several times, what will we expect her to want the next time?
"

# ╔═╡ 3dc36e8c-a494-4022-8c72-fd1b737d0301

# ╔═╡ Cell order:
# ╠═ad4f98b4-4d34-40e6-a546-1013badd310a
# ╟─798f9234-39f5-42dc-9fe7-48f7c77cf2d2
# ╠═f9c9e3cf-cb7b-4591-94ce-a4aae8679c92
# ╠═9c0214f4-2c2d-4828-a77b-9b0d144f7bbc
# ╠═4b961d20-f630-4b34-8809-884a4bdf26bc
# ╠═dae8d29a-fa2f-4322-8cf1-23acabea44e3
# ╠═01c5334f-ff4b-4087-96ac-57131e65886b
# ╠═68979be3-f3c4-4926-90d4-b095d42ccbb1
# ╠═2c88252c-ba83-423e-a340-c963f8cae334
# ╠═772edb3b-9d81-4cf1-bc8c-7f8373334a6c
# ╠═7734b963-5cd6-42e6-9c09-fa1f4691109f
# ╠═7d411137-313e-473a-b9ef-4c82ccb3e2d8
# ╟─38576ae3-8cee-4ba0-9eb4-d81c167c808b
# ╠═b90bf4ab-7a19-439d-bae8-5d6a834f0d9a
# ╟─9ea0aab7-1228-43d5-9a68-47b4d117036c
# ╠═1d1b7607-e0f5-4234-b05b-40daba65bc36
# ╠═96fbba9c-516e-47cc-a3bc-7eefcbfdfc17
# ╠═db8999c9-6293-4e3c-9a79-6b5c29770261
# ╠═94a623bd-0c48-478c-802a-21288ce6c9a5
# ╠═58ec1fc3-2ac8-4ed3-9028-e09dbc10cfc7
# ╠═1fab7ab9-258d-4526-ac80-52b107058fe2
# ╠═deaa4d92-1863-4963-8a54-25b925240c5e
# ╠═64c941d8-900e-41a8-b5ab-3539b45dccfb
# ╠═041bd9be-7161-4e2c-8ea1-f1e93c63836e
# ╠═95cc8a2b-3280-4ed0-9d64-488062b09a10
# ╟─6f619867-9af2-47d9-9621-759b65593a84
# ╠═344b85fc-d799-43a0-856c-297679084156
# ╠═1f10beb3-9483-4619-83ba-49f743135898
# ╠═fcf06a26-35e4-47ac-878b-dae9b34182c8
# ╠═7d2f1234-84c3-4a63-87b9-f36b492fd3e0
# ╠═402d913d-8788-4ade-a11e-d0f09853e5f5
# ╠═0a4ac84e-f51a-41f3-88a2-65cb695b45cf
# ╟─3130e952-63ed-4236-bca6-f6629c38440c
# ╠═3dc36e8c-a494-4022-8c72-fd1b737d0301

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Computes Stein discrepancy between a weighted sample and a target with gradlogdensity
#
# Args:
#   points - n x d array of sample points
#   weights - n x 1 array of real-valued weights associated with sample points
#     (default: equal weights)
#   gradlogdensity - a function taking point x in Rd and outputting grad log p(x)
#   operator - string in {"langevin", "riemannian-langevin"} indicating the Stein operator
#     to use (default: "langevin")
#   method - One of {"graph", "classical", "kernel"}, this uses different
#     methods in order to construct the different discrepancies.

function stein_discrepancy(; points=[],
                           gradlogdensity=nothing,
                           method="graph",
                           operator="langevin",
                           kwargs...)
    # Check arguments
    isempty(points) && error("Must provide non-empty array of sample points")
    isa(gradlogdensity, Function) ||
        error("Must specify gradlogdensity of type Function")

    if method in ["graph", "classical"]
        # call appropriate graph discrepancy
        if method == "graph"
            if operator == "langevin"
                langevin_graph_discrepancy(;points=points, gradlogdensity=gradlogdensity, kwargs...)
            elseif operator == "riemannian-langevin"
                riemannian_langevin_graph_discrepancy(;points=points, gradlogdensity=gradlogdensity, kwargs...)
            else
                error("unrecognized operator: $(operator)")
            end
        elseif method == "classical"
            if operator == "langevin"
                langevin_classical_discrepancy(;points=points, gradlogdensity=gradlogdensity, kwargs...)
            else
                error("unrecognized operator: $(operator)")
            end
        end
    elseif method == "kernel"
        if operator == "langevin"
            langevin_kernel_discrepancy(;points=points, gradlogdensity=gradlogdensity, kwargs...)
        else
            error("unrecognized operator: $(operator)")
        end
    else
        error("unrecognized method: $(method)")
    end
end

# stein_discrepancy aliases
function ksd(;kwargs...)
    stein_discrepancy(; method="kernel", kwargs...)
end

function gsd(;kwargs...)
    stein_discrepancy(; method="graph", kwargs...)
end


text:
#
## Needed for this script:
# - Energy calibration done
# - A/E calibration done

function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this = parse(Int64, this)
# num_workers, this = 10,1
println("Worker " * string(this) * " of " * string(num_workers))

#
# Load packages and functions
include("../src/init.jl")
include("../src/fct.jl")
include("../src/worker_fct.jl")

#
## Choose "2vbb" or "0vbb"
bb = "0vbb"

if bb == "2vbb"
    base_path = "../../waveforms/sim/2vbb_check/"
    output_path = "../../waveforms/sim/2vbb_AoE/"
else
    base_path = "../../waveforms/sim/0vbb_check/"
    output_path = "../../waveforms/sim/0vbb_AoE/"
end
include("../src/fct.jl")
include("../src/fitting-fct.jl")
set = "sim"

calib_filepath = "../dicts/calib.json"
calib = JSON.parsefile(calib_filepath)
cut_lib_filepath = "../dicts/cut_lib.json"
cut_lib = JSON.parsefile(cut_lib_filepath)
AE_cal_filepath = "../dicts/AE_cal.json"
AE_cal = JSON.parsefile(AE_cal_filepath)
sf_lib_filepath = "../dicts/sf_lib.json"
sf_lib = JSON.parsefile(sf_lib_filepath);

channels = eval_ch

channels = get_share_for_worker(channels, num_workers, this)

for ch in channels
    ch_str = lpad(ch, 2, "0")

    Base.run(`clear`)
    @info("This worker has: " * string(channels))
    @info("---")
    @info("Start with Ch$ch_str")

    files = glob(joinpath(base_path, ch_str * "-" * ctn[ch] * "/*.h5"));

    numofele = 5
    BackDelta5 = div(numofele,2)
    ForwardDelta5 = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1
    numofele = 201
    BackDelta201 = div(numofele,2)
    ForwardDelta201 = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1

    data = Table(A=[], E=[], E_unsmeared=[], energy=[], evtno=[], pos=[])

    @showprogress "Calculating A & E " for file in files
        tmp = HDF5.h5open(file, "r") do h5f
            LegendHDF5IO.readdata(h5f, "data")
        end

        A = []
        E = []
        for wf in tmp.waveform
            push!(A, maximum(movingaverage(diff(wf.value),5,BackDelta5,ForwardDelta5,3)))
            push!(E, maximum(movingaverage(diff(wf.value),201,BackDelta201,ForwardDelta201,13)))
        end
        E_cal = linmodel(E, calib[ctn[ch]][set]["lin_cal"][1])
        energy = []
        for e in sum.(tmp.energy)
            push!(energy, e.val)
        end
        append!(data, Table(A=A, E=E, E_unsmeared=E_cal, energy=energy, evtno=tmp.evtno, pos=tmp.pos))
    end

    smearing_dict = JSON.parsefile("../dicts/smearing.json")["Gauss"]
    p0 = [smearing_dict[string(ch)]["params"]["0"], smearing_dict[string(ch)]["params"]["1"], smearing_dict[string(ch)]["params"]["2"]]
    E = []
    @showprogress "Smearing energy " for e in data.E_unsmeared
        d = Normal(0, sqrt_fct(e,p0))
        push!(E, e + rand(d))
    end
    AoE = data.A ./ E
    AoE ./= calib[ctn[ch]][set]["AE_norm"]
    AoE ./= linmodel(E, AE_cal[ctn[ch]][set]["lin_fit"])
    AoE ./= AE_cal[ctn[ch]][set]["DEP_norm"]
    AoE .-= 1
    AoE ./= hypmodel(E, AE_cal[ctn[ch]][set]["sig_fit"])

    filename = joinpath(output_path, ch_str * "-" * ctn[ch] * "-AE_smeared.h5")
    !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
    HDF5.h5open(filename, "w") do h5f
        LegendHDF5IO.writedata(h5f, "data", Table(A = Array{Float64,1}(data.A),
                                                E = Array{Float64,1}(E),
                                                AoE = Array{Float64,1}(AoE),
                                                E_unsmeared = Array{Float64,1}(data.E_unsmeared),
                                                energy = Array{Float64,1}(data.energy),
                                                evtno = Array{Int64}(data.evtno),
                                                pos = VectorOfArrays(Array{typeof(data.pos[1])}(data.pos))
            )
        )
    end
end