text:
@testset "evaluate" begin

    @test evaluate(GPNode(Const(1)), []) == 1
    @test evaluate(GPNode(Const(1.25)), []) == 1.25

    @test evaluate(GPNode(Variable(1)), [3.14]) == 3.14
    @test evaluate(GPNode(Variable(1)), [3.14]) == 3.14

    add = PrimitiveFunction(+, 2)
    mul = PrimitiveFunction(*, 2)
    sub = PrimitiveFunction(-, 2)

    t = GPNode(mul, GPNode[
        GPNode(sub, GPNode[
            GPNode(Const(10)),
            GPNode(Variable(1))
        ]),
        GPNode(add, GPNode[
            GPNode(Variable(1))
            GPNode(Variable(2))
        ])
    ])
    t_func(vars) = (10 - vars[1]) * (vars[1] + vars[2])

    @test evaluate(t, [0, 1]) == t_func([0, 1])
    @test evaluate(t, [123, 456]) == t_func([123, 456])

    @test evaluate(t, [0., 1.]) == t_func([0., 1.])

end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
using Images

colors = [colorant"black", colorant"red", colorant"green", colorant"darkblue",
          colorant"purple", colorant"blue", colorant"yellow", colorant"white"]
wcol = 60 # width of each color bar
h, w = 150, wcol * length(colors) + 1
img = Matrix{RGB{N0f8}}(h, w);
for (j, col) in zip(1:wcol:w, colors)
    img[:, j:j+wcol] = col
end
save("data/colourbars.jpg", img)

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#!/usr/bin/env julia

using JLD2
using FileIO
using Cumulants

function main()
  d = try load("data/datafortests.jld2")
  catch
    println("please run gendata.jl")
    return ()
  end
  c = cumulants(d["data"], 6)
  save("data/cumulants.jld2", Dict("cumulants" => c))
end

main()


text:
export AbstractStateSpace, ndim, volume, isuniform, encode, decode, BooleanSpace

abstract type AbstractStateSpace end

@unimpl Base.size(space::AbstractStateSpace)

ndim(space::AbstractStateSpace) = space.ndim

volume(space::AbstractStateSpace) = space.volume
Base.length(space::AbstractStateSpace) = volume(space)

isuniform(space::AbstractStateSpace) = space.isuniform

@unimpl encode(space::AbstractStateSpace, state)

@unimpl decode(space::AbstractStateSpace, code)

@unimpl Base.iterate(space::AbstractStateSpace)
@unimpl Base.iterate(space::AbstractStateSpace, state)

struct BooleanSpace <: AbstractStateSpace
    ndim::Int
    volume::Int
    isuniform::Bool
    function BooleanSpace(ndim)
        if ndim < 1
            error("too few dimension ($ndim), should be at least 1")
        else
            new(ndim, 2^ndim, true)
        end
    end
end

Base.size(space::BooleanSpace) = tuple(fill(2, space.ndim)...)

volume(space::BooleanSpace) = space.volume

function encode(space::BooleanSpace, state::AbstractVector)
    code = 0
    for x in state
        code = (code << 1) | x
    end
    code
end

function decode(space::BooleanSpace, code::Int)
    state = Array{Int}(undef, space.ndim)
    for i in space.ndim:-1:1
        state[i] = code & 1
        code >>= 1
    end
    state
end

function Base.iterate(space::BooleanSpace)
    state = (1, zeros(Int, space.ndim))
    copy(state[2]), state
end

function Base.iterate(space::BooleanSpace, state::Tuple{Int, Vector{Int}})
    n, s = state
    if n != space.volume
        for i in space.ndim:-1:1
            if s[i] == 0
                s[i] = 1
                @views fill!(s[i+1:end], 0)
                break
            end
        end
        copy(s), (n + 1, s)
    end
end


text:
"""
    SeisBandPass(in; <keyword arguments>)

Apply a bandpass filter to seismic data. Input and output data is 2D in tx domain. Filter is applied in fx domain.

# Arguments
- `in`: Input 2D data array in tx domain. Time is first dimension.

# Keyword arguments
- `dt=0.004`: time sampling interval in secs.
- `fa=0`: lower frequency cut [Hz]
- `fb=0`: lower frequency banpass [Hz]
- `fc=60`: upper frequency bandpass [Hz]
- `fd=80`: upper frequency cut [Hz]

# Example
```
julia> d = SeisLinearEvents(); SeisPlotAmplitude(d,125,0.004);
julia> d_filter = SeisBandPass(d;dt=0.004,fa=2,fb=8,fc=12,fd=20); SeisPlotAmplitude(d_filter,125,0.004)
julia> SeisPlot([d d_filter],title="Data and Filtered data")

```
"""
function SeisBandPass(d;dt=0.004,fa=0,fb=0,fc=60,fd=80)

    nt =size(d,1)
	dn = reshape(d,nt,:)
    v = size(dn)
	ntr = v[2]
    nf = iseven(nt) ? nt : nt + 1
	df = 1/nf/dt
	nw = round(Int,nf/2) + 1

	if(fd*dt*nf < nw)
		iw_max = round(Int,floor(fd*dt*nf))
	else
		iw_max = round(Int,floor(0.5/dt))
	end

	dp = pad_first_axis(dn,nf)
	m = fft(dp,1)/sqrt(size(dp,1))
	if fa > 0.
		m[1,:] *= 0.
	end
	for iw=2:iw_max
		f = df*(iw-1)
		if (f<fa)
			m[iw,:] *= 0.
		elseif (f >= fa && f < fb)
			m[iw,:] *= (f-fa)/(fb-fa)
		elseif (f >= fb && f <= fc)
			m[iw,:] *= 1.
		elseif (f > fc && f <= fd)
			m[iw,:] *= 1. - (f-fc)/(fd-fc)
		else
			m[iw,:] *= 0.
		end
	end
	m[iw_max:end,:] .= 0.

	# symmetries
	for iw=nw+1:nf
		m[iw,:] = conj(m[nf-iw+2,:])
	end
	dn = real(bfft(m,1)/sqrt(size(m,1)))
	dout = dn[1:nt,1:ntr]
	return reshape(dout,v);
