text:
module BoxChopper

export volfrac  # function

const N, P = 1, 2  # negative, positive directions
const NP = (N,P)  # tuple allocation is more efficient than array allocation

const X, Y, Z = 1, 2, 3  # x, y, z coordinates
const XYZ, YZX, ZXY = (X,Y,Z), (Y,Z,X), (Z,X,Y)
const CYC_AXES = (XYZ, YZX, ZXY)

typealias Tuple2{T} NTuple{2,T}
typealias Tuple3{T} NTuple{3,T}

function calc_vcbits(box::Tuple3{Tuple2{Real}}, nout::Tuple3{Real}, nr₀::Real)
    # Calculate the bit array that indicates vertex contained-ness.

    # The bit array is 8 bits (= 1 byte).  The kth bit is 1 if the kth verex is
    # contained in the half space defined by the plane, an 0 otherwise.  The 8 vertices
    # of the box are indexed 1 through 8 in the order of (---), (+--), (-+-), (++-),
    # (--+), (+-+), (-++), (+++), where, e.g., (+--) indicates the vertex at the
    # +x, -y, -z corner.
    nx, ny, nz = nout
    vcbits = 0x00
    bit = 0x01
    for sz = NP, sy = NP, sx = NP
        if nx*box[X][sx] + ny*box[Y][sy] + nz*box[Z][sz] ≤ nr₀
            vcbits |= bit
        end
        bit <<= 1
    end

    return vcbits
end

function isquadsect(vcbits::UInt8)
    # Determine if the vertex contained-ness bit array corresponds to the case where
    # the plane crosses one of the three sets of four parallel edges of the box.

    return vcbits==0x0F || vcbits==~0x0F || vcbits==0x33 || vcbits==~0x33 || vcbits==0x55 || vcbits==~0x55

    # Equivalent to
    # n = count_ones(vcbits)
    # m = count_ones(vcbits & 0x0F)
    # return n==4 && iseven(m)
end

function edgedir_quadsect(vcbits::UInt8)
    # For the cases where the plane crosses a set of four parallel edges of the
    # box, determine which direction those edges lie.

    if vcbits==0x0F || vcbits==~0x0F
        dir = Z
    elseif vcbits==0x33 || vcbits==~0x33
        dir = Y
    else
        assert(vcbits==0x55 || vcbits==~0x55)
        dir = X
    end

    return dir
end

function rvol_quadsect(box::Tuple3{Tuple2{Real}}, nout::Tuple3{Real}, nr₀::Real, vcbits::UInt8)
    # Return the volume fraction for the case where the plane croses a set of four
    # parallel edges.

    w = edgedir_quadsect(vcbits)
    ∆w = box[w][P] - box[w][N]

    ~, u, v = CYC_AXES[w]
    nw, nu, nv = nout[w], nout[u], nout[v]
    mean_cepts = 4nr₀
    for sv = NP, su = NP
        mean_cepts -= nu*box[u][su] + nv*box[v][sv]
    end
    mean_cepts /=  nw * 4∆w

    sw = nw>0 ? N:P  # nw cannot be 0
    return abs(mean_cepts - box[w][sw]/∆w)
end

function rvol_gensect(box::Tuple3{Tuple2{Real}}, nout::Tuple3{Real}, nr₀::Real, vcbits::UInt8)
    # Calculate the volume fraction of most general cases, by cutting out corners
    # from a triangular pyramid.
    # Assume count_ones(vcbits) ≤ 4.  Othewise, call this function with flipped
    # nout, nr₀, vcbits.

    nx, ny, nz = nout
    sx, sy, sz = (nx≥0 ? N:P), (ny≥0 ? N:P), (nz≥0 ? N:P)  # signs of corner
    cx, cy, cz = box[X][sx], box[Y][sy], box[Z][sz]  # corner coordinates
    ∆x, ∆y, ∆z = box[X][P]-box[X][N], box[Y][P]-box[Y][N], box[Z][P]-box[Z][N]  # box edges
    nxcx, nycy, nzcz = nx*cx, ny*cy, nz*cz

    rmax, rmid, rmin =  # (lengths from corner to intercetps) / (box edges)
        abs((nr₀-nycy-nzcz)/nx-cx)/∆x, abs((nr₀-nzcz-nxcx)/ny-cy)/∆y, abs((nr₀-nxcx-nycy)/nz-cz)/∆z

    # Sort rmax, rmin, rmin properly.
    if rmax < rmid; rmax, rmid = rmid, rmax; end
    if rmid < rmin; rmid, rmin = rmin, rmid; end
    if rmax < rmid; rmax, rmid = rmid, rmax; end

    # Calculate the volume of the triangular pyramid, and cut off appropriate corners.
    tmax = 1 - 1/rmax
    rvol_core = 1 + tmax + tmax^2
    if rmid > 1
        tmid = 1 - 1/rmid
        rvol_core -= rmax * tmid^3
    end
    if rmin > 1
        tmin = 1 - 1/rmin
        rvol_core -= rmax * tmin^3
    end

    return rvol_core * rmin * rmid / 6