end

"""
    SeisBandPass(in,out,parameters; <keyword arguments>)

# Arguments
- `in::String`: Input file - Seis format
- `out::String`: Output file - Seis format.
- `parameters` : list of the keyword arguments for the function SeisBandPass.

# Keyword arguments
- `group="gather"` : Options are all, some or gather
- `key=["imx","imy"]` : Defines type of gather
- `itrace=1` : Initial trace number
- `ntrace=10000` : Total number of traces to process at once

# Example
```
julia> param = Dict(:fa>=0,:fb=>0,:fc=>60,:fd=>80)
julia> SeisBandPass(filein,fileout, param,group="all")
```

"""
function SeisBandPass(in::String,out::String,parameters;group="gather",key=["imx","imy"],itrace=1,ntrace=10000)

    fa = get(parameters,:fa,0)
    fb = get(parameters,:fb,0)
    fc = get(parameters,:fc,60)
    fd = get(parameters,:fd,80)

    ext = SeisMain.ReadTextHeader(in);
    dt = ext.d1
    parameters = Dict(:dt=>dt,:fa=>fa,:fb=>fb,:fc=>fc,:fd=>fd )

	SeisProcessFile(in,out,[SeisBandPass],[parameters];group=group,key=key,itrace=itrace,ntrace=ntrace)
end

function pad_first_axis(a,N1)
	n1 = size(a,1)
	nx = size(a,2)
	b = zeros(N1,nx)
	for ix = 1 : nx
		b[1:n1,ix] = a[:,ix]
	end
	return b
end


text:
module SoundPressureLevel

using Distributed
using Dates
using Random
using Soundcard
using Libaudio
using DeviceUnderTest

struct Instrument
    calibrator::String
    dbspl::Float64
    dba::Float64
    lastcal::Date
    mic::String
    preamp::String
    gainpreamp::Float64
    soundcard::String
end

function inst2str(x::Instrument)
    s = string(x)
    s = replace(s, "SoundPressureLevel.Instrument("=>"")
    s = replace(s, "\"" => "")
    s = replace(s, ")"=>"")
    s = replace(s, ':'=>'-')
    s = replace(s, ", "=>"_")
end

function addlatest(mm::Matrix, t, fs, root, id=Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"))
    r = Soundcard.record(round(Int, t * fs), mm, fs)
    p = replace(string(now()), [':','.']=>'-')
    Libaudio.wavwrite(joinpath(root, p * "+" * inst2str(id) * ".wav"), r, Int(fs), 32)
    return r
end

function getlatest(root, id=Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"))
    loc = ""
    tspan = Vector{DateTime}([now(), now()])
    archive = [(DateTime(String(split(basename(i),"+")[1]), DateFormat("y-m-dTH-M-S-s")), i) for i in Libaudio.list(root, ".wav")]
    sort!(archive, by=x->x[1], rev=true)
    for i in archive
        if String(split(basename(i[2]),"+")[2]) == inst2str(id) * ".wav"
            tspan[1] = i[1]
            loc = i[2]
            break
        end
    end
    loc, diff(tspan)[1]
end

"""
- 'fs': is the canonical sample rate in [8000, 16000, 44100, 48000, 96000, 192000]
        so to convert from precise to canonical use [findmin(abs.(av.-fs))[2]]
"""
function recording(f, y, ms::Matrix, mm::Matrix, fs, synchronous=true)
    if synchronous
        r = Soundcard.playrecord(y, ms, mm, fs)
    else
        out = "_splout.wav"
        Libaudio.wavwrite(out, DeviceUnderTest.mixer(y, ms), Int(fs), 32)
        try
            f[:init]()
            f[:readyplay](out)
            done = remotecall(f[:play], workers()[1])
            r = Soundcard.record(size(y,1), mm, fs)
            fetch(done)
        finally
            # rm(out, force=true)
        end
    end
    return convert(Matrix{Float64},r)
end

"""
    getdbaspread() set tartget level roughly to 'dbasetting' and measures the new level several times to find out the spreadout.
    if 'dbasetting' is NaN then the spreadout is measured at its 'gaininit' level
"""
function getdbaspread(
    f,
    symbol::Vector,
    rep::Int,
    gaininit,
    ms::Matrix,
    mm::Matrix,
    fs,
    dbasetting,
    nspread::Int,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0)

    @assert nprocs() > 1
    @assert size(ms, 1) == 1
    @assert size(mm, 2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pstnl, pstnd = getlatest(root, piston)
    pezol, pezod = getlatest(root, piezo)

    rootlog = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::vector): use latest calibration files:")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::vector): $pstnl")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::vector): $pezol")

    if maxdayadd > 0
        pstnd ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::vector): calibration is too old")
        pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::vector): calibration is too old")
    end
    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pstn, sr = Libaudio.wavread(pstnl, Float64)
    pezo, sr = Libaudio.wavread(pezol, Float64)

    m = length(symbol)
    n = round(Int, td*fs)
    x = zeros(m+n,1)
    x[1:m,1] = symbol * 10^(gaininit/20)
    y = recording(f, [zeros(round(Int,tcs*fs),1); repeat(x,rep,1)], ms, mm, fs, synchronous)

    nspread = max(1, nspread)
    val = zeros(Bool, 3+nspread)
    val[1], pstnspl = Libaudio.spl(pstn[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piston.dbspl+barocorrection)
    val[2], pezospl = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dbspl)
    val[3], pezodba = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::vector): calibration deviation $(abs(pstnspl[1]-pezospl[1])) dB")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::vector): levels = $(pstnspl[1]) dB(SPL), $(pezospl[1]) dB(SPL), $(pezodba[1]) dB(A)")

    gainadj = isnan(dbasetting) ? gaininit : (gaininit+(dbasetting-pezodba[1]))
    x[1:m,1] = symbol * 10^(gainadj/20)

    dbaspread = zeros(eltype(pezodba), nspread)
    for i = 1:nspread
        y = recording(f, [zeros(round(Int,tcs*fs),1); repeat(x,rep,1)], ms, mm, fs, synchronous)
        val[3+i], dbaspreads = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
        dbaspread[i] = dbaspreads[1]
    end
    return (all(val), gainadj, dbaspread)
end

function getdbaspread(
    f,
    source::Matrix,
    gaininit,
    ms::Matrix,
    mm::Matrix,
    fs,
    dbasetting,
    nspread::Int,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    syncatten = -12,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0,
    fm = fs)

    @assert nprocs() > 1
    @assert size(mm,2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pstnl, pstnd = getlatest(root, piston)
    pezol, pezod = getlatest(root, piezo)

    rootlog = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::matrix): use latest calibration files:")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::matrix): $pstnl")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::matrix): $pezol")

    if maxdayadd > 0
        pstnd ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::matrix): calibration is too old")
        pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::matrix): calibration is too old")
    end
    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pstn, sr = Libaudio.wavread(pstnl, Float64)
    pezo, sr = Libaudio.wavread(pezol, Float64)

    nspread = max(1, nspread)
    val = zeros(Bool, 4+2nspread)

    rate = synchronous ? fs : fm
    s = Libaudio.symbol_expsinesweep(800, 2000, 0.5, rate)
    x = Libaudio.encode_syncsymbol(tcs, s, td, 10^(gaininit/20) * source, rate, 1, syncatten)
    y = recording(f, x, ms, mm, fs, synchronous)
    val[1], symloc = Libaudio.decode_syncsymbol(y, s, td, size(source,1)/rate, fs)
    bl = symloc[1]
    br = symloc[1] + round(Int,size(source,1)/rate*fs) - 1

    val[2], pstnspl = Libaudio.spl(pstn[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piston.dbspl+barocorrection)
    val[3], pezospl = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dbspl)
    val[4], pezodba = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::matrix): calibration deviation $(abs(pstnspl[1]-pezospl[1])) dB")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdbaspread(::matrix): levels = $(pstnspl[1]) dB(SPL), $(pezospl[1]) dB(SPL), $(pezodba[1]) dB(A)")

    gainadj = isnan(dbasetting) ? gaininit : (gaininit+(dbasetting-pezodba[1]))
    x = Libaudio.encode_syncsymbol(tcs, s, td, 10^(gainadj/20) * source, rate, 1, syncatten)

    dbaspread = zeros(eltype(pezodba), nspread)
    for i = 1:nspread
        y = recording(f, x, ms, mm, fs, synchronous)
        val[4+2i-1], symloc = Libaudio.decode_syncsymbol(y, s, td, size(source,1)/rate, fs)
        bl = symloc[1]
        br = symloc[1] + round(Int,size(source,1)/rate*fs) - 1
        val[4+2i], dbaspreads = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
        dbaspread[i] = dbaspreads[1]
    end

    return (all(val), gainadj, dbaspread)
end

function getdba(
    f,
    symbol::Vector,
    rep::Int,
    gaininit,
    ms::Matrix,
    mm::Matrix,
    fs,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0)

    @assert nprocs() > 1
    @assert size(ms, 1) == 1
    @assert size(mm, 2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pstnl, pstnd = getlatest(root, piston)
    pezol, pezod = getlatest(root, piezo)

    rootlog = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): use latest calibration files:")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): $pstnl")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): $pezol")

    if maxdayadd > 0
        pstnd ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): calibration is too old")
        pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): calibration is too old")
    end
    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pstn, sr = Libaudio.wavread(pstnl, Float64)
    pezo, sr = Libaudio.wavread(pezol, Float64)

    m = length(symbol)
    n = round(Int, td*fs)
    x = zeros(m+n,1)

    x[1:m,1] = symbol * 10^(gaininit/20)
    y = recording(f, [zeros(round(Int,tcs*fs),1); repeat(x,rep,1)], ms, mm, fs, synchronous)

    val = Bool[0,0,0]
    val[1], pstnspl = Libaudio.spl(pstn[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piston.dbspl+barocorrection)
    val[2], pezospl = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dbspl)
    val[3], dbainit = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): calibration deviation $(abs(pstnspl[1]-pezospl[1])) dB")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): levels = $(pstnspl[1]) dB(SPL), $(pezospl[1]) dB(SPL), $(dbainit[1]) dB(A)")

    return (all(val), dbainit[1])
end