end

"""
    volfrac(box, nout, r₀)

Returns the volume fraction `rvol = vol(box ⋂ half-space) / vol(box)` that indicates
how much portion of the volume of a box is included the half-space defined by a
plane.  The result is between 0.0 and 1.0, inclusive.

The box needs to be aligned with the Cartesian axes, and is described by `box =
[xn xp; yn yp; zn zp]` that specifies the x-, y-, z-boundaries of the box.  The
half-space is described by the boundary plane.  The boundary plane is described by
the outward normal vector `nout = [nx, ny, nz]` and a point `r₀ = [rx, ry, rz]`
on the plane.  `nout` does not have to be normalized.
"""
function volfrac(box::Tuple3{Tuple2{Real}}, nout::Tuple3{Real}, r₀::Tuple3{Real})
    # Return the volume fraction rvol = vol(box ⋂ half-space) / vol(box).

    nr₀ = nout[X]⋅r₀[X] + nout[Y]⋅r₀[Y] + nout[Z]⋅r₀[Z]
    vcbits = calc_vcbits(box, nout, nr₀)
    nvc = count_ones(vcbits)  # number of vertices contained

    if nvc == 8  # box is inside half-space
        rvol = 1.
    elseif nvc == 0  # box is outside half-space
        rvol = 0.
    elseif isquadsect(vcbits) # plane crosses a set of four parallel edges of box
        rvol = rvol_quadsect(box, nout, nr₀, vcbits)
    elseif nvc ≤ 4 # general cases with nvc ≤ 4
        rvol = rvol_gensect(box, nout, nr₀, vcbits)
    else  # general cases with nvc ≥ 5
        assert(nvc ≥ 5)
        rvol = 1. - rvol_gensect(box, (-nout[1],-nout[2],-nout[3]), -nr₀, ~vcbits)
    end

    return rvol
end

end  # module PlanarBoxDivider


text:
include("./exportData.jl")
using .exportData, DifferentialEquations, DiffEqParamEstim, PyPlot, Optim, DelimitedFiles

const N = Int(round(3/0.015 * sqrt(2)))
const D1= 0.0010/(0.015*0.015)
const D2= 0.0002/(0.015*0.015)

# define fenkar system
function H(x;k=100.0)
	return 0.5*(1.0 + tanh(k*x))
end
function fenkar!(dx, x, p, t)

	# parameters
	@views tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi,θ,t0,TI,IA = p[1:17]

	@inbounds for n=1:N

		# fenkar dynamics
		dx[n,1] =-(x[n,1]*H(uc-x[n,1])/to + H(x[n,1]-uc)/tr - x[n,2]*H(x[n,1]-uc)*(1.0-x[n,1])*(x[n,1]-uc)/td - x[n,3]*H(x[n,1]-ucsi;k=xk)/tsi)
		dx[n,2] = H(uc-x[n,1])*(1.0-x[n,2])/(tv1m*H(uv-x[n,1]) + tv2m*H(x[n,1]-uv)) - H(x[n,1]-uc)*x[n,2]/tvp
		dx[n,3] = H(uc-x[n,1])*(1.0-x[n,3])/twm - H(x[n,1]-uc)*x[n,3]/twp

		# no-flux with angular dependence on two diffusions
		if n==1
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(2*x[2,1] - 2*x[n,1])
		elseif n==N
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(2*x[N-1,1] - 2*x[n,1])
		else
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(x[n-1,1] + x[n+1,1] - 2*x[n,1])
		end
	end

	@inbounds for n=1:10
		dx[n,1] = dx[n,1] + H(t-t0;k=1.0)*IA*sin(pi*(t-t0)/TI)^500
	end

	return nothing
end