function setdba(
    f,
    symbol::Vector,
    rep::Int,
    gaininit,
    dbainit,
    dbasetting,
    dbatol,
    ms::Matrix,
    mm::Matrix,
    fs,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0)

    @assert nprocs() > 1
    @assert size(ms, 1) == 1
    @assert size(mm, 2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pezol, pezod = getlatest(root, piezo)
    maxdayadd > 0 && pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(joinpath(Libaudio.folder(), Libaudio.logfile()), :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): calibration is too old")

    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pezo, sr = Libaudio.wavread(pezol, Float64)

    m = length(symbol)
    n = round(Int, td*fs)
    x = zeros(m+n,1)

    gainadj = gaininit
    dbaadj = dbainit

    while true
        gainadj = gainadj+(dbasetting-dbaadj)
        x[1:m,1] = symbol * 10^(gainadj/20)
        y = recording(f, [zeros(round(Int,tcs*fs),1); repeat(x,rep,1)], ms, mm, fs, synchronous)
        val, dbanow = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
        dbaadj = dbanow[1]

        if val && (dbasetting + dbatol ≥ dbaadj ≥ dbasetting - dbatol)
            Libaudio.printl(joinpath(Libaudio.folder(), Libaudio.logfile()), :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): tolerance satisfied")
            break
        else
            Libaudio.printl(joinpath(Libaudio.folder(), Libaudio.logfile()), :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): deviation = $(abs(dbaadj-dbasetting)), retry...")
        end
    end
    return (gainadj, dbaadj)
end

function getdba(
    f,
    source::Matrix,
    gaininit,
    ms::Matrix,
    mm::Matrix,
    fs,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    syncatten = -12,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0,
    fm = fs)

    @assert nprocs() > 1
    @assert size(mm,2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pstnl, pstnd = getlatest(root, piston)
    pezol, pezod = getlatest(root, piezo)

    rootlog = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::matrix): use latest calibration files:")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::matrix): $pstnl")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::matrix): $pezol")

    if maxdayadd > 0
        pstnd ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdba(::matrix): calibration is too old")
        pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.getdba(::matrix): calibration is too old")
    end
    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pstn, sr = Libaudio.wavread(pstnl, Float64)
    pezo, sr = Libaudio.wavread(pezol, Float64)

    val = Bool[0,0,0,0]

    rate = synchronous ? fs : fm
    s = Libaudio.symbol_expsinesweep(800, 2000, 0.5, rate)
    x = Libaudio.encode_syncsymbol(tcs, s, td, 10^(gaininit/20) * source, rate, 1, syncatten)
    y = recording(f, x, ms, mm, fs, synchronous)
    val[1], symloc = Libaudio.decode_syncsymbol(y, s, td, size(source,1)/rate, fs)
    bl = symloc[1]
    br = symloc[1] + round(Int,size(source,1)/rate*fs) - 1

    val[2], pstnspl = Libaudio.spl(pstn[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piston.dbspl+barocorrection)
    val[3], pezospl = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dbspl)
    val[4], dbainit = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")

    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::matrix): calibration deviation $(abs(pstnspl[1]-pezospl[1])) dB")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.getdba(::vector): levels = $(pstnspl[1]) dB(SPL), $(pezospl[1]) dB(SPL), $(dbainit[1]) dB(A)")

    return (all(val), dbainit[1])
end

function setdba(
    f,
    source::Matrix,
    gaininit,
    dbainit,
    dbasetting,
    dbatol,
    ms::Matrix,
    mm::Matrix,
    fs,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    syncatten = -12,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0,
    fm = fs)

    @assert nprocs() > 1
    @assert size(mm,2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pezol, pezod = getlatest(root, piezo)
    maxdayadd > 0 && pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(joinpath(Libaudio.folder(), Libaudio.logfile()), :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): calibration is too old")

    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pezo, sr = Libaudio.wavread(pezol, Float64)
    rate = synchronous ? fs : fm

    gainadj = gaininit
    dbaadj = dbainit

    while true
        gainadj = gainadj+(dbasetting-dbaadj)
        s = Libaudio.symbol_expsinesweep(800, 2000, 0.5, rate)
        x = Libaudio.encode_syncsymbol(tcs, s, td, 10^(gainadj/20) * source, rate, 1, syncatten)
        y = recording(f, x, ms, mm, fs, synchronous)
        vald, symloc = Libaudio.decode_syncsymbol(y, s, td, size(source,1)/rate, fs)
        bl = symloc[1]
        br = symloc[1] + round(Int,size(source,1)/rate*fs) - 1
        vals, dbanow = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
        dbaadj = dbanow[1]

        if vald && vals && (dbasetting + dbatol ≥ dbaadj ≥ dbasetting - dbatol)
            Libaudio.printl(joinpath(Libaudio.folder(), Libaudio.logfile()), :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): tolerance satisfied")
            break
        else
            Libaudio.printl(joinpath(Libaudio.folder(), Libaudio.logfile()), :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): deviation = $(abs(dbaadj-dbasetting)), retry...")
        end
    end
    return (gainadj, dbaadj)
end

#DEPRECATE

"""

# Note
 1) 'symbol' is the segment of signal for level measurement
 2) 'rep' if for multiple trial --- t_context + (symbol + decay) x repeat
 3) 'root' is the path for reference mic recordings of the calibrators (piston and piezo etc...)
 4) validation method 1: compare against the spl meter
    validation method 2: 200hz -> 10dB lower than dBSPL, 1kHz-> the same, 6kHz-> almost the same, 7kHz-> 0.8 dB lower than dBSPL
 5) 'dbasetting' if given NaN then no gain adjustment is applied and the dba measurement corresponds to 'gaininit'
"""
function setdba(
    f,
    symbol::Vector,
    rep::Int,
    gaininit,
    ms::Matrix,
    mm::Matrix,
    fs,
    dbasetting,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0)

    @assert nprocs() > 1
    @assert size(ms, 1) == 1
    @assert size(mm, 2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pstnl, pstnd = getlatest(root, piston)
    pezol, pezod = getlatest(root, piezo)

    rootlog = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): use latest calibration files:")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): $pstnl")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): $pezol")

    if maxdayadd > 0
        pstnd ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): calibration is too old")
        pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): calibration is too old")
    end
    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pstn, sr = Libaudio.wavread(pstnl, Float64)
    pezo, sr = Libaudio.wavread(pezol, Float64)

    m = length(symbol)
    n = round(Int, td*fs)
    x = zeros(m+n,1)
    x[1:m,1] = symbol * 10^(gaininit/20)
    y = recording(f, [zeros(round(Int,tcs*fs),1); repeat(x,rep,1)], ms, mm, fs, synchronous)

    val = Bool[0,0,0,0]
    val[1], pstnspl = Libaudio.spl(pstn[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piston.dbspl+barocorrection)
    val[2], pezospl = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dbspl)
    if abs(pstnspl[1]-pezospl[1]) > 0.5
        Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): calibration deviation > 0.5 dB")
    else
        Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): calibration deviation $(abs(pstnspl[1]-pezospl[1])) dB")
    end

    val[3], pstndba = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
    gainadj = isnan(dbasetting) ? gaininit : (gaininit+(dbasetting-pstndba[1]))

    x[1:m,1] = symbol * 10^(gainadj/20)
    y = recording(f, [zeros(round(Int,tcs*fs),1); repeat(x,rep,1)], ms, mm, fs, synchronous)
    val[4], pstndba = Libaudio.spl(pezo[:,1], y, symbol, rep, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")

    return (all(val), gainadj, pstndba[1])
end

# note: source is multichannel sound tracks for spl measurement, it is based on async method, therefore no need for parameter repeat
function setdba(
    f,
    source::Matrix,
    gaininit,
    ms::Matrix,
    mm::Matrix,
    fs,
    dbasetting,
    root,
    piston = Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    piezo = Instrument("42AB",114,114,Date("2018-07-24"),"26XX","12AA",0,"UFX"),
    synchronous = true,
    syncatten = -12,
    tcs = 3.0,
    td = 2.0,
    maxdayadd = 1,
    maxdaycal = 180,
    barocorrection = 0.0,
    fm = fs)

    @assert nprocs() > 1
    @assert size(mm,2) == 1
    @assert now() - DateTime(piston.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))
    @assert now() - DateTime(piezo.lastcal) ≤ Dates.Millisecond(Dates.Day(maxdaycal))

    pstnl, pstnd = getlatest(root, piston)
    pezol, pezod = getlatest(root, piezo)

    rootlog = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): use latest calibration files:")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): $pstnl")
    Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): $pezol")

    if maxdayadd > 0
        pstnd ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): calibration is too old")
        pezod ≥ Dates.Millisecond(Dates.Day(maxdayadd)) && Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): calibration is too old")
    end
    wf = Libaudio.WindowFrame(fs,16384,16384÷4)
    pstn, sr = Libaudio.wavread(pstnl, Float64)
    pezo, sr = Libaudio.wavread(pezol, Float64)

    val = Bool[0,0,0,0,0,0]
    rate = synchronous ? fs : fm
    s = Libaudio.symbol_expsinesweep(800, 2000, 0.5, rate)
    x = Libaudio.encode_syncsymbol(tcs, s, td, 10^(gaininit/20) * source, rate, 1, syncatten)
    y = recording(f, x, ms, mm, fs, synchronous)
    val[1], symloc = Libaudio.decode_syncsymbol(y, s, td, size(source,1)/rate, fs)
    bl = symloc[1]
    br = symloc[1] + round(Int,size(source,1)/rate*fs) - 1

    val[2], pstnspl = Libaudio.spl(pstn[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piston.dbspl+barocorrection)
    val[3], pezospl = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dbspl)
    if abs(pstnspl[1]-pezospl[1]) > 0.5
        Libaudio.printl(rootlog, :light_red, Libaudio.nows() * " | SoundPressureLevel.setdba(::vector): calibration deviation > 0.5 dB")
    else
        Libaudio.printl(rootlog, :light_yellow, Libaudio.nows() * " | SoundPressureLevel.setdba(::matrix): calibration deviation $(abs(pstnspl[1]-pezospl[1])) dB")
    end

    val[4], pstndba = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
    gainadj = isnan(dbasetting) ? gaininit : (gaininit+(dbasetting-pstndba[1]))

    x = Libaudio.encode_syncsymbol(tcs, s, td, 10^(gainadj/20) * source, rate, 1, syncatten)
    y = recording(f, x, ms, mm, fs, synchronous)
    val[5], symloc = Libaudio.decode_syncsymbol(y, s, td, size(source,1)/rate, fs)
    bl = symloc[1]
    br = symloc[1] + round(Int,size(source,1)/rate*fs) - 1

    val[6], pstndba = Libaudio.spl(pezo[:,1], y[bl:br,:], y[bl:br,1], 1, wf, 0, 0, 100, 12000, piezo.dba, weighting="A")
    return (all(val), gainadj, pstndba[1])
end

end # module


text:
using Coverage, AzManagers
coverage = process_folder(pkgdir(AzManagers))
LCOV.writefile("lcov.info", coverage)


text:
using LinearAlgebra: sqrt
using Printf: Threads
using Plots: Threads
using ProgressMeter
using Printf
include("PolygonBase.jl")
include("ESBase.jl")
include("silhouette&pieces.jl")

# X[1つめのピースのx座標, 1つめのピースのy座標, １つめの回転角Θ, 2つめのx座標, ...]
function loss_poly(X::Array, loss_args::Array{Float64, 1})
    global pices, silhouette

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    tmp = LibGEOS.intersection(unioned, MYPolygon2LibGEOS(silhouette))
    A = LibGEOS.area(tmp) / LibGEOS.area(MYPolygon2LibGEOS(silhouette))

    # 各2pieceの共通部分を求める
    inter = 1.0 - LibGEOS.area(unioned) / sum_pieces
    outer = 1.0 - LibGEOS.area(tmp)/sum_pieces

    # 頂点の誤差成分 (シルエットベース 緩い)
    Δv = 0.0
    for v_s in 1:silhouette.n
        Δmin = 100.0
        for j in 1:length(pieces)
            tmp_p = move(pieces[j], X[3j - 2], X[3j - 1])
            rotate!(tmp_p, X[3j] * 2π)
            for v_p in 1:tmp_p.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - tmp_p.vertexes[:, v_p]).^2))
            end
        end
        Δv += Δmin / silhouette.n
        # Δv = max(Δv, Δmin)
    end

    # 頂点の誤差成分 (ピースベース 厳しい)
    """
    Δv = 0.0
    for i in 1:length(pieces)
        p = move(pieces[i], X[3i - 2], X[3i - 1])
        rotate!(p, X[3i] * 2π)
        Δmax = 0.0
        # 頂点ごとに誤差を見る
        for v_p in 1:p.n
            Δmin = 100.0
            # シルエットについて見る
            for v_s in 1:silhouette.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - p.vertexes[:, v_p]).^2))
            end
            # その他のピースについて見る
            for j in 1:length(pieces)
                if j == i continue end
                q = move(pieces[j], X[3j - 2], X[3j - 1])
                rotate!(p, X[3j] * 2π)
                for v_q in 1:q.n
                    Δmin = min(Δmin, sum((q.vertexes[:,v_q] - p.vertexes[:, v_p]).^2))
                end
            end
            Δmax = max(Δmax, Δmin)
        end
        # Δv += 2tanh(2Δmin) / (π*silhouette.n)
        # Δv = max(Δv, Δmin)
        Δv += Δmax
    end
    """

    # 辺の誤差成分
    E = 0.0
    for i in 1:length(pieces)
        P = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(P, X[3i] * 2π)

        cosθ = 0.0

        # シルエットの各頂点と距離が最小な点を探す
        min_d = 100.0
        for p in 0:P.n-1, s in 0:silhouette.n-1
            d = sum((P.vertexes[:, p + 1] - silhouette.vertexes[:, s + 1]).^2)
            # 最小の大きさの角度を探す
            if d < min_d
                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ1 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ2 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ3 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ4 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                cosθ = max(cosθ1, cosθ2, cosθ3, cosθ4)
            end
        end
        E += cosθ / length(pieces)
    end

    return -loss_args[1] * A + loss_args[2] * Δv + loss_args[3] * outer + loss_args[4] * (E - 1)
end

# X[1つめのピースのx座標, 1つめのピースのy座標, １つめの回転角Θ, 2つめのx座標, ...]
function loss_poly_light(X::Array)
    global pices, silhouette

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    tmp = LibGEOS.intersection(unioned, MYPolygon2LibGEOS(silhouette))
    A = LibGEOS.area(tmp) / LibGEOS.area(MYPolygon2LibGEOS(silhouette))

    return -A * 100
end

silhouette = hexagon_5
pieces = [tri_m, tri_s, tri_s, square_s, parallelogram]

# loss_args = [100.0, 0.0, 0.0, 0.0]
loss_args = [100.0, 280.0, 0.0, 40.0]

# 初期化
max_gen = 168
sample_num = 128

# 最終結果を保存するか否か
save_results = true
if save_results
    rm("outputs", force=true, recursive=true)
    mkdir("outputs")

    best_fitnesses = Array{Float64, 1}()
    correct_time = 0
    p = Progress(sample_num)
    for t in 1:sample_num
        global best_fitnesses, correct_time
        local cmaes, rng
        cmaes = init_CMAES(zeros(3 * length(pieces)), 0.8, 18)
        fitnesses = zeros(cmaes.dim)
        X = zeros((cmaes.dim, cmaes.λ))
        rng = MersenneTwister(t)

        for gen in 1:max_gen
            # 個体生成
            X = samplePopulation(cmaes, rng=rng)

            # 回転角のパラメータを -π ~ π までに正規化
            for j in 1:size(X)[2]
                for i in 3:3:size(X)[1]
                    X[i, j] = mod(X[i, j], 1.0)
                end
            end

            # 評価値
            # loss(x) = loss_poly_light(x)
            loss(x) = loss_poly(x, loss_args)
            fitnesses = get_fitness(X, loss)

            # 更新
            update!(cmaes, X, fitnesses, gen)
        end
        push!(best_fitnesses, maximum(fitnesses))

        if -best_fitnesses[end] > 80
            correct_time += 1
        end

        # 表示用
        best_arg = argmin(fitnesses)
        best_pieces = Array{MYPolygon, 1}()
        for i in 1:length(pieces)
            todisplay = move(pieces[i], X[3i-2, best_arg], X[3i-1, best_arg])
            rotate!(todisplay, X[3i, best_arg] * 2π)
            push!(best_pieces, todisplay)
        end
        display(silhouette, best_pieces...)
        savefig(@sprintf "outputs\\%02.0f_%d_results" -best_fitnesses[end] t)

        next!(p)
    end
else
    best_fitnesses = Array{Float64, 1}()
    correct_time = Threads.Atomic{Int}(0)
    p = Progress(sample_num÷Threads.nthreads())
    Threads.@threads for t in 1:sample_num
        global best_fitnesses, correct_time
        local cmaes
        cmaes = init_CMAES(zeros(3 * length(pieces)), 1.0, 0)
        fitnesses = zeros(cmaes.dim)
        X = zeros((cmaes.dim, cmaes.λ))
        rng = MersenneTwister(t)

        for gen in 1:max_gen
            # 個体生成
            X = samplePopulation(cmaes, rng=rng)

            # 回転角のパラメータを -π ~ π までに正規化
            for j in 1:size(X)[2]
                for i in 3:3:size(X)[1]
                    X[i, j] %= mod(X[i, j], 1.0)
                end
            end

            # 評価値
            # loss(x) = loss_poly(x, loss_args)
            loss(x) = loss_poly_light(x)
            fitnesses = get_fitness(X, loss)

            # 更新
            update!(cmaes, X, fitnesses, gen)
        end
        push!(best_fitnesses, maximum(fitnesses))

        if -best_fitnesses[end] > 94
            Threads.atomic_add!(correct_time, 1)
        end

        if Threads.threadid() == 1
            next!(p)
        end
    end
end

# pl = histogram(-best_fitnesses, bins = 40:5:100, xlabel="Fitness", ylabel="The number of answers", lab=false)

# println("correct prob: $(correct_time[])/$(sample_num) = $(100correct_time[]/sample_num)%")

# savefig(pl, "fitness_distribusion_$(sample_num)samples_$(correct_time[])-corrects.png")
# Base.display(pl)

max_stars_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_issues_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_forks_repo_name: toshiki-kosen/SolvingTangramPuzzle


text:
abstract type SequentialSamplingModel <: ContinuousUnivariateDistribution end

function Base.show(io::IO, ::MIME"text/plain", model::SequentialSamplingModel)
    values = [getfield(model, f) for f in fieldnames(typeof(model))]
    values = map(x->typeof(x)== Bool ? string(x) : x, values)
    T = typeof(model)
    model_name = string(T.name.name)
    return pretty_table(
        values;
        title=model_name,
        row_name_column_title="Parameter",
        compact_printing=false,
        header=["Value"],
        row_name_alignment=:l,
        row_names=[fieldnames(typeof(model))...],
        formatters=ft_printf("%5.2f"),
        alignment=:l,
    )
end


text:
include("test_common.jl")
using Random

@testset "arrays" begin
    new_context()

    arr1 = collect(1.0:20.0)
    arr2 = reinterpret(Float32, arr1)
    bytes = reinterpret(UInt8, arr2)
    # JavaScript's arrays start with 0.
    buffer = node"(arr1, arr2) => {
        arr1[0] = 5
        arr2[3] = 4
        return arr1.buffer
    }"(arr1, arr2)
    @test buffer == bytes
    @test pointer(buffer) == pointer(arr2)
    @test arr1[1:2] == [5, 512]
    set_copy_array(true)
    js_arr = node_value(arr1)
    @test length(js_arr) == 20
    # `in` is checking indices.
    @test 0 in js_arr
    @test js_arr[0] == 5
    arr3 = Array(js_arr)
    set_copy_array(false)
    @test pointer(arr1) != pointer(arr3)
    js_arraybuffer = get(js_arr, "buffer", result=RESULT_NODE)
    @test length(js_arraybuffer) == 160
    dataview = @new node"DataView"(js_arraybuffer)
    @test length(dataview) == 160

    js_buffer = node"Buffer.alloc(5)"o
    @test value(Vector{UInt8}, js_buffer) == zeros(UInt8, 5)
    @test length(js_buffer) == 5
    @test Tuple(js_buffer) == Tuple(zeros(UInt8, 5))
    @test value(Tuple{Int, Int}, node"[1,2]"o) == (1, 2)

    delete_context()
end


text:
# include the functions of the module
using LinearAlgebra

include("simula01.jl")
include("mac.jl")
include("modalparam.jl")
include("modalparami.jl")


text:
const dx=0.05
const dt = 0.02;

function Fmatrix2{T}(F::AbstractArray{T,2},params::Vector, bin_centers,absorbing=0)

    stim=params[1]
    c2=params[2]
    c4=params[3]
    c6=params[4]
    sigma2=params[5]
    bias_d=params[6]
    hb_d=params[7]
    # println("sigmaaaa ",sigma2)
    # sigma2_sbin = convert(typeof(sigma2), sigma2)
    sigma2_sbin = convert(Float64, sigma2)

    # println(typeof(sigma2))
    # println(size(sigma2))
    # println(sigma2)
    # println("Converted ", sigma2_sbin)
    # println(typeof(sigma2_sbin))

    n_sbins = max(70, ceil(10*sqrt(sigma2_sbin)/dx))
    #println(n_sbins,sigma2_sbin)
    if absorbing==1
        F[1,1] = 1;
        F[end,end] = 1;
    else
        F[2,1] = 1;
        F[end-1,end] = 1;
    end
    swidth = 5.*sqrt(sigma2_sbin)

    sbinsize = swidth/n_sbins #sbins[2] - sbins[1]
    base_sbins    = collect(-swidth:sbinsize:swidth)

    ps       = exp.(-base_sbins.^2/(2*sigma2)) # exp(Array) -> exp.(x)
    ps       = ps/sum(ps);

    sbin_length = length(base_sbins)
    binN = length(bin_centers)

    mu = 0.
    for j in 2:binN-1
        mu = bin_centers[j] - (-(stim+bias_d+hb_d)-c2*bin_centers[j]+c4*bin_centers[j]^3+c6*bin_centers[j]^5)*dt#(exp(lam*dt))
        # mu=bin_centers[j] +stim*dt

        for k in 1:sbin_length
            sbin = (k-1)*sbinsize + mu - swidth

            if sbin <= bin_centers[1] #(bin_centers[1] + bin_centers[2])/2
                F[1,j] = F[1,j] + ps[k]
            elseif bin_centers[end] <= sbin#(bin_centers[end]+bin_centers[end-1])/2 <= sbins[k]
                F[end,j] = F[end,j] + ps[k]
            else # more condition
                if (sbin > bin_centers[1] && sbin < bin_centers[2])
                    lp = 1; hp = 2;
                elseif (sbin > bin_centers[end-1] && sbin < bin_centers[end])
                    lp = binN-1; hp = binN;
                else
                    lp = floor(Int,((sbin-bin_centers[2])/dx)) + 2#find(bin_centers .<= sbins[k])[end]
                    hp = ceil(Int,((sbin-bin_centers[2])/dx)) + 2#lp+1#Int(ceil((sbins[k]-bin_centers[2])/dx) + 1);
                end

                if lp == hp
                    F[lp,j] = F[lp,j] + ps[k]
                else
                    F[hp,j] = F[hp,j] + ps[k]*(sbin - bin_centers[lp])/(bin_centers[hp] - bin_centers[lp])
                    F[lp,j] = F[lp,j] + ps[k]*(bin_centers[hp] - sbin)/(bin_centers[hp] - bin_centers[lp])
                end
            end
        end
    end
end

function make_bins{T}(bins, B, dx::T, binN)
    cnt = 1
    for i=-binN:binN
        #println(cnt)
        bins[cnt] = i*dx
        cnt = cnt+1
    end

    if binN*dx == B
        bins[end] = B + dx
        bins[1] = -B - dx
    else
        bins[end] = 2*B - (binN-1)*dx
        bins[1] = -2*B + (binN-1)*dx
    end
end

function pdf_x_time_1stim(s,Tframe,param;B=0,hb=0)

    Nframes=length(s)
    stim_trans=zeros(Nframes)
    sigma_a=param["sigma_a"][1]

    Nsteps=convert(Int,floor(Tframe/dt))

    c2=zeros(typeof(sigma_a),Nframes)
    c4=zeros(typeof(sigma_a),Nframes)
    c6=zeros(typeof(sigma_a),Nframes)
    bias_d=zeros(typeof(sigma_a),Nframes)
    hbias_d=zeros(typeof(sigma_a),Nframes)

    binN=Int(ceil(B/dx))
    Nbins=2*binN+1
    bins=zeros(Float64,Nbins)
    make_bins(bins,B,dx,binN)

    F=zeros(typeof(sigma_a),Nbins,Nbins)
    sigma_effective=dt*sigma_a^2   #here I transform std to variance. The energy landscape is symtric for pos and neg sitgma_a
    binN=convert(Int,(Nbins-1)/2)

    params2=zeros(typeof(sigma_a),7)
    params2[1]=s[1]
    params2[2]=c2[1]
    params2[3]=c4[1]
    params2[4]=c6[1]
    params2[5]=sigma_effective
    params2[6]=bias_d[1]
    params2[7]=hbias_d[1]*hb

    x=zeros(typeof(sigma_a),Nsteps*Nframes,Nbins)
    j=1
    absorbing=1
    Fmatrix2(F,params2,bins,absorbing)
    x0=zeros(typeof(sigma_a),Nbins)
    x0[Int( (Nbins+1)/2)]=1
    println("sizeF:",size(F), " size x0",size(x0))
    x[1,:]=F*x0

    for iframe in 1:Nframes

        params2[1]=s[iframe]
        params2[2]=c2[iframe]
        params2[3]=c4[iframe]
        params2[4]=c6[iframe]
        F=zeros(typeof(sigma_a),Nbins,Nbins)
        Fmatrix2(F,params2,bins,absorbing)
        if iframe==1
            N_aux=Nsteps-1
        else
            N_aux=Nsteps
        end
        for iN in 1:N_aux
            #println(sum(x[j,:]))
            x[j+1,:]=F*x[j,:]
            j+=1
        end

    end
    return x,bins
end

function simulation_Nframe_general_trace(param,stim,internal_noise,T,tau)
    #println("hola")
    Nstim,Nframes=size(stim)
    aux2=size(internal_noise)
    NT=aux2[2]
    NTframes=NT/Nframes
    deltat=T/NT
    dt_sim=deltat/tau
    dt_sqrt=sqrt(dt_sim)
    d=zeros(Float64,Nstim)
    x=zeros(Float64,Nstim,NT+1)
    c2=zeros(Nframes)
    c4=zeros(Nframes)
    c6=zeros(Nframes)
    println("hola ",dt_sim," ", deltat," ",T," ", NT)
    hbias_d=param["hbias_d"][1]
    bias_d=param["bias_d"][1]

    s=stim
    #println(hb,hbias_d)
    for itrial in 1:Nstim
        it=2
        x[itrial,1]=param["x0"][1]

        #println(s)
        for iframe in 1:Nframes
            for i in 1:NTframes

                # println("it: ",it-1, " itrial ", itrial," iframe ", iframe)
                x[itrial,it]=x[itrial,it-1]-(-(s[iframe])-c2[iframe]*x[itrial,it-1]+c4[iframe]*x[itrial,it-1]^3+c6[iframe]*x[itrial,it-1]^5)*dt_sim+dt_sqrt*internal_noise[itrial,it-1]
                it=it+1
            end
        end

        d[itrial]=sign(x[itrial,end])

    end
    return d,x
end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
# "volume" derivatives

function ∂x!(∇Q, Q, D, nrealelem, N, ξ1x1)
    for e in 1:nrealelem
        for i in 1:(N[1]+1)
            for j in 1:(N[2]+1)
                for k in 1:(N[3]+1)
                    for n in 1:(N[1]+1)
                        ∇Q[i,j,k,e,1] += ξ1x1[i,j,k,e] * D[i, n] * Q[n,j,k,e]
                    end
                end
            end
        end
    end
    return nothing
end

function ∂y!(∇Q, Q, D, nrealelem, N, ξ2x2)
    for e in 1:nrealelem
        for i in 1:(N[1]+1)
            for j in 1:(N[2]+1)
                for k in 1:(N[3]+1)
                    for n in 1:(N[1]+1)
                        ∇Q[i,j,k,e,2] += ξ2x2[i,j,k,e] * D[j, n] * Q[i,n,k,e]
                    end
                end
            end
        end
    end
    return nothing
end

function ∂z!(∇Q, Q, D, nrealelem, N, ξ3x3)
    for e in 1:nrealelem
        for i in 1:(N[1]+1)
            for j in 1:(N[2]+1)
                for k in 1:(N[3]+1)
                    for n in 1:(N[1]+1)
                        ∇Q[i,j,k,e,3] += ξ3x3[i,j,k,e] * D[k, n] * Q[i,j,n,e]
                    end
                end
            end
        end
    end
    return nothing
end

function ∇!(∇Q, Q, grid)
    N = polynomialorders(grid)
    nrealelem = Int(length(grid.vgeo[:, grid.x1id, : ,: ]) / prod(N .+ 1))
    r_∇Q = reshape(∇Q, ((N .+1)..., nrealelem, 3))
    r_Q  = reshape(Q, ((N .+1)..., nrealelem))
    r_∇Q .= -0.0
    r_ξ1x1  = reshape(grid.vgeo[:, grid.ξ1x1id, :, : ], ((N .+1)..., nrealelem))
    r_ξ2x2  = reshape(grid.vgeo[:, grid.ξ2x2id, :, : ], ((N .+1)..., nrealelem))
    r_ξ3x3  = reshape(grid.vgeo[:, grid.ξ3x3id, :, : ], ((N .+1)..., nrealelem))
    D = grid.D[1]
    ∂x!(r_∇Q, r_Q, D, nrealelem, N, r_ξ1x1)
    ∂y!(r_∇Q, r_Q, D, nrealelem, N, r_ξ2x2)
    ∂z!(r_∇Q, r_Q, D, nrealelem, N, r_ξ3x3)
    return nothing
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using CMAP
using Test
using DataFrames

@testset "CMAP.jl" begin
    status, response = query(API(), "select * from tblMakes")
    @test nrow(response) == 3
end

@testset "Gradients 2" begin
    @testset "Nutrients" begin
        @test get_dataset("tblMGL1704_Gradients2_Nutrients") isa DataFrame
    end
end