# fenkar system parameters in p[1:13] (BR; from https://doi.org/10.1063/1.166311)
# and then p[14:14] = [fiber_angle]
# and then p[15:17] = [offset_time, stimulus_period, stimulus_amplitude]
   p = [29.0, 19.6, 1250.0, 3.33, 41.0, 870.0, 0.25, 12.5, 33.3, 10.0, 0.13, 0.04, 0.85, 0.0,  1.0,   125.0, 0.1]

tspan = (0.0, 3000.0)
u0 = zeros(Float64,N,3)

#solve and take data
prob = ODEProblem(fenkar!, u0, tspan, p)
sol = solve(prob, Tsit5(); saveat=tspan[begin]:5.0:tspan[end])

x = collect(0:(N-1)).*0.015;
plt.pcolormesh(sol.t, x, sol[:,1,:], rasterized=true, snap=true, vmin=0.0, vmax=1.0)
plt.colorbar()
plt.savefig("./pde/optim/baseline.pdf")
plt.close()

# getExpData
ind=78
t, data, p[15] = getExpData(ind; tInds=1:1000)

# remake problem
prob = remake(prob, tspan=(0.0, t[end]), p=p)

# use the save_idxs to select the mid-point of the u-variable
sol = solve(prob, Tsit5(); saveat=t, save_idxs=Int(round(N//2)):Int(round(N//2)))

fig = plt.figure(figsize=(4,2))
plt.plot(t, data[:], "-", label="Exp. Sample")
plt.plot(sol.t, sol[1,:], "--", label="Initial Parameters")
plt.xlabel("\$ t \$")

# define the cost function
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t,data), maxiters=100000, verbose_opt=true, verbose_steps=5, saveat=t, save_idxs=Int(round(N//2)):Int(round(N//2)));

# define limits for parameters
#          tsi,  tv1m,   tv2m,   tvp,   twm,    twp,  td,    to,    tr,   xk,  uc,    uv, ucsi,     θ,   t0,    TI, IA
#   p = [ 33.4,  25.7, 1162.6,   4.00, 26.7,  662.8, 0.26, 18.7,  29.0, 10.0, 0.26, 0.13,  0.41,  0.0,  1.0,   BCL, 0.1]
lower = [  1.0,  10.0,  400.0,   2.33, 10.0,  100.0, 0.10, 10.0,  10.0,  9.0, 0.10, 0.01,  0.10, -0.3,  0.0, 100.0, 0.0]
upper = [ 50.0,  50.0, 2000.0,   4.33, 50.0, 1000.0, 0.50, 50.0,  50.0, 11.0, 0.90, 0.90,  0.90,  0.3, t[1], 150.0, 1.0]

@assert all(p .>= lower) && all(p .<= upper)

# fit
result = optimize(cost_function, lower, upper, p, Fminbox(NelderMead()))
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.minimizer) \n"
prob = remake(prob, p=result.minimizer)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=Int(round(N//2)):Int(round(N//2)))
plt.plot(sol.t, sol[1,:], "--", label="Nelder-Mead Parameters")

result = optimize(cost_function, lower, upper, result.minimizer, Fminbox(BFGS()))
@info "Initial p = $(p) \n"
@info "Minimizer = $(result.minimizer) \n"
prob = remake(prob, p=result.minimizer)
sol = solve(prob, Tsit5(); saveat=t, save_idxs=Int(round(N//2)):Int(round(N//2)))
plt.plot(sol.t, sol[1,:], "--", label="BFGS Parameters")

plt.savefig("./pde/optim/$(ind)_fittings.pdf")
plt.close()

open("./pde/optim/$(ind)_params.txt", "w") do io
	writedlm(io, ["tsi" "tv1m" "tv2m" "tvp" "twm" "twp" "td" "to" "tr" "xk" "uc" "uv" "ucsi" "θ" "t0" "TI" "IA"])
	writedlm(io, "\n")
	writedlm(io, transpose(p))
	writedlm(io, "\n")
	writedlm(io, transpose(result.minimizer))
end

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "rocminfo"
version = v"4.0.0"

# Collection of sources required to build
sources = [
    ArchiveSource("https://github.com/RadeonOpenCompute/rocminfo/archive/rocm-$(version).tar.gz",
                  "0b3d692959dd4bc2d1665ab3a838592fcd08d2b5e373593b9192ca369e2c4aa7"), # 4.0.0
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/rocminfo*/

mkdir build && cd build
cmake \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_PREFIX_PATH=${prefix} \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("rocminfo", :rocminfo),
    ExecutableProduct("rocm_agent_enumerator", :rocm_agent_enumerator),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("hsa_rocr_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil