text:
# BTE parameter setup

# Constants #################################################################
# These numbers were adapted from constants_module.f90 from the original code

# Fundamental
const dirac = 1.054571628e-34;
const boltzman = 1.3806503e-23;
const hobol = 7.63822401661014e-12; # dirac/boltzman

# Material
# NOTE: many of these have a 3-letter suffix like "_LAS"
#   1st: L = longitudinal polarization, T = transverse polarization
#   2nd: A = accoustic phonons, O = optical phonons
#   3rd: S = silicon, G = germanium
# For this test we are only considering TAS. Later LAS will be added.
const freq_min_TAS = 0;
const freq_max_TAS = 2.97927417405992e13;
const c_TAS = -2.26e-7; # c and vs are coefficients in the wave vector equation:
const vs_TAS= 5230;     # wk = wo + vs*|K| + c*|K|^2   where wo=0 for Transverse.

const freq_min_LAS = freq_max_TAS; # In the original code this was 0??
const freq_max_LAS = 7.728337675901222e13;
const c_LAS = -2.0e-7;
const vs_LAS= 9010;

# Others to be added as needed

# For relaxation time (non-Broido)
const wmax_half = 2.417e13;
const lattice_const = 5.43e-10;
const rho = 2.33e3;
const btu = 5.5e-18;
const bl = 2.0e-24;
const btn =9.3e-13;

# For Broido relaxation time
const debye_temp = 636.0;
const a_n_TAS = 10.9e-20;
const a_u_TAS = 37.8e-47;
const a_n_LAS = 7.10e-20;
const a_u_LAS = 9.51e-47;

# 5-point gaussian quadrature
g5xi = [-0.906179845938664, -0.538469310105683, 0.0, 0.538469310105683, 0.906179845938664]; # gaussian quadrature points
g5wi = [0.23692688505618908, 0.47862867049936647, 0.5688888888888889, 0.47862867049936647, 0.23692688505618908]; # gaussian weights

# 20-point gaussian quadrature
g20xi = zeros(20);
g20wi = zeros(20);
g20xi[1] = 0.076526521133497; g20wi[1] = 0.152753387130725;
g20xi[2] = 0.227785851141645; g20wi[2] = 0.149172986472603;
g20xi[3] = 0.373706088715419; g20wi[3] = 0.142096109318382;
g20xi[4] = 0.510867001950827; g20wi[4] = 0.131688638449176;
g20xi[5] = 0.636053680726515; g20wi[5] = 0.118194531961518;
g20xi[6] = 0.746331906460150; g20wi[6] = 0.101930119817240;
g20xi[7] = 0.839116971822218; g20wi[7] = 0.083276741576704;
g20xi[8] = 0.912234428251325; g20wi[8] = 0.062672048334109;
g20xi[9] = 0.963971927277913; g20wi[9] = 0.040601429800386;
g20xi[10] = 0.993128599185094; g20wi[10] = 0.017614007139152;
for i = 1:10
    g20xi[10+i] = -g20xi[i];
    g20wi[10+i] = g20wi[i];
end

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

# Direction vectors
# input: number of directions
# output: (x_i, y_i) direction vector components
function get_directions_2d(n)
    dir_x = zeros(n);
    dir_y = zeros(n);
    reflect = zeros(Int,n); # The index of the reflecting direction
    dphi = 2*pi/n;
    phione = dphi/2;
    for di=1:n
        # dir_x[di] = cos(2*pi*(di-1)/n);
        # dir_y[di] = sin(2*pi*(di-1)/n);

        # modified to match original
        thisphi = phione + (di-1)*dphi;
        dir_x[di] = sin(thisphi);
        dir_y[di] = cos(thisphi);
    end
    if mod(n,2) == 0
        halfn = Int(n/2);
        for di=1:halfn
            reflect[di] = n-di+1;
            reflect[n-di+1] = di;
        end
    else
        println("Please use an even number of directions for reflection purposes. Reflections will be incorrect.")
        halfn = Int(floor(n/2));
        for di=1:halfn
            reflect[di] = n-di+1;
            reflect[n-di+1] = di;
        end
        reflect[halfn+1] = 1;
    end

    return (dir_x, dir_y, reflect);
end

# Frequency bands
# input: number of bands
# output: centers of frequency bands, width of bands
function get_band_frequencies(n, low=0, high=0; polarization="T")
    if polarization=="T"
        wmin = freq_min_TAS;
        wmax = freq_max_TAS;
    else # L
        wmin = freq_min_LAS;
        wmax = freq_max_LAS;
    end

    if low==0
        low = 1;
    end
    if high==0
        high = n;
    end

    dw = (wmax-wmin)/n;
    centers = zeros(high-low+1);
    for i=1:length(centers)
        centers[i] = wmin + (i - 1)*dw + (low-1)*dw + dw/2;
    end
    return (centers, dw);
end

# Group speed
# input: band center frequencies
# output: group speeds for each band
function get_group_speeds(freq; polarization="T")
    if polarization=="T"
        vs = vs_TAS;
        c = c_TAS;
    else # L
        vs = vs_LAS;
        c = c_LAS;
    end

    n = length(freq);
    gv = zeros(n);
    for i=1:n
        gv[i] = sqrt(vs*vs + 4 * c * freq[i]);
    end
    return gv;
end

# Scattering time scale from relaxtime.f90 (NOT Broido version)
# input: beta array to update, band center frequencies, temperature
# output: band average inverse time scales (beta = 1/tau)
function get_time_scale!(beta::Array, freq::Array, temp::Array; polarization="T")
    n = length(temp); # number of cells
    m = length(freq); # number of bands
    for j=1:n # loop over cells
        for i=1:m # loop over bands
            if polarization == "T"
                beta[i, j] = btn * (temp[j]^4) * freq[i]
                if freq[i] >= wmax_half
                    beta[i, j] += btu * (freq[i]*freq[i]) / sinh(hobol * freq[i] / temp[j]);
                end
            else
                beta[i, j] = bl * (temp[j]^3) * (freq[i]*freq[i]);
            end
        end
    end
end
# Similar to above, but for a single band and temp
function get_time_scale(freq::Number, temp::Number; polarization="T")
    if polarization == "T"
        beta = btn * (temp^4) * freq
        if freq >= wmax_half
            beta += btu * (freq*freq) / sinh(hobol * freq / temp);
        end
    else
        beta = bl * (temp^3) * (freq*freq);
    end

    return beta;
end

# Equilibrium intensity
# Io = integral_dw(dirac*freq/(8*pi^3) * 1/(a) * (b)^2 * dw)
#   where a = exp(hobol*freq/temp) - 1
#         b = (-vs + sqrt(vs^2 + 4*freq*c)) / (2*c)
#
# input: Intensity array to update, band center frequency, band width, temperature of each cell
# output: equilibrium intensity in each cell
function equilibrium_intensity!(intensity::Array, freq::Array, dw, temp::Array; polarization="T")
    if polarization=="T"
        vs = vs_TAS;
        c = c_TAS;
        extra_factor = 2;
    else # L
        vs = vs_LAS;
        c = c_LAS;
        extra_factor = 1;
    end

    n = length(temp); # should be number of cells
    # dirac/(32*pi^3) = 1.062861036647414e-37
    const_part = 1.062861036647414e-37 / (c*c) * dw/2; # constants to pull out of integral
    for ci=1:n # loop over cells
        for i=1:length(freq) # loop over bands
            tmp = 0.0;
            for gi=1:20
                fi = freq[i] + dw/2 * g20xi[gi]; # frequency at gauss point
                tmp += (fi * (-vs + sqrt(vs*vs + 4*fi*c))^2 / (exp(hobol*fi/temp[ci]) - 1)) * g20wi[gi] * extra_factor;
            end
            intensity[i,ci] = tmp * const_part;
        end
    end
end
# Similar to above, but for a single frequency band and temperature
function equilibrium_intensity(freq::Number, dw, temp::Number; polarization="T")
    if polarization=="T"
        vs = vs_TAS;
        c = c_TAS;
        extra_factor = 2;
    else # L
        vs = vs_LAS;
        c = c_LAS;
        extra_factor = 1;
    end

    # dirac/(32*pi^3) = 1.062861036647414e-37
    const_part = 1.062861036647414e-37 * dw/2 / (c*c); # constants to pull out of integral
    intensity = 0.0;
    for gi=1:20
        fi = freq + dw/2 * g20xi[gi]; # frequency at gauss point
        K2 = (-vs + sqrt(vs*vs + 4*fi*c))^2; # K^2 * (2*c)^2   the (2*c)^2 is put in the const_part
        intensity += (fi * K2 / (exp(hobol*fi/temp) - 1)) * g20wi[gi] * extra_factor;
    end
    intensity *= const_part;

    return intensity;
end

# The integrated intensity in each band in each cell (integrated over directions)
function get_integrated_intensity!(int_intensity, intensity, ndirs, nbands)
    n = size(intensity,2); # num cells
    omega = 2*pi/ndirs; # angle per direction
    # Integrate over local bands and cells
    for i=1:n
        for j=1:nbands
            int_intensity[j,i] = 0.0;
            for k=1:ndirs
                int_intensity[j,i] += intensity[(j-1)*ndirs+k, i] * omega;
            end
        end
    end
end

# Derivative of equilibrium intensity with respect to temperature
# input: frequency bands, bandwidth, temperature at each cell
# output: dI/dT at each cell
# NOTE: This produces an array of dIdT for every cell and band. Single cell/band version is below
function dIdT(freq, dw, temp; polarization="T")
    if polarization=="T"
        vs = vs_TAS;
        c = c_TAS;
    else # L
        vs = vs_LAS;
        c = c_LAS;
    end

    n = length(temp); # number of cells
    m = length(freq); # number of bands
    didt = zeros(m, n);
    const_part = dirac * hobol / (8*pi^3);

    for i=1:n # loop over cells
        for j=1:m # loop over bands
            tmp = 0.0;
            for gi=1:29 # gaussian quadrature
                fi = freq[j] + dw/2 * g20xi[gi]; # frequency at gauss point
                # K2 = ((-vs + sqrt(vs*vs + 4*fi*c)) / (2*c))^2; # K^2
                tmpK = (-vs + sqrt(vs*vs + 4*fi*c)) / (2*c); # K updated to match ipcalc
                tmp2 = exp(hobol*fi/temp[i]);
                # tmp += (tmp2 * fi * K2 / (tmp2 - 1)^2) * wi[gi];
                tmp += (tmp2 * (fi * tmpK)^2 / (tmp2 - 1)^2) * g20wi[gi]; # updated to match ipcalc
            end
            didt[j, i] = const_part * dw * tmp / temp[i]^2;
        end
    end

    return didt;
end

# Derivative of equilibrium intensity with respect to temperature
# input: center frequency, bandwidth, temperature at one cell
# output: dI/dT for one cell and band
# NOTE: This is for one cell and band. Full cells/bands version above.
function dIdT_single(freq, dw, temp; polarization="T")
    if polarization=="T"
        vs = vs_TAS;
        c = c_TAS;
    else # L
        vs = vs_LAS;
        c = c_LAS;
    end

    tmp = 0.0;
    for gi=1:20 # gaussian quadrature
        fi = freq + dw/2 * g20xi[gi]; # frequency at gauss point
        # K2 = ((-vs + sqrt(vs*vs + 4*fi*c)) / (2*c))^2; # K^2
        tmpK = (-vs + sqrt(vs*vs + 4*fi*c)) / (2*c); # K updated to match ipcalc
        tmp2 = exp(hobol*fi/temp);
        # tmp += (tmp2 * fi * K2 / (tmp2 - 1)^2) * wi[gi];
        tmp += (tmp2 * (fi * tmpK)^2 / (tmp2 - 1)^2) * g20wi[gi]; # updated to match ipcalc
    end
    didt = tmp * dw * 3.2473482785757725e-48 / (temp * temp); # dirac * hobol / (8*pi^3) = 3.2473482785757725e-48

    return didt;
end

# Change temperature for one time step
# input: temperature array to update, previous step temperature, intensity from this and the previous time step,
#        temperature from previous time step, frequency bands, bandwidth
# output: temperature for next step for each cell
function get_next_temp!(temp_next, temp_last, I_next, freq, dw; polarization="T")
    debug = false;
    n = length(temp_last); # number of cells
    m = length(freq); # number of bands

    # G_last = get_integrated_intensity(I_last, ndirs, nbands);
    # G_next = get_integrated_intensity(I_next, ndirs, nbands);
    get_integrated_intensity!(G_next.values, I_next, ndirs, nbands);
    # didt = dIdT(freq, dw, temp_last, polarization=polarization); # Use single version in loop instead

    dt = Finch.time_stepper.dt;
    idt = 1/dt;

    if debug
        max_iters = 0;
        ave_iters = 0;
    end

    tol = 1e-10;
    maxiters = 50;
    for i=1:n # loop over cells
        # old values are not updated
        uold = 0.0; # These strange names are taken from the fortran code
        gnb = 0.0;  #
        for j=1:m # loop over bands
            uold += Io.values[j,i] * idt / group_v[j];
            gnb += G_last.values[j,i] * idt / group_v[j];
            G_last.values[j,i] = G_next.values[j,i]; # That was the only place we need G_last, so update it here.
        end
        uold = 4*pi*uold;

        delta_T = 0;
        for iter=1:maxiters # iteratively refine delta_T
            unew = 0.0; # These strange names are taken from the fortran code
            gna = 0.0;  #
            uprime = 0.0;
            for j=1:m # loop over bands
                beta = 1 / get_time_scale(freq[j], temp_next[i], polarization=polarization);
                didt = dIdT_single(freq[j], dw, temp_last[i], polarization=polarization);

                unew += equilibrium_intensity(freq[j], dw, temp_next[i], polarization=polarization) * (beta + idt) / group_v[j];
                gna += G_next.values[j,i] * (beta + idt) / group_v[j];
                uprime += 4*pi * didt * (beta + idt) / group_v[j];
            end
            unew = 4*pi*unew;

            delta_T = (uold + gna - gnb - unew) / uprime;

            if debug println("cell "*string(i)*" ("*string(uold)*", "*string(unew)*", "*string(gna)*", "*string(gnb)*", "*string(uprime)*") : "*string(delta_T)) end

            temp_next[i] = temp_next[i] + delta_T;

            if abs(delta_T) < tol
                if debug
                    max_iters = max(max_iters, iter);
                    ave_iters += iter;
                end
                break;
            else
                if iter==maxiters
                    println("Temperature change didn't converge. Last delta_T = "*string(delta_T));
                end
            end
        end# iterative refinement
    end# cell loop
    if debug
        println("ave iterations: "*string(ave_iters/n)*"  max: "*string(max_iters))
    end

    return temp_next;
end

# Band-based parallel version of above
# Change temperature for one time step
# input: temperature array to update, previous step temperature, intensity from this and the previous time step,
#        temperature from previous time step, frequency bands, bandwidth
# output: temperature for next step for each cell
function get_next_temp_par!(temp_next, temp_last, I_next, freq, dw; polarization="T")
    debug = false;
    n = length(temp_last); # number of cells
    m = length(freq); # number of bands

    get_integrated_intensity!(G_next.values, I_next, ndirs, nbands);
    # didt = dIdT(freq, dw, temp_last, polarization=polarization); # Use single version in loop instead

    dt = Finch.time_stepper.dt;
    idt = 1/dt;

    if debug
        max_iters = 0;
        ave_iters = 0;
    end

    tol = 1e-10;
    maxiters = 50;

    # First find uold and gnb from the previous step
    cell_comm_rank = MPI.Comm_rank(cell_comm);
    offset = cell_comm_rank * n;
    for i=1:n # loop over cells
        uold[offset + i] = 0.0;
        gnb[offset + i] = 0.0;
        for j=1:m # loop over local bands
            uold[offset + i] += Io.values[j,i] * idt / group_v[j];
            gnb[offset + i] += G_last.values[j,i] * idt / group_v[j];
            G_last.values[j,i] = G_next.values[j,i]; # That was the only place we need G_last, so update it here.
        end
        uold[offset + i] = 4*pi*uold[offset + i];
    end

    # Allgather across this mesh communicator to combine bands
    p_data1 = MPI.UBuffer(uold, n, num_band_partitions, MPI.Datatype(Float64));
    MPI.Allgather!(p_data1, cell_comm);
    p_data2 = MPI.UBuffer(gnb, n, num_band_partitions, MPI.Datatype(Float64));
    MPI.Allgather!(p_data2, cell_comm);

    # reduce
    for i=1:n # loop over cells
        for j=2:num_band_partitions # loop over band partitions
            uold[i] += uold[i + (j-1)*n];
            gnb[i] += gnb[i + (j-1)*n];
        end
    end

    # Then iteratively refine delta_T
    converged = fill(false, n);
    num_not_converged = n;
    for iter=1:maxiters
        for i=1:n # loop over cells
            if converged[i]
                continue;
            end

            unew[offset + i] = 0.0;
            gna[offset + i] = 0.0;
            uprime[offset + i] = 0.0;

            for j=1:m # loop over local bands
                beta = 1 / get_time_scale(freq[j], temp_next[i], polarization=polarization);
                didt = dIdT_single(freq[j], dw, temp_last[i], polarization=polarization);

                unew[offset + i] += equilibrium_intensity(freq[j], dw, temp_next[i], polarization=polarization) * (beta + idt) / group_v[j];
                gna[offset + i] += G_next.values[j,i] * (beta + idt) / group_v[j];
                uprime[offset + i] += 4*pi * didt * (beta + idt) / group_v[j];
            end
            unew[offset + i] = 4*pi*unew[offset + i];
        end

        # Allgather across this mesh communicator to combine bands
        p_data1 = MPI.UBuffer(unew, n, num_band_partitions, MPI.Datatype(Float64));
        MPI.Allgather!(p_data1, cell_comm);
        p_data2 = MPI.UBuffer(gna, n, num_band_partitions, MPI.Datatype(Float64));
        MPI.Allgather!(p_data2, cell_comm);
        p_data3 = MPI.UBuffer(uprime, n, num_band_partitions, MPI.Datatype(Float64));
        MPI.Allgather!(p_data3, cell_comm);

        # reduce
        for i=1:n # loop over cells
            if converged[i]
                continue;
            end
            for j=2:num_band_partitions # loop over band partitions
                unew[i] += unew[i + (j-1)*n];
                gna[i] += gna[i + (j-1)*n];
                uprime[i] += uprime[i + (j-1)*n];
            end

            delta_T = (uold[i] + gna[i] - gnb[i] - unew[i]) / uprime[i];

            if debug println("cell "*string(i)*" ("*string(uold[i])*", "*string(unew[i])*", "*string(gna[i])*", "*string(gnb[i])*", "*string(uprime[i])*") : "*string(delta_T)) end

            temp_next[i] = temp_next[i] + delta_T;

            if abs(delta_T) < tol
                if debug
                    max_iters = max(max_iters, iter);
                    ave_iters += iter;
                end
                converged[i] = true;
                num_not_converged = num_not_converged - 1;
            else
                if iter==maxiters
                    println("Temperature change didn't converge. Last delta_T = "*string(delta_T));
                end
            end
        end# cell loop
        if num_not_converged == 0
            break;
        end
    end# iterative refinement
    if debug
        println("ave iterations: "*string(ave_iters/n)*"  max: "*string(max_iters))
    end

    return temp_next;
end

# Update temperature, equilibrium I, and time scale
#
function update_temperature(temp, I_next, freq, dw, band_parallel=false)

    temp_last = deepcopy(temp);

    if band_parallel
        get_next_temp_par!(temp, temp_last, I_next, freq, dw);
    else
        get_next_temp!(temp, temp_last, I_next, freq, dw);
    end

    equilibrium_intensity!(Io.values, freq, dw, temp);
    get_time_scale!(beta.values, freq, temp);

    # check for problems
    for i=1:length(I_next)
        if I_next[i] === NaN
            println("Error: NaN found in I");
            exit(0);
        end
    end
end

# # Temporary until I figure out where to do this automatically
# function communicate_bands()
#     rank = Finch.config.proc_rank;
#     np = Finch.config.num_procs;
#     nel = Finch.mesh_data.nel;

#     # This is my band range
#     band_range = (rank*40/np + 1):((rank+1)*40/np);
#     # Part of Intensity that I did
#     my_count = Int(nel*16*40/np);
#     my_intensity = zeros(my_count);
#     next = 1;
#     for i=1:nel
#         for k=1:band_range
#             for j=1:16
#                 m_intensity[next] = I.values[(k-1)*16 + j, i];
#                 next += 1;
#             end
#         end
#     end

#     # allgather the intensity
#     p_data = zeros(length(I.values));
#     p_data_in = Finch.MPI.UBuffer(p_data, my_count, config.num_procs, MPI.Datatype(Float64));
#     p_data_out = my_intensity;
#     Finch.MPI.Allgather!(p_data_out, p_data_in, MPI.COMM_WORLD);

#     # reorganize it
#     next = 1;
#     for i=1:nel
#         for k=1:40
#             for j=1:16
#                 I.values[(k-1)*16 + j, i] = p_data_in[next];
#                 next += 1;
#             end
#         end
#     end

# end

#############################################################################
## alternative models.

# # Scattering time scale using Broido data from relaxtime_broido.f90
# # input: tau array to update, band center frequencies, temperature
# # output: band average time scales
# function get_time_scale!(tau::Array, freq::Array, temp::Array; polarization="T")
#     if polarization=="T"
#         a_n = a_n_TAS;
#         a_u = a_u_TAS;
#     else # L
#         a_n = a_n_LAS;
#         a_u = a_u_LAS;
#     end

#     n = length(temp); # number of cells
#     m = length(freq); # number of bands
#     for j=1:n # loop over cells
#         for i=1:m # loop over bands
#             tmp = temp[j]*(1 - exp(-3*temp[j]/debye_temp));
#             tau[i, j] = 1 / (a_n * freq[i]^2 * tmp + a_u * freq[i]^4 * tmp);
#         end
#     end
# end
# # Similar to above, but for a single band and temp
# function get_time_scale(freq::Number, temp::Number; polarization="T")
#     if polarization=="T"
#         a_n = a_n_TAS;
#         a_u = a_u_TAS;
#     else # L
#         a_n = a_n_LAS;
#         a_u = a_u_LAS;
#     end

#     tmp = temp*(1 - exp(-3*temp/debye_temp));
#     tau = 1 / (a_n * freq^2 * tmp + a_u * freq^4 * tmp);

#     return tau;
# end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
### Implementation of a nodal polynomial basis for elements in 1 and 2 dimensions
#
# Polynomials are written in coordinates on reference elements.
###

# 1D nodal basis

using Polynomials

# Chebyshev nodes
nodes1(p) = reverse([cos((2i-1)/(2p)*pi) for i = 1:p])
nodes(p::DG1D) = nodes1(nbf(p))
#nodes(p) = [-1,-1/3,1/3,1]

phi1(nodes::Vector) = (x = Poly([0.0,1.0]); Poly{Float64}[Poly([1.0])*prod([k == i ? 1 : (x-nodes[i])/(nodes[k]-nodes[i]) for i=1:length(nodes)]) for k=1:length(nodes)])
phi(p::DG1D) = phi1(nodes(p))
dphi(p::DG1D) = polyder(phi(p))

# 2D nodal basis

#
# Consider the folling sketch of a reference element. We label
# reference coordinates as ξ and η (not x and y as those refer to real
# coordinates in the mesh)
#
#
#               η
#               ▲
#               |
#             1 -      v3
#               |      |\
#               |      |  \
#               |      |    \
#               |      |      \
#               |      |        \
#               |      |          \
#             0 -      |____________\
#               |      v1           v2
#               |
#               |
#               +------|-------------|----▶ ξ
#                      0             1
#
#
#
# Note the choice of counterclockwise orientation when numbering
# The vertices.
#

## p = 0 ##
#
# One node and (1/3,1/3) with constant basis functions
#
###########

p0ϕ1(ξ,η) = 1.
∇p0ϕ1{T}(ξ::T,η::T) = Vector2{T}(0.,0.)
const P0 = [p0ϕ1]
const DP0 = [∇p0ϕ1]
const N0 = Vector2{Float64}[Vector2{Float64}(1//3,1//3)]

## p = 1 ##
#
# Three nodes at the vertices of the reference element with linear
# basis functions
#
###########

p1ϕ1(ξ,η) = 1 - ξ - η
p1ϕ2(ξ,η) = ξ
p1ϕ3(ξ,η) = η

∇p1ϕ1{T}(ξ::T,η::T) = Vector2{T}(-1,-1)
∇p1ϕ2{T}(ξ::T,η::T) = Vector2{T}(1,0)
∇p1ϕ3{T}(ξ::T,η::T) = Vector2{T}(0,1)

const P1 = [p1ϕ1,p1ϕ2,p1ϕ3]
const DP1 = [∇p1ϕ1,∇p1ϕ2,∇p1ϕ3]
const N1 = Vector2{Float64}[Vector2{Float64}(0,0),Vector2{Float64}(1,0),Vector2{Float64}(0,1)]

## p = 2 ##
#
# Six nodes at the vertices + midpoints of the reference element with
# quadratic functions
#
###########

p2ϕ1(ξ,η) = 2.0*(0.5 - ξ - η)*(1.0 - ξ - η)
p2ϕ2(ξ,η) = 2.0ξ*(ξ - 0.5)
p2ϕ3(ξ,η) = 2.0η*(η - 0.5)
p2ϕ4(ξ,η) = 4.0ξ*(1. - ξ - η)
p2ϕ5(ξ,η) = 4.0η*ξ
p2ϕ6(ξ,η) = 4.0η*(1. - ξ - η)

# TODO: Make Vector2?
∇p2ϕ1{T}(ξ::T,η::T) = Vector2{T}(4ξ+4η-3,4ξ+4η-3)
∇p2ϕ2{T}(ξ::T,η::T) = Vector2{T}(4ξ - 1,0)
∇p2ϕ3{T}(ξ::T,η::T) = Vector2{T}(0,4η - 1)
∇p2ϕ4{T}(ξ::T,η::T) = Vector2{T}(-8ξ-4η+4,-4ξ)
∇p2ϕ5{T}(ξ::T,η::T) = Vector2{T}(4η,4ξ)
∇p2ϕ6{T}(ξ::T,η::T) = Vector2{T}(-4η,-8η-4ξ+4)

const P2 = [p2ϕ1,p2ϕ2,p2ϕ3,p2ϕ4,p2ϕ5,p2ϕ6]
const DP2 = [∇p2ϕ1,∇p2ϕ2,∇p2ϕ3,∇p2ϕ4,∇p2ϕ5,∇p2ϕ6]
const N2 = Vector2{Float64}[Vector2{Float64}(0.,0.),Vector2{Float64}(1.,0.),Vector2{Float64}(0.,1.),
            Vector2{Float64}(0.5,0),Vector2{Float64}(0.5,0.5),Vector2{Float64}(0,0.5)]

### End of P=2 Basis functions

for f in [P0,P1,P2,DP0,DP1,DP2]
    @eval $(f.env.name)(v::Vertex2) = $(f.env.name)(v.coords.e1,v.coords.e2)
    @eval $(f.env.name)(v::Vector2{Float64}) = $(f.env.name)(v.e1,v.e2)
end

getPhi(p) = (@assert 0 <= p <= 2; p==0 ? P0 : p == 1 ? P1 : P2)
getDPhi(p) = (@assert 0 <= p <= 2; p==0 ? DP0 : p == 1 ? DP1 : DP2)
nodes2(p) = (@assert 0 <= p <= 2; p==0 ? N0 : p == 1 ? N1 : N2)
phi(p::Galerkin2D) = getPhi(porder(p))
dphi(p::Galerkin2D) = getDPhi(porder(p))
nodes(p::Galerkin2D) = nodes2(porder(p))

# Evaluating coefficients on a basis
function evaluate_ref{N,T}(coeffs::Coeffs{N,T},basis)
    order(coeffs) == order(basis) || error("Basis and vector must agree")
    sum([a*f for (f,a) in zip(coeffs,basis)])
end

order(x::Array) = length(x)

function evaluate_ref(coeffs,basis,p)
    order(coeffs) == order(basis) || error("Basis and vector must agree (got $basis and $coeffs)")
    result = coeffs[1]*fapply(basis[1],p)
    for i = 2:length(basis)
        result += coeffs[i]*fapply(basis[i],p)
    end
    result
end

# Below are the hottest functions in the system. We unroll them manually
# For maximum performce. This is usually a bad idea, but because they are
# so important, we really have no other choice.

macro unroll_eval(call, coefff, basis)
    coefff = eval(coefff)
    funcs = eval(basis)
    adds = Expr(:block)
    for i = 2:length(funcs)
        addend = coefff(i,:($(Base.function_name(funcs[i]))(p)))
        push!(adds.args,:(@inbounds result += $addend))
    end
    Expr(:function,esc(call),quote
        @inbounds result = $(coefff(1,:($(Base.function_name(funcs[1]))(p))))
        $adds
        result
    end)
end

evalrefcoefff = (k,e)->:(coeffs[$k]*$e)
@unroll_eval evaluate_ref2dp0(coeffs,p) evalrefcoefff P0
@unroll_eval evaluate_ref2dp1(coeffs,p) evalrefcoefff P1
@unroll_eval evaluate_ref2dp2(coeffs,p) evalrefcoefff P2

cidevalrefcoefff = (k,e)->:(coeffs[1,cid,$k]*$e)
@unroll_eval evaluate_ref2dp0(coeffs,cid,p) cidevalrefcoefff P0
@unroll_eval evaluate_ref2dp1(coeffs,cid,p) cidevalrefcoefff P1
@unroll_eval evaluate_ref2dp2(coeffs,cid,p) cidevalrefcoefff P2

function evaluate_ref2d{T}(p,coeffs::Array{T},a,b)
    po = porder(p)
    if po == 0
        return evaluate_ref2dp0(coeffs,a,b)
    elseif po == 1
        return evaluate_ref2dp1(coeffs,a,b)
    elseif po == 2
        return evaluate_ref2dp2(coeffs,a,b)
    else
        error("Not specialized for this case")
    end
end
function evaluate_ref2d{T}(p,coeffs::Array{T},a)
    po = porder(p)
    if po == 0
        return evaluate_ref2dp0(coeffs,a)
    elseif po == 1
        return evaluate_ref2dp1(coeffs,a)
    elseif po == 2
        return evaluate_ref2dp2(coeffs,a)
    else
        error("Not specialized for this case")
    end
end
evaluate_ref2d(p,coeffs::Array{Coeffs},point) = evaluate_ref2d(p,coeffs.coeffs,point)

macro unroll_map(call, basis)
    funcs = eval(basis)
    cats = Expr(:vcat)
    for i = 1:length(funcs)
        push!(cats.args,:($(Base.function_name(funcs[i]))(p)))
    end
    Expr(:function,esc(call),quote
        $cats
    end)
end

@unroll_map evaluate_basis2dp0(p) P0
@unroll_map evaluate_basis2dp1(p) P1
@unroll_map evaluate_basis2dp2(p) P2

function evaluate_basis2d(p,point)
    po = porder(p)
    if po == 0
        return evaluate_basis2dp0(point)
    elseif po == 1
        return evaluate_basis2dp1(point)
    elseif po == 2
        return evaluate_basis2dp2(point)
    else
        error("Not specialized for this case")
    end
end

evalbasiscoefff = (k,e)->:(vRHS[$k] += mul*$e)
@unroll_eval evaluate_basis2dp0!(vRHS,mul,p) evalbasiscoefff P0
@unroll_eval evaluate_basis2dp1!(vRHS,mul,p) evalbasiscoefff P1
@unroll_eval evaluate_basis2dp2!(vRHS,mul,p) evalbasiscoefff P2

function evaluate_basis2d!(vRHS,p,mul,point)
    po = porder(p)
    if po == 0
        return evaluate_basis2dp0!(vRHS,mul,point)
    elseif po == 1
        return evaluate_basis2dp1!(vRHS,mul,point)
    elseif po == 2
        return evaluate_basis2dp2!(vRHS,mul,point)
    else
        error("Not specialized for this case")
    end
end

const evaluate = evaluate_ref

ext: jl


text:
"""
    geneMutation(genotype::GEPGenotype,
                 gpExperimentInfo::GEPInfo,
                 rng::Random.AbstractRNG)

documentation
"""
function geneMutation(genotype::GEPGenotype,
                      gpExperimentInfo::GEPInfo,
                      rng::Random.AbstractRNG)

    representation = genotype._representation
    functionSet = gpExperimentInfo._functionSet
    terminalSet = gpExperimentInfo._terminalSet
    headSize = gpExperimentInfo._head
    tailSize = gpExperimentInfo._tail
    lenRep = headSize + tailSize
    prob = lenRep < 5 ? 1/lenRep : 4/lenRep
    probTer = length(terminalSet) / (length(terminalSet) + length(functionSet))

    mutatedGenotype = Array{Node}(undef, lenRep)

    for i=1:headSize
        if (rand(rng) - (lenRep-i+1)*prob/(2*lenRep)) < prob
            if rand(rng) < probTer
                mutatedGenotype[i] = rand(rng, terminalSet)
            else
                mutatedGenotype[i] = rand(rng, functionSet)
            end
        else
            mutatedGenotype[i] = representation[i]
        end
    end

    for i=headSize+1:lenRep
        if (rand(rng) - (lenRep-i+1)*prob/(2*lenRep)) < prob
            mutatedGenotype[i] = rand(rng, terminalSet)
        else
            mutatedGenotype[i] = representation[i]
        end
    end

    return GEPGenotype(mutatedGenotype)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
import StatsBase.StatisticalModel
import Base.rand
using Optim

import Distributions: loglikelihood, dim

export coef!, coef_jac!, rand, loglikelihood, loglikelihood_jac, dim, mle

coef!(model::StatisticalModel,parameter)=error("Method coef!(model::$(typeof(model)),parameter) must be implemented for $(typeof(model)).")
coef_jac!(model::StatisticalModel,parameter)=error("Method coef_jac!(model::$(typeof(model)),parameter) must be implemented for $(typeof(model)).")
rand(model::StatisticalModel,T::Int)=error("Method rand(model::$(typeof(model)),T::Int) must be implemented for $(typeof(model)).")
loglikelihood(model::StatisticalModel,data)=error("Method loglikelihood(model::$(typeof(model)),data) must be implemented for $(typeof(model)).")
loglikelihood_jac(model::StatisticalModel,data)=error("Method loglikelihood_jac(model::$(typeof(model)),data) must be implemented for $(typeof(model)).")

#return the dimension of the parameter: used to provide a random starting value when calling mle()
dim(model::StatisticalModel)=error("If you want to use the default MLE method on a $(typeof(model)), you must define dim(model::$(typeof(model))).")

function loglikelihood(model::StatisticalModel,data,parameter)
	coef!(model,parameter)
	loglikelihood(model,data)
end

function loglikelihood_jac(model::StatisticalModel,data,parameter)
	coef_jac!(model,parameter)
	loglikelihood_jac(model,data)
end

function mle_nojac(model::StatisticalModel,data,thetai::Array{Float64,1}=rand(dim(model)),L=15000)
	fcalls=0
	function ff(theta)
		fcalls+=1
		llk=-loglikelihood(model,data,theta)
		println(fcalls,": ",round(theta,3),", ",round(llk,3))
		llk
	end
	df=Optim.DifferentiableFunction(ff)
	ret = Optim.optimize(df, thetai,method=:cg,iterations=L)
	println("  no jacobian, $(ret.iterations) iterations, $fcalls ff evaluations, final log-likelihood: $(round(-ret.f_minimum,4))")
	ret.minimum
end

function mle_jac(model::StatisticalModel,data,thetai::Array{Float64,1}=rand(dim(model)),L=15000)
	fcalls=0
	fjcalls=0
	ffjcalls=0
	function ff(theta)
		fcalls+=1
		-loglikelihood(model,data,theta)
	end
	function fj!(theta,jac)
		fjcalls+=1
		jac[:]=-loglikelihood_jac(model,data,theta)[2]
	end
	function ffj!(theta,jac)
		ffjcalls+=1
		res=loglikelihood_jac(model,data,theta)
		jac[:]=-res[2]
		-res[1]
	end
	df=Optim.DifferentiableFunction(ff,fj!,ffj!)
	ret = Optim.optimize(df, thetai,method=:cg,iterations=L)
	println("  with jacobian, $(ret.iterations) iterations, ($fcalls,$fjcalls,$ffjcalls) (ff,fj,ffj) evaluations, final log-likelihood: $(round(-ret.f_minimum,4))")
	ret.minimum
end

function mle(model::StatisticalModel,data,thetai::Array{Float64,1}=rand(dim(model)),L=15000)
	println()
	println(" generic optimization of the likelihood...")
	try
		mle_jac(model,data,thetai,L)
	catch
		mle_nojac(model,data,thetai,L)
	end
end

function fit!(model::StatisticalModel,data)
	thetahat=mle(model,data)
	coef!(model,thetahat)
end


text:
# Synchronization Functions

export barrier

const cl_mem_fence_flags = Cuint
const CLK_LOCAL_MEM_FENCE = cl_mem_fence_flags(1)
const CLK_GLOBAL_MEM_FENCE = cl_mem_fence_flags(2)

#barrier(flags=0) = @builtin_ccall("barrier", Cvoid, (Cuint,), flags)
barrier(flags=0) = Base.llvmcall(("""
        declare void @_Z7barrierj(i32) #0
        define void @entry(i32 %0) #1 {
            call void @_Z7barrierj(i32 %0)
            ret void
        }
        attributes #0 = { convergent }
        attributes #1 = { alwaysinline }
        """, "entry"),
    Cvoid, Tuple{Int32}, convert(Int32, flags))
push!(opencl_builtins, "_Z7barrierj")
# TODO: add support for attributes to @builting_ccall/LLVM.@typed_ccall


text:
using LCMCore
using CaesarLCMTypes
using Rotations, CoordinateTransformations
using MeshCat
using TransformUtils
using AprilTags

tagData = Dict{Int,apriltag_t}()

function tag_handler(c,m,td)
    td[m.id] = deepcopy(m)
    nothing
end

lcm=LCM()

subscribe(lcm, "TEST_CHANNEL", (m,c)->tag_handler(m,c,tagData), apriltag_t)

@async begin
    while true
        tnow = Int64(time()*1000_000)
        @show keys(tagData)
        for tag in tagData
            if abs(tnow - tag[2].utime) <= 100_000
                msg = tag_landmark_t()
                msg.id = tag[2].id
                msg.bearing = atan2(tag[2].pose[3,2], tag[2].pose[3,3])
                msg.range = norm(tag[2].pose[1:3,4])
                msg.utime = tag[2].utime
                publish(lcm,"Landmarks",encode(msg))
            elseif 100_000 < abs(tnow - tag[2].utime) < 1000_000
                @warn "AprilTag pipeline is taking more than 100ms: $(round(abs(tnow - tag[2].utime)/1000.0, 2))"
            end
            delete!(tagData, tag[1])
        end
        sleep(2)
    end
end

while true
    handle(lcm)
    sleep(0.01)
end

#


text:
# !!! this script needs Julia v 0.4 or higher !!!
# constructing:
start_time = time()
# long computation
time_elapsed = time() - start_time
println("Time elapsed: $time_elapsed") #> 0.0010001659393310547

d = Date(2014,9,1) #> 2014-09-01
dt = DateTime(2014,9,1,12,30,59,1) #> 2014-09-01T12:30:59.001
# accessors:
Dates.year(d)
Dates.month(d)
Dates.week(d)
Dates.day(d)
# functions:
Dates.isleapyear(d)
Dates.dayofyear(d)
Dates.monthname(d)
Dates.daysinmonth(d)

max_stars_repo_path: Module 1/Chapter02/dates.jl

max_issues_repo_path: Module 1/Chapter02/dates.jl

max_forks_repo_path: Module 1/Chapter02/dates.jl


text:
module Deprecated

using Test
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "@based_on" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    m = [100, 200, 300, 400, 500]

    gq = :g
    iq = :i
    tq = :t
    yq = :y
    cq = :c

    gr = "g"
    ir = "i"
    tr = "t"
    yr = "y"
    cr = "c"

    gd = groupby(df, :g)

    n_str = "new_column"
    n_sym = :new_column
    n_space = "new column"

    @test @based_on(gd, n = mean(:i)).n == [2.0, 4.5]
    @test @based_on(gd, n = mean(:i) + mean(:g)).n == [3.0, 6.5]
    @test @based_on(gd, n = first(:t .* string.(:y))).n == ["av", "cy"]
    @test @based_on(gd, n = first(Symbol.(:y, ^(:t)))).n == [:vt, :yt]
    @test @based_on(gd, n = first(Symbol.(:y, ^(:body)))).n == [:vbody, :ybody]
    @test @based_on(gd, body = :i).body == df.i
    @test @based_on(gd, transform = :i).transform == df.i
    @test @based_on(gd, (n1 = [first(:i)], n2 = [first(:y)])).n1 == [1, 4]

    @test @based_on(gd, n = mean(cols(iq))).n == [2.0, 4.5]
    @test @based_on(gd, n = mean(cols(iq)) + mean(cols(gq))).n == [3.0, 6.5]
    @test @based_on(gd, n = first(cols(tq) .* string.(cols(yq)))).n == ["av", "cy"]
    @test @based_on(gd, n = first(Symbol.(cols(yq), ^(:t)))).n == [:vt, :yt]
    @test @based_on(gd, n = first(Symbol.(cols(yq), ^(:body)))).n == [:vbody, :ybody]
    @test @based_on(gd, body = cols(iq)).body == df.i
    @test @based_on(gd, transform = cols(iq)).transform == df.i
    @test @based_on(gd, (n1 = [first(cols(iq))], n2 = [first(cols(yq))])).n1 == [1, 4]

    @test @based_on(gd, n = mean(cols(ir))).n == [2.0, 4.5]
    @test @based_on(gd, n = mean(cols(ir)) + mean(cols(gr))).n == [3.0, 6.5]
    @test @based_on(gd, n = first(cols(tr) .* string.(cols(yr)))).n == ["av", "cy"]
    @test @based_on(gd, n = first(Symbol.(cols(yr), ^(:t)))).n == [:vt, :yt]
    @test @based_on(gd, n = first(Symbol.(cols(yr), ^(:body)))).n == [:vbody, :ybody]
    @test @based_on(gd, body = cols(ir)).body == df.i
    @test @based_on(gd, transform = cols(ir)).transform == df.i
    @test @based_on(gd, (n1 = [first(cols(ir))], n2 = [first(cols(yr))])).n1 == [1, 4]
    @test @based_on(gd, n = mean(cols("i")) + 0 * first(cols(:g))).n == [2.0, 4.5]
    @test @based_on(gd, n = mean(cols(2)) + first(cols(1))).n == [3.0, 6.5]

    @test @based_on(gd, :i) == select(df, :g, :i)
    @test @based_on(gd, :i, :g) ≅ select(df, :g, :i)

    @test @based_on(gd, :i, n = 1).n == fill(1, nrow(df))

    @test @based_on(gd, cols("new_column") = 2).new_column == [2, 2]
    @test @based_on(gd, cols(n_str) = 2).new_column == [2, 2]
    @test @based_on(gd, cols(n_sym) = 2).new_column == [2, 2]
    @test @based_on(gd, cols(n_space) = 2)."new column" == [2, 2]
    @test @based_on(gd, cols("new" * "_" * "column") = 2)."new_column" == [2, 2]
end

# `y` needs to be in global scope here because the testset relies on `y`
y = 0
@testset "byrow" begin
    df = DataFrame(A = 1:3, B = [2, 1, 2])

    @test @byrow(df, if :A > :B; :A = 0 end) == DataFrame(A = [1, 0, 0], B = [2, 1, 2])

    @test  df == DataFrame(A = [1, 2, 3], B = [2, 1, 2])

    df = DataFrame(A = 1:3, B = [2, 1, 2])  # Restore df
    @byrow(df, if :A + :B == 3; global y += 1 end)
    @test  y == 2

    df = DataFrame(A = 1:3, B = [2, 1, 2])
    df2 = @byrow df begin
        @newcol colX::Array{Float64}
        @newcol colY::Array{Float64}
        :colX = :B == 2 ? pi * :A : :B
        if :A > 1
            :colY = :A * :B
        end
    end

    @test  df2.colX == [pi, 1.0, 3pi]
    @test  df2[2, :colY] == 2

    df = DataFrame(b = [5], a = [(n1 = 1, n2 = 2)])
    df2 = @byrow df begin
        :b = :a.n1
    end

    @test df2.b == [1]

    _DF = 5
    _N = 0
    @byrow df begin
        # nothing
    end
    @test _DF == 5
    @test _N == 0

    df3 = @byrow df begin
        :b = row
    end
    @test df3.b == [1]

    x = (a = 400, b = 600)
    df4 = @byrow df begin
        :b = x.a
    end
    @test df4.b == [400]
end

@testset "cols with @byrow" begin
    df = DataFrame(A = 1:3, B = [2, 1, 2])
    n = :A
    df2 = @byrow df begin
        :B = cols(n)
    end
    @test df2 == DataFrame(A = 1:3, B = 1:3)

    n = "A"
    df2 = @byrow df begin
        :B = cols(n)
    end
    @test df2 == DataFrame(A = 1:3, B = 1:3)

    df2 = @byrow df begin
        :A = cols(:A) + cols("B")
    end
    @test df2.A == df.A + df.B

    n = :A
    df2 = @byrow df begin
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [2, 1, 2], B = [2, 1, 2])

    n = "A"
    df2 = @byrow df begin
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [2, 1, 2], B = [2, 1, 2])

    n = :C
    df2 = @byrow df begin
        @newcol cols(n)::Vector{Int}
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [1, 2, 3], B = [2, 1, 2], C = [2, 1, 2])

    n = "C"
    df2 = @byrow df begin
        @newcol cols(n)::Vector{Int}
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [1, 2, 3], B = [2, 1, 2], C = [2, 1, 2])
end

df = DataFrame(A = 1:3, B = [2, 1, 2])

# `y` needs to be in global scope here because the testset relies on `y`
y = 0
@testset "byrow!" begin
    df = DataFrame(A = 1:3, B = [2, 1, 2])

    @test @byrow!(df, if :A > :B; :A = 0 end) == DataFrame(A = [1, 0, 0], B = [2, 1, 2])

    @test  df == DataFrame(A = [1, 2, 3], B = [2, 1, 2])

    df = DataFrame(A = 1:3, B = [2, 1, 2])  # Restore df
    @byrow!(df, if :A + :B == 3; global y += 1 end)
    @test  y == 2

    df = DataFrame(A = 1:3, B = [2, 1, 2])
    df2 = @byrow! df begin
        @newcol colX::Array{Float64}
        @newcol colY::Array{Float64}
        :colX = :B == 2 ? pi * :A : :B
        if :A > 1
            :colY = :A * :B
        end
    end

    @test  df2.colX == [pi, 1.0, 3pi]
    @test  df2[2, :colY] == 2

    df = DataFrame(b = [5], a = [(n1 = 1, n2 = 2)])
    df2 = @byrow! df begin
        :b = :a.n1
    end

    @test df2.b == [1]

    _DF = 5
    _N = 0
    @byrow! df begin
        # nothing
    end
    @test _DF == 5
    @test _N == 0

    df3 = @byrow! df begin
        :b = row
    end
    @test df3.b == [1]

    x = (a = 400, b = 600)
    df4 = @byrow! df begin
        :b = x.a
    end
    @test df4.b == [400]
end

@testset "cols with @byrow!" begin
    df = DataFrame(A = 1:3, B = [2, 1, 2])
    n = :A
    df2 = @byrow! df begin
        :B = cols(n)
    end
    @test df2 == DataFrame(A = 1:3, B = 1:3)

    n = "A"
    df2 = @byrow! df begin
        :B = cols(n)
    end
    @test df2 == DataFrame(A = 1:3, B = 1:3)

    df2 = @byrow! df begin
        :A = cols(:A) + cols("B")
    end
    @test df2.A == df.A + df.B

    n = :A
    df2 = @byrow! df begin
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [2, 1, 2], B = [2, 1, 2])

    n = "A"
    df2 = @byrow! df begin
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [2, 1, 2], B = [2, 1, 2])

    n = :C
    df2 = @byrow! df begin
        @newcol cols(n)::Vector{Int}
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [1, 2, 3], B = [2, 1, 2], C = [2, 1, 2])

    n = "C"
    df2 = @byrow! df begin
        @newcol cols(n)::Vector{Int}
        cols(n) = :B
    end
    @test df2 == DataFrame(A = [1, 2, 3], B = [2, 1, 2], C = [2, 1, 2])
end

df = DataFrame(A = 1:3, B = [2, 1, 2])

end # module


text:
module NasaSpice

using StrPack

const libNasaSpice = Pkg.dir("NasaSpice", "deps", "libcspice")

# replacement for SpiceEllipse
#The members are:
#    center:     Vector defining ellipse's
#                center.
#    semiMajor:  Vector defining ellipse's
#                semi-major axis.
#    semiMinor:  Vector defining ellipse's
#                semi-minor axis.
#
# The ellipse is the set of points
#   {X:  X =                  center
#              + cos(theta) * semiMajor
#              + sin(theta) * semiMinor,
#    theta in [0, 2*Pi) }
@struct type Ellipse
    center::Array{Float64,1}(3)
    semiMajor::Array{Float64,1}(3)
    semiMinor::Array{Float64,1}(3)
end
Base.convert(::Type{Ellipse}) = Ellipse([0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0])

function codeForEllipse( phase::AbstractString, name::AbstractString, ctype::AbstractString )
    return code( phase, "Ellipse", name, ctype );
end

#e = NasaSpice.Ellipse( [0.,0.,0.], [0.,0.,0.], [0.,0.,0.] )
#io = IOBuffer()
#StrPack.pack( io, e )
#convert( Ptr{Float64}, convert( Ptr{Void}, pointer(io.data) ) )
#ccall...
#seek(io,0)
#f = StrPack.unpack( io, NasaSpice.Ellipse )

# replacement for SpicePlane
# Plane =
#               {X: <normal,X> = constant}
#  normal:     Vector normal to plane.
#  constant:   Constant of plane equation
@struct type Plane
    normal::Array{Float64,1}(3)
    constant::Float64;
end
Base.convert(::Type{Plane}) = Plane([0.0,0.0,0.0],0.0)
function codeForPlane( phase::AbstractString, name::AbstractString, ctype::AbstractString )
    return code( phase, "Plane", name, ctype );
end

function codeForAbstractString( phase::AbstractString, name::AbstractString, maxLength::Int64 )
    str = Array(AbstractString,0)
    # Array{UInt8}(32)
    if phase == "pack" # convert input structure to stream and pointer
        error("codeForAbstractString does not implement 'pack'")
    elseif phase == "allocate" # allocate output structure, convert to stream and pointer
        str = vcat( str, name * "_array = Array{UInt8}(" * string(maxLength) * ")" )
        str = vcat( str, name * "_ptr = pointer( " * name * "_array )" )
# warning on 63 is really in StrPack.jl @34 STRUCT_REGISTRY
    elseif phase == "unpack"   # convert output structure from stream
        name = replace( name, "_ptr", "")
        str = vcat( str, name * " = bytestring( " * name * "_ptr )" )
    end
    return str
end

function codeForAbstractString256( phase::AbstractString, name::AbstractString, ctype::AbstractString )
    codeForAbstractString( phase, name, 256 )
end

function code( phase::AbstractString, structure::AbstractString, name::AbstractString, ctype::AbstractString )
    str = Array(AbstractString,0)
    if phase == "pack" # convert input structure to stream and pointer
        str = vcat( str, "io" * name * " = IOBuffer()" )
        str = vcat( str, "StrPack.pack( io" * name * ", " * name * ")" )
        str = vcat( str, name * "_ptr = convert( Ptr{" * ctype * "}, convert( Ptr{Void}, pointer(io" * name * ".data) ) )" )
    elseif phase == "allocate" # allocate output structure, convert to stream and pointer
        str = vcat( str, name * " = " * structure * "()" )
        str = vcat( str, "io" * name * " = IOBuffer()" )
        str = vcat( str, "StrPack.pack( io" * name * ", " * name * ")" )
        str = vcat( str, name * "_ptr = convert( Ptr{" * ctype * "}, convert( Ptr{Void}, pointer(io" * name * ".data) ) )" )
    elseif phase == "unpack"   # convert output structure from stream
        name = replace( name, "_ptr", "")
        str = vcat( str, "seek(io" * name * ",0)" )
        str = vcat( str, name * " = StrPack.unpack( io" * name * ", " * structure * " )" )
    end
    return str
end

include("NasaSpiceApi.jl")
include("NasaSpiceApi_vector.jl")
include("NasaSpiceApi_stringOut.jl")

# Hand Conversions

function resetError()
    ccall((:reset_c,libNasaSpice),Void,())
end

function setErrorAction( action::AbstractString )
    ccall((:erract_c,libNasaSpice),Void,(Ptr{UInt8},Int64,Ptr{UInt8}),
    "SET",length(action),action)
end

include("cspice.jl")

## bitstype 64 dflFloat <: FloatingPoint
## Base.convert(::Type{NasaSpice.dflFloat}, x::Integer) = (real(x))
## Base.convert(::Type{NasaSpice.dflFloat}, x::FloatingPoint) = (x)

#Array_6_SpiceDouble(v::SpiceDouble) = Array_6_SpiceDouble(v,v,v,v,v,v)
#Array_6_Array_6_SpiceDouble(v::SpiceDouble) = Array_6_Array_6_SpiceDouble(
#    Array_6_SpiceDouble(v),
#    Array_6_SpiceDouble(v),
#    Array_6_SpiceDouble(v),
#    Array_6_SpiceDouble(v),
#    Array_6_SpiceDouble(v),
#    Array_6_SpiceDouble(v))

## include("SpiceZdf.jl")

## include("SpiceUsr.jl")

##include("/home/don/xf2eul.jl")

function furnsh_c(file::AbstractString)
    ccall((:furnsh_c,libNasaSpice),Void,(Ptr{UInt8},),file)
end

#function # axisar
##       (r::Array{Float64}(3,3)) =   # Array_3_Array_3_SpiceDouble
#    axisar(
#        axis::Array{Float64}, # Array_3_ConstSpiceDouble
#        angle::Float64) # SpiceDouble

##       enforce input array sizes
#    if length(axis) != 3
#        error("Incorrect size for parameter 1")
#    end
##       allocate the output parameters
#    r = Array{Float64}(3,3);  # Array_3_Array_3_SpiceDouble
#    r_ptr = pointer(r)

##       make transposed copies of all input arrays and their pointers
#    axis_t = axis'
#    axis_ptr = pointer(axis_t)

#    ccall((:axisar_c,libNasaSpice),Void,
#        (Ptr{Float64},Float64,Ptr{Float64}),
#        axis_ptr,angle,r_ptr)

#    return (r')
#end

#function str2et( date::AbstractString )
#    et = Array(SpiceDouble,1)
#    ccall((:str2et_c,libNasaSpice),Void,(Ptr{ConstSpiceChar},Ptr{SpiceDouble}),
#                date,
#                et )
#    return convert(Float64, et[1] )
#end

#function sxform( from::AbstractString, to::AbstractString, et::Float64 )
#    xform = zeros(SpiceDouble, 6, 6);
#    xformPtr = convert(Ptr{SpiceDouble}, pointer(xform))
#    et_c = convert(SpiceDouble, et );

#    ccall((:sxform_c,libNasaSpice),Void,(Ptr{ConstSpiceChar},Ptr{ConstSpiceChar},SpiceDouble,Ptr{SpiceDouble}),
#        from,
#        to,
#        et_c,
#        xformPtr )
#    return xform' # TODO remove this an transposes in C wrappers
#end

function __init__()
    setErrorAction( "REPORT" )
end

end

max_stars_repo_name: lintondf/NasaSpice

max_issues_repo_name: lintondf/NasaSpice

max_forks_repo_name: lintondf/NasaSpice


text:
# This file contains a list of the currently available covariance functions

abstract type Kernel end

"""
    KernelData

Data to be used with a kernel object to calculate a covariance matrix,
which is independent of kernel hyperparameters.

See also [`EmptyData`](@ref).
"""
abstract type KernelData end

"""
    EmptyData <: KernelData

Default empty `KernelData`.
"""
struct EmptyData <: KernelData end

KernelData(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix) = EmptyData()
kernel_data_key(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix) = "EmptyData"

"""
    cov(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix)

Create covariance matrix from kernel `k` and matrices of observations `X1` and `X2`, where
each column is an observation.
"""
function cov(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData=EmptyData())
    dim1, nobs1 = size(X1)
    dim2, nobs2 = size(X2)
    dim1==dim2 || throw(ArgumentError("X1 and X2 must have same dimension"))
    cK = Array{promote_type(eltype(X1), eltype(X2))}(undef, nobs1, nobs2)
    cov!(cK, k, X1, X2, data)
end

function cov!(cK::AbstractMatrix, k::Kernel, X::AbstractMatrix, data::KernelData=EmptyData())
    dim, nobs = size(X)
    (nobs,nobs) == size(cK) || throw(ArgumentError("cK has size $(size(cK)) and X has size $(size(X))"))
    @inbounds for j in 1:nobs
        cK[j,j] = cov_ij(k, X, X, data, j, j, dim)
        for i in 1:j-1
            cK[i,j] = cov_ij(k, X, X, data, i, j, dim)
            cK[j,i] = cK[i,j]
        end
    end
    return cK
end
"""
    cov!(cK::AbstractMatrix, k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix)

Like [`cov(k, X1, X2)`](@ref), but stores the result in `cK` rather than a new matrix.
"""
function cov!(cK::AbstractMatrix, k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData=EmptyData())
    if X1 === X2
        return cov!(cK, k, X1, data)
    end
    dim1, nobs1 = size(X1)
    dim2, nobs2 = size(X2)
    dim1==dim2 || throw(ArgumentError("X1 and X2 must have same dimension"))
    dim = size(X1, 1)
    (nobs1,nobs2) == size(cK) || throw(ArgumentError("cK has size $(size(cK)) X1 $(size(X1)) and X2 $(size(X2))"))
    @inbounds for i in 1:nobs1
        for j in 1:nobs2
            cK[i,j] = cov_ij(k, X1, X2, data, i, j, dim)
        end
    end
    return cK
end

"""
    cov(k::Kernel, X::AbstractMatrix[, data::KernelData = EmptyData()])

Create covariance matrix from kernel `k`, matrix of observations `X`, where each column is
an observation, and kernel data `data` constructed from input observations.
"""
cov(k::Kernel, X::AbstractMatrix, data::KernelData=EmptyData()) = cov(k, X, X, data)

@inline cov_ij(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, dim::Int) = cov(k, @view(X1[:,i]), @view(X2[:,j]))
# the default is to drop the KernelData
@inline cov_ij(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, dim::Int) = cov_ij(k, X1, X2, i, j, dim)

############################
##### Kernel Gradients #####
############################
@inline @inbounds function dKij_dθ!(dK::AbstractVector, kern::Kernel, X1::AbstractMatrix, X2::AbstractMatrix,
                                    data::KernelData, i::Int, j::Int, dim::Int, npars::Int)
    for iparam in 1:npars
        dK[iparam] = dKij_dθp(kern, X1, X2, data, i, j, iparam, dim)
    end
end

function grad_slice!(dK::AbstractMatrix, k::Kernel, X::AbstractMatrix, data::KernelData, p::Int)
    dim, nobs = size(X)
    (nobs,nobs) == size(dK) || throw(ArgumentError("dK has size $(size(dK)) and X has size $(size(X))"))
    @inbounds for j in 1:nobs
        dK[j,j] = dKij_dθp(k,X,X,data,j,j,p,dim)
        @simd for i in 1:(j-1)
            dK[i,j] = dKij_dθp(k,X,X,data,i,j,p,dim)
            dK[j,i] = dK[i,j]
        end
    end
    return dK
end
function grad_slice!(dK::AbstractMatrix, k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, p::Int)
    if X1 === X2
        return grad_slice!(dK, k, X1, data, p)
    end
    dim1, nobs1 = size(X1)
    dim2, nobs2 = size(X2)
    dim1==dim2 || throw(ArgumentError("X1 and X2 must have same dimension"))
    (nobs1,nobs2) == size(dK) || throw(ArgumentError("dK has size $(size(dK)) X1 $(size(X1)) and X2 $(size(X2))"))
    dim=dim1
    @inbounds for i in 1:nobs1
        @simd for j in 1:nobs2
            dK[i,j] = dKij_dθp(k,X1,X2,data,i,j,p,dim)
        end
    end
    return dK
end

# Calculates the stack [dk / dθᵢ] of kernel matrix gradients
function grad_stack!(stack::AbstractArray, k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData)
    @inbounds for p in 1:num_params(k)
        grad_slice!(view(stack, :, :, p), k, X1, X2, data, p)
    end
    stack
end

grad_stack!(stack::AbstractArray, k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix) =
    grad_stack!(stack, k, X1, X2, KernelData(k, X1, X2))

grad_stack(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix) = grad_stack(k, X1, X2, KernelData(k, X1, X2))

function grad_stack(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData)
    nobs1 = size(X1, 2)
    nobs2 = size(X2, 2)
    stack = Array{eltype(X)}(undef, nobs1, nobs2, num_params(k))
    grad_stack!(stack, k, X1, X2, data)
end

@inline function dKij_dθp(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(k, X1, X2, i, j, p, dim)
end

include("stationary.jl")
include("distance.jl")
include("lin.jl")               # Linear covariance function
include("se.jl")                # Squared exponential covariance function
include("rq.jl")                # Rational quadratic covariance function
include("mat.jl")               # Matern covariance function
include("periodic.jl")          # Periodic covariance function
include("poly.jl")              # Polnomial covariance function
include("noise.jl")             # White noise covariance function
include("const.jl")             # Constant (bias) covariance function

# Wrapped kernels
include("masked_kernel.jl")     # Masked kernels (apply to subset of X dims)
include("fixed_kernel.jl")      # Fixed kernels (fix some hyperparameters)

# Composite kernels
include("composite_kernel.jl")
include("pair_kernel.jl")
include("sum_kernel.jl")
include("prod_kernel.jl")

# Autodifferentiation
include("autodiff.jl")


text:
abstract type Sampler{T<:AbstractFloat} end

struct SingleSampler{T} <: Sampler{T} end
SingleSampler(::Type{T}) where T = SingleSampler{T}()

(::SingleSampler{T})(i::I, j::I) where {T,I} = [(convert(T,i),convert(T,j))]

weight(::SingleSampler) = 1

struct JitteredSampler{T,I<:Integer} <: Sampler{T}
    samples::I
    n::I
    d::T
    w::T
end

function JitteredSampler(::Type{T}, samples::I) where {T<:AbstractFloat, I<:Integer}
    n = isqrt(samples)
    JitteredSampler(n^2, n, one(T)/n, one(T)/samples)
end

function (s::JitteredSampler{T,I})(i::I, j::I) where {T,I}
    hcat([[(i - one(T)/2 + (ii - 1 + rand(T))*s.d,
            j - one(T)/2 + (jj - 1 + rand(T))*s.d)
           for ii in 1:s.n]
          for jj in 1:s.n]...)
end

weight(s::JitteredSampler) = s.w

export Sampler, SingleSampler, JitteredSampler, weight


text:
mutable struct VkPhysicalDeviceFeatures
    robustBufferAccess::vk.VkBool32
    fullDrawIndexUint32::vk.VkBool32
    imageCubeArray::vk.VkBool32
    independentBlend::vk.VkBool32
    geometryShader::vk.VkBool32
    tessellationShader::vk.VkBool32
    sampleRateShading::vk.VkBool32
    dualSrcBlend::vk.VkBool32
    logicOp::vk.VkBool32
    multiDrawIndirect::vk.VkBool32
    drawIndirectFirstInstance::vk.VkBool32
    depthClamp::vk.VkBool32
    depthBiasClamp::vk.VkBool32
    fillModeNonSolid::vk.VkBool32
    depthBounds::vk.VkBool32
    wideLines::vk.VkBool32
    largePoints::vk.VkBool32
    alphaToOne::vk.VkBool32
    multiViewport::vk.VkBool32
    samplerAnisotropy::vk.VkBool32
    textureCompressionETC2::vk.VkBool32
    textureCompressionASTC_LDR::vk.VkBool32
    textureCompressionBC::vk.VkBool32
    occlusionQueryPrecise::vk.VkBool32
    pipelineStatisticsQuery::vk.VkBool32
    vertexPipelineStoresAndAtomics::vk.VkBool32
    fragmentStoresAndAtomics::vk.VkBool32
    shaderTessellationAndGeometryPointSize::vk.VkBool32
    shaderImageGatherExtended::vk.VkBool32
    shaderStorageImageExtendedFormats::vk.VkBool32
    shaderStorageImageMultisample::vk.VkBool32
    shaderStorageImageReadWithoutFormat::vk.VkBool32
    shaderStorageImageWriteWithoutFormat::vk.VkBool32
    shaderUniformBufferArrayDynamicIndexing::vk.VkBool32
    shaderSampledImageArrayDynamicIndexing::vk.VkBool32
    shaderStorageBufferArrayDynamicIndexing::vk.VkBool32
    shaderStorageImageArrayDynamicIndexing::vk.VkBool32
    shaderClipDistance::vk.VkBool32
    shaderCullDistance::vk.VkBool32
    shaderFloat64::vk.VkBool32
    shaderInt64::vk.VkBool32
    shaderInt16::vk.VkBool32
    shaderResourceResidency::vk.VkBool32
    shaderResourceMinLod::vk.VkBool32
    sparseBinding::vk.VkBool32
    sparseResidencyBuffer::vk.VkBool32
    sparseResidencyImage2D::vk.VkBool32
    sparseResidencyImage3D::vk.VkBool32
    sparseResidency2Samples::vk.VkBool32
    sparseResidency4Samples::vk.VkBool32
    sparseResidency8Samples::vk.VkBool32
    sparseResidency16Samples::vk.VkBool32
    sparseResidencyAliased::vk.VkBool32
    variableMultisampleRate::vk.VkBool32
    inheritedQueries::vk.VkBool32

    VkPhysicalDeviceFeatures() =  new(false, false, false, false, false, false, false, false, false, false,
                                    false, false, false, false, false, false, false, false, false, false,
                                    false, false, false, false, false, false, false, false, false, false,
                                    false, false, false, false, false, false, false, false, false, false,
                                    false, false, false, false, false, false, false, false, false, false,
                                    false, false, false, false, false
                                    )
end

function build(this::VkPhysicalDeviceFeatures)::vk.VkPhysicalDeviceFeatures
    return vk.VkPhysicalDeviceFeatures(
        this.robustBufferAccess,
        this.fullDrawIndexUint32,
        this.imageCubeArray,
        this.independentBlend,
        this.geometryShader,
        this.tessellationShader,
        this.sampleRateShading,
        this.dualSrcBlend,
        this.logicOp,
        this.multiDrawIndirect,
        this.drawIndirectFirstInstance,
        this.depthClamp,
        this.depthBiasClamp,
        this.fillModeNonSolid,
        this.depthBounds,
        this.wideLines,
        this.largePoints,
        this.alphaToOne,
        this.multiViewport,
        this.samplerAnisotropy,
        this.textureCompressionETC2,
        this.textureCompressionASTC_LDR,
        this.textureCompressionBC,
        this.occlusionQueryPrecise,
        this.pipelineStatisticsQuery,
        this.vertexPipelineStoresAndAtomics,
        this.fragmentStoresAndAtomics,
        this.shaderTessellationAndGeometryPointSize,
        this.shaderImageGatherExtended,
        this.shaderStorageImageExtendedFormats,
        this.shaderStorageImageMultisample,
        this.shaderStorageImageReadWithoutFormat,
        this.shaderStorageImageWriteWithoutFormat,
        this.shaderUniformBufferArrayDynamicIndexing,
        this.shaderSampledImageArrayDynamicIndexing,
        this.shaderStorageBufferArrayDynamicIndexing,
        this.shaderStorageImageArrayDynamicIndexing,
        this.shaderClipDistance,
        this.shaderCullDistance,
        this.shaderFloat64,
        this.shaderInt64,
        this.shaderInt16,
        this.shaderResourceResidency,
        this.shaderResourceMinLod,
        this.sparseBinding,
        this.sparseResidencyBuffer,
        this.sparseResidencyImage2D,
        this.sparseResidencyImage3D,
        this.sparseResidency2Samples,
        this.sparseResidency4Samples,
        this.sparseResidency8Samples,
        this.sparseResidency16Samples,
        this.sparseResidencyAliased,
        this.variableMultisampleRate,
        this.inheritedQueries
    )
end
    # robustBufferAccess::vk.VkBool32
    # fullDrawIndexUint32::vk.VkBool32
    # imageCubeArray::vk.VkBool32
    # independentBlend::vk.VkBool32

function setGeometryShader(this::VkPhysicalDeviceFeatures, value::vk.VkBool32)::VkPhysicalDeviceFeatures
    this.geometryShader = value
    return this
end

function setTessellationShader(this::VkPhysicalDeviceFeatures, value::vk.VkBool32)::VkPhysicalDeviceFeatures
    this.tessellationShader = value
    return this
end

    # sampleRateShading::vk.VkBool32
    # dualSrcBlend::vk.VkBool32
    # logicOp::vk.VkBool32
    # multiDrawIndirect::vk.VkBool32
    # drawIndirectFirstInstance::vk.VkBool32
    # depthClamp::vk.VkBool32
    # depthBiasClamp::vk.VkBool32
    # fillModeNonSolid::vk.VkBool32
    # depthBounds::vk.VkBool32
    # wideLines::vk.VkBool32
    # largePoints::vk.VkBool32
    # alphaToOne::vk.VkBool32
    # multiViewport::vk.VkBool32
function setSamplerAnisotropy(this::VkPhysicalDeviceFeatures, value::vk.VkBool32)::VkPhysicalDeviceFeatures
    this.samplerAnisotropy = value
    return this
end
    # textureCompressionETC2::vk.VkBool32
    # textureCompressionASTC_LDR::vk.VkBool32
    # textureCompressionBC::vk.VkBool32
    # occlusionQueryPrecise::vk.VkBool32
    # pipelineStatisticsQuery::vk.VkBool32

function setVertexPipelineStoresAndAtomics(this::VkPhysicalDeviceFeatures, value::vk.VkBool32)::VkPhysicalDeviceFeatures
    this.vertexPipelineStoresAndAtomics = value
    return this
end

function setFragmentStoresAndAtomics(this::VkPhysicalDeviceFeatures, value::vk.VkBool32)::VkPhysicalDeviceFeatures
    this.fragmentStoresAndAtomics = value
    return this
end
    # shaderTessellationAndGeometryPointSize::vk.VkBool32
    # shaderImageGatherExtended::vk.VkBool32
    # shaderStorageImageExtendedFormats::vk.VkBool32
    # shaderStorageImageMultisample::vk.VkBool32
    # shaderStorageImageReadWithoutFormat::vk.VkBool32
    # shaderStorageImageWriteWithoutFormat::vk.VkBool32
    # shaderUniformBufferArrayDynamicIndexing::vk.VkBool32
    # shaderSampledImageArrayDynamicIndexing::vk.VkBool32
    # shaderStorageBufferArrayDynamicIndexing::vk.VkBool32
    # shaderStorageImageArrayDynamicIndexing::vk.VkBool32
    # shaderClipDistance::vk.VkBool32
    # shaderCullDistance::vk.VkBool32
    # shaderFloat64::vk.VkBool32
    # shaderInt64::vk.VkBool32
    # shaderInt16::vk.VkBool32
    # shaderResourceResidency::vk.VkBool32
    # shaderResourceMinLod::vk.VkBool32
    # sparseBinding::vk.VkBool32
    # sparseResidencyBuffer::vk.VkBool32
    # sparseResidencyImage2D::vk.VkBool32
    # sparseResidencyImage3D::vk.VkBool32
    # sparseResidency2Samples::vk.VkBool32
    # sparseResidency4Samples::vk.VkBool32
    # sparseResidency8Samples::vk.VkBool32
    # sparseResidency16Samples::vk.VkBool32
    # sparseResidencyAliased::vk.VkBool32
    # variableMultisampleRate::vk.VkBool32
    # inheritedQueries::vk.VkBool32


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -20.345339284625858 <= q <= 16.479136897808104)

                     add_NL_constraint(m, :(swish(0.24719801624822457 + 0.29438950671798914*swish(0.7291821855933653 + 0.701331929085816*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.7293079760660532*swish(0.05096633991424904 + 0.31501139336769945*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8601808654748551*swish(0.18255055897839467 + -0.7759686413687503*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.9720568561360068*swish(0.5199334546165302 + 0.4081309459110076*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + swish(-0.49402133773229506 + 0.2465550492028581*swish(0.7291821855933653 + 0.701331929085816*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.17626257051164496*swish(0.05096633991424904 + 0.31501139336769945*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8688281780619751*swish(0.18255055897839467 + -0.7759686413687503*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.3904178182819744*swish(0.5199334546165302 + 0.4081309459110076*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + swish(-0.6759660684871824 + 0.9311431633625067*swish(0.7291821855933653 + 0.701331929085816*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.2015839225926066*swish(0.05096633991424904 + 0.31501139336769945*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.01676614484829253*swish(0.18255055897839467 + -0.7759686413687503*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.07861065087286123*swish(0.5199334546165302 + 0.4081309459110076*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + swish(-0.08101830092609674 + 0.5708385401734484*swish(0.7291821855933653 + 0.701331929085816*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.9629021826429458*swish(0.05096633991424904 + 0.31501139336769945*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.4780289946412233*swish(0.18255055897839467 + -0.7759686413687503*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.25224290978348307*swish(0.5199334546165302 + 0.4081309459110076*swish(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*swish(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*swish(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*swish(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# You should not include or use this file directly
# it contains basic functionality of MPS.
export MPS,MPO,getindex,setindex!,overlap,multiply!,multiply,compress!,*,transpose,transpose!,dag,dag!,disp

abstract type TensorNetwork end
abstract type TensorArray <: TensorNetwork end
abstract type AbstractMPS <: TensorArray end
abstract type abstractMPO <: TensorArray end
getindex(tensors::TensorArray,target::Int) = getindex(tensors.tensor,target)
setindex!(tensors::TensorArray,tensor::AbstractArray,target::Int) = setindex!(tensors.tensor,tensor,target)

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
#########
#File: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\new_application_\wiki_load_change_stats_and_estimate.jl
#Project: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\new_application_
#Created Date: Friday April 23rd 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Thursday June 3rd 2021 5:18:20 pm
#Modified By:  Domenico Di Gangi
#-----
#Description:
#-----
########

#%%
using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"
DrWatson.greet()
using Distributed

begin
using TableView
using Blink
viewtab(df) = body!(Window(), showtable(df))
using Logging
using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets,DynNets, Utilities, ErgmRcall
using LinearAlgebra
using CSV
using DataFrames
using Dates
using PyPlot
pygui(true)

end

#%%
#load data
begin
data_path = "./data/wiki_tk/raw_data/"

df = CSV.read(data_path * "wiki_talk_daily_edges.csv", DataFrame)

df.date = Dates.Date.(df.date)

dates = unique(df.date)

# get sequence of edge lists at a given frequency
edge_list_T = [Matrix(gdf[[:source, :target]]) for gdf in groupby(df,[:date])]
edge_list_pres_T = [Utilities.edge_list_pres(l) for l in edge_list_T]

n_links_T = [size(e)[1] for e in edge_list_T]
n_nodes_T = [ length(unique(reshape(e, :))) for e in edge_list_T]

end

# t_ref = argmax(n_nodes_T)
# n_links_sub = [100, 500, 1000, 2000, 3000, 5000, 10000]
# time = [@time DynNets.stats_from_mat(model1, edge_list_pres(edge_list_pres_T[t_ref][1:n_links, :])) for n_links in n_links_sub]
# plot(n_links_sub, time)
# time = [@time DynNets.stats_from_mat(model_rec_p_star, edge_list_pres(edge_list_pres_T[t_ref][1:n_links, :])) for n_links in n_links_sub]

@elapsed dfEst = collect_results( datadir("wiki_tk", "ch_stats_present"))
viewtab(dfEst[Not(:ch_stats)])

row = dfEst[dfEst.ergmString .==  "edges + gwidegree(decay = 1.5, fixed = TRUE, cutoff=10) + gwodegree(decay = 1.5, fixed = TRUE, cutoff=10)",:][1,:]

row = dfEst[dfEst.ergmString .==  "edges + mutual",:][1,:]

row = dfEst[dfEst.ergmString .==  "edges + gwesp(decay = 0.25, fixed = TRUE, cutoff=10)",:][1,:]
begin
row = dfEst[dfEst.ergmString .==  "edges + mutual + gwidegree(decay = 1.5, fixed = TRUE, cutoff=10) + gwodegree(decay = 1.5, fixed = TRUE, cutoff=10)",:][1,:]

model = row.model
T = row.Tf
t0_train = 1
tend_train = 365
N = n_nodes_T[t0_train:tend_train]
obsT = row.ch_stats[t0_train:tend_train]

ENV["JULIA_DEBUG"] = ScoreDrivenERGM

res_est = DynNets.estimate_and_filter(model, N, obsT; show_trace = true)
fig, ax = DynNets.plot_filtered(model, N, res_est[3])
end

# res_conf = DynNets.conf_bands_given_SD_estimates(model, N, obsT, DynNets.unrestrict_all_par(model, res_est.vEstSdResParAll), res_est.ftot_0, [[0.975, 0.025]];  offset=0, plotFlag=true, parUncMethod = "WHITE-MLE", xval = dates[t0_train:tend_train], winsorProp = 0)
# res_conf.ax[1].set_title("wiki talk $(model.staticModel.ergmTermsString)")

est_SS = DynNets.estimate_single_snap_sequence(model, obsT)
# fig, ax = DynNets.plot_filtered(model, N, est_SS, ax=res_conf.ax, lineType = ".", lineColor="r")

dates[365]
using Statistics
var(res_est.fVecT_filt[:, 1:365], dims = 2)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
function explicit_crank_nicolson(mu, beta, p0, h, t0, t1, a0, a1)
    a = a0:h:a1
    t = t0:h:t1
    m = length(a)
    n = length(t)

    u = zeros(h, m, n)
    u[:, 1] = p0(a, t0)

    for i in 1:(n - 1)
        h2 = 0.5 * (t[i] + t[i + 1])
        u[:, i + 1] = (1 .- 0.5 * h * mu(a, h2)) ./ (1 .+ 0.5 * h * mu(a, h2)) * u[:, i]
        u[1, i + 1] = h * sum(beta(a, t[i + 1]) .* u[1:m, i + 1]) + 0.5 * h * u[1, i]
    end
    u
end


text:
# Copied from linalg/lu.jl
function ipiv2perm(v::PList{T}, maxi::Integer) where T
    p = T[1:maxi]
    @inbounds for i in 1:length(v)
        p[i], p[v[i]] = p[v[i]], p[i]
    end
    return p
end

ipiv2perm(v::PList) = ipiv2perm(v, length(v))

# A naive algorithm, maybe there exists a faster one.
# Build the pivot while permuting 1:n. Make
# inverse perm at the same time, so we know where to
# find p[i]
function perm2ipiv(p::PList)
    n = length(p)
    m = [1:n]
    v = similar(m)
    x = copy(m)
    @inbounds for i in 1:n
        j = x[p[i]]
        mi = m[i]
        mj = m[j]
        m[i],m[j] = mj,mi
        v[i] = j
        x[mj] = i
        x[mi] = j
    end
    return v
end


text:
# v0.12 deprecations

function ones(dims...)
  Base.depwarn("Flux.ones(size...) is deprecated, please use Flux.ones32(size...) or Base.ones(Float32, size...)", :ones)
  Base.ones(Float32, dims...)
end
ones(T::Type, dims...) = Base.ones(T, dims...)

function zeros(dims...)
  Base.depwarn("Flux.zeros(size...) is deprecated, please use Flux.zeros32(size...) or Base.zeros(Float32, size...)", :zeros)
  Base.zeros(Float32, dims...)
end
zeros(T::Type, dims...) = Base.zeros(T, dims...)

ones32(::Type, dims...) = throw(ArgumentError("Flux.ones32 is always Float32, use Base.ones to specify the element type"))
zeros32(::Type, dims...) = throw(ArgumentError("Flux.zeros32 is always Float32, use Base.zeros to specify the element type"))

# v0.13 deprecations


text:
module Command

using Random
using Base.Iterators

sequence(l::Int) = String(rand(['A','C','G','T'], l))
function line(s; col=80)
    n = length(s) ÷ col
    l =  ( s[(i-1)*col+1:(i*col)] for i in 1:n )

    length(s) % col == 0 && return join(l,"\n")
    return join(l,"\n") * "\n" * s[(n*col)+1:end]
end

function usage()
    println("usage: julia script/ancestors.jl -n [int] -l [int]")
    exit(2)
end

function main(args)
    N,L = 100, Int(1e5)
    i = 1
    while i ≤ length(args)
        if args[i] == "-N"
            N = parse(Int,args[i+1])
            i += 2
            continue
        end
        if args[i] == "-L"
            L = parse(Int,args[i+1])
            i += 2
            continue
        end
        break
    end

    for i in 1:N
        s = sequence(L)
        name = "isolate_$(i)"
        write(stdout, ">$(name)", '\n')
        write(stdout, line(s), '\n')
    end
    flush(stdout)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
import Snapshots
import Test

git = Snapshots.Utils._get_git_binary_path()

previous_directory = pwd()

temp_directory_1 = joinpath(mktempdir(), "TEMPGITREPOLOCAL")
mkpath(temp_directory_1)

temp_directory_2 = joinpath(mktempdir(), "TEMPGITREPOREMOTE")
mkpath(temp_directory_2)

cd(temp_directory_2)
run(`$(git) init --bare`)

cd(temp_directory_1)
run(`$(git) init`)
Snapshots.Utils.git_add_all!()
Snapshots.Utils.git_commit!(
    ;
    message="test commit 1",
    allow_empty=true,
    committer_name="test name",
    committer_email="test email",
    )
run(`git branch branch1`)
run(`git branch branch2`)
run(`git branch branch3`)

run(`git checkout master`)
Test.@test( typeof(Snapshots.Utils.git_version()) <: VersionNumber )
Test.@test( typeof(Snapshots.Utils.get_all_branches_local()) <: Vector{String} )
Test.@test( typeof(Snapshots.Utils.get_all_branches_local_and_remote()) <: Vector{String} )
Test.@test( typeof(Snapshots.Utils.get_current_branch()) <: String )
Test.@test( Snapshots.Utils.branch_exists("branch1") )
Test.@test( !Snapshots.Utils.branch_exists("non-existent-branch") )
Test.@test( !Snapshots.Utils.branch_exists("non-existent-but-create-me") )
Test.@test( typeof(Snapshots.Utils.checkout_branch!("branch1")) <: Nothing )
# Test.@test_throws(
#     ErrorException,
#     Snapshots.Utils.checkout_branch!("non-existent-branch"),
#     )
Test.@test_warn(
    "",
    Snapshots.Utils.checkout_branch!("non-existent-branch"; error_on_failure=false,),
    )
Test.@test(
    typeof(
        Snapshots.Utils.checkout_branch!("non-existent-but-create-me"; create=true)
        ) <: Nothing
    )
Snapshots.Utils.git_add_all!()
Snapshots.Utils.git_commit!(
    ;
    message="test commit 2",
    allow_empty=true,
    committer_name="test name",
    committer_email="test email",
    )
run(`git checkout master`)
Test.@test( Snapshots.Utils.branch_exists("branch1") )
Test.@test( !Snapshots.Utils.branch_exists("non-existent-branch") )
Test.@test( Snapshots.Utils.branch_exists("non-existent-but-create-me") )

run(`$(git) remote add origin $(temp_directory_2)`)
Test.@test( typeof(Snapshots.Utils.git_push_upstream_all!()) <: Nothing )

run(`git checkout master`)
include_patterns = Regex[
    r"^bRANCh1$"i,
    r"^bRanCh3$"i,
    ]
exclude_patterns = Regex[
    r"^brANcH3$"i,
    ]
branches_to_snapshot = Snapshots.Utils.make_list_of_branches_to_snapshot(
    ;
    default_branch = "maSTeR",
    include = include_patterns,
    exclude = exclude_patterns,
    )
Test.@test( length(branches_to_snapshot) == 2 )
Test.@test( length(unique(branches_to_snapshot)) == 2 )
Test.@test(
    length(branches_to_snapshot) == length(unique(branches_to_snapshot))
    )
Test.@test( branches_to_snapshot[1] == "branch1" )
Test.@test( branches_to_snapshot[2] == "master" )

cd(previous_directory)

Snapshots.Utils.delete_everything_except_dot_git!(temp_directory_1)
Snapshots.Utils.delete_only_dot_git!(temp_directory_2)

rm(temp_directory_1; recursive=true, force=true)
rm(temp_directory_2; recursive=true, force=true)


text:
function testerror_heatmap(test::Array{StatisticalMoments{Float64}},fm::FeatureMap,  ml::ML_model,
      options, output_type::AbstractString;
      strdir="default")

     pairs = unique([ td.label[1]/(td.label[2]^2), td.label[2]] for td in test);
    #pairs = unique([ td.volfrac/(td.particles[1].r^2), td.particles[1].r] for td in tds);
    zerohypo = mean(p.^2 for p in pairs);
    Ferrvol = mom -> Ferr_internal(mom, ml,fm,zerohypo,options, output_type)
    plotmom_heatmap(test,Ferrvol)
    if strdir == "default"
      strdir = isdir("../figures/") ? "../figures" : ""
    end
    savefig("$(strdir)heatmap-$(ml.name).pdf");
end

function Ferr_internal(mom::StatisticalMoments, ml::ML_model,fm::FeatureMap, zerohypo::Array{Float64},options, output_type::AbstractString)
    ypred = predict_linear(mom, fm, ml, options); ypred= ypred[1];
    rd = mom.label[2];
    if(output_type =="radius")
        ferror =  (ypred-rd)^2/zerohypo[2];
    else
        vl  =mom.label[1] ;
        vlf = vl/rd^2;
        ferror = (ypred-vlf)^2/zerohypo[1];
     end# +abs(ypred[2] -rd)/rd;   abs(ypred[1]-vlf)/vlf
    return ferror;

end

function plotmom_heatmap(moms::Array{StatisticalMoments{Float64}}, contract_tds = mom ->  mom.label[2] )
  label_f = mom -> [mom.label[1],mom.label[2]]
  labels = union([ label_f(td)  for td in moms])
  #tds_arr = [filter(td -> label_f(td) == l, moms) for l in labels]
  F_arr = [contract_tds(tds) for tds in moms];
  vols = sort([l[1] for l in labels])
  as = sort([l[2] for l in labels])
  z_matrix = NaN*zeros(length(as),length(vols))
  for i in indices(F_arr,1)
    n = find(labels[i][1] .== vols)
    m = find(labels[i][2] .== as)
    z_matrix[m,n] = F_arr[i]
  end
  heatmap(vols,as,z_matrix, xlabel="numfrac", ylabel="radius")
end


text:
using NumericalMethodsforEngineers

A = Float64[16 4 8;4 5 -4;8 -4 22]
b = Float64[4, 2, 5]
n = size(A, 1)
d = zeros(n)

ldlt!(A, d)
lower = zeros(3, 3)
for i in 1:n
  for j in 1:i
    lower[i, j] = A[i, j] / d[j]
  end
end

println("\nLower Triangular Factors: \n $lower")
println("\nDiagonal Terms: \n $d")

ldlfor!(A, b)
for i in 1:n
  A[i, :] = A[i, :] / d[i]
end

subbac!(A, b)
println("\nSolution Vector: \n $b")


text:
module Config

using AbInitioSoftwareBase.Commands: CommandConfig, MpiexecConfig
using Configurations: from_dict
using Express: myuparse
using Express.PhononWorkflow: Dfpt, RealSpaceForceConstants, PhononDispersion, VDos
using QuantumESPRESSO.Commands:
    QuantumESPRESSOConfig, PwxConfig, PhxConfig, Q2rxConfig, MatdynxConfig
using QuantumESPRESSO.Inputs.PWscf: PWInput
using QuantumESPRESSO.Inputs.PHonon: PhInput, Q2rInput, MatdynInput
using Express.PhononWorkflow.Config: RuntimeConfig, Template

import Configurations: convert_to_option
import Express.PhononWorkflow.Config: ExpandConfig

function (::ExpandConfig)(template::Template)
    inputs = map(
        (:scf, :dfpt, :q2r, :disp),
        (PWInput, PhInput, Q2rInput, MatdynInput),
    ) do field, T
        str = read(getproperty(template, field), String)
        parse(T, str)
    end
    return (; zip((:scf, :dfpt, :q2r, :disp), inputs)...)
end

convert_to_option(::Type{RuntimeConfig}, ::Type{CommandConfig}, dict) =
    QuantumESPRESSOConfig(;
        mpi = get(dict, "mpi", MpiexecConfig()),
        pw = get(dict, "pw", PwxConfig()),
        ph = get(dict, "ph", PhxConfig()),
        q2r = get(dict, "q2r", Q2rxConfig()),
        matdyn = get(dict, "matdyn", MatdynxConfig()),
    )

end


text:
module XGrad

export
    xdiff,
    xgrad,
    kgrad,
    @diffrule,
    VectorCodeGen,
    BufCodeGen,
    CuCodeGen,
    # re-exports
    #  from Espresso
    Espresso,
    squeeze_sum,
    @get_or_create,
    __construct,
    #  from LinearAlgebra
    mul!,
    # (other) helpers
    ungetindex,
    sum_grad

include("core.jl")

end # module


text:
# module StringsGJL

export readbytes, readlinebytes, filter_non_ascii, count_non_ascii, countsubstrings,
    compactstring

"""
    compactstring(d::AbstractDict)

Return a compact string representation of `d`.

Example:
```
compactstring(Dict(Void=>46,Float64=>21,Int64=>18))
"Void:46,Float64:21,Int64:18"
```
"""
compactstring(d::AbstractDict) = join([string(x[1] , ":", x[2]) for x in d], ",")

"""
    compactstring(::Array{String})

Return a compact string representation of `a`.
This just joins the strings with a comma.
"""
compactstring(a::Array{T}) where T <: AbstractString = join(a,",")

"""
    readbytes(filename::AbstractString, numbytes=0)

Read the contents of `filename`, but return a byte array, rather than a
string. If numbytes is greater than `0` return at most `numbytes` bytes.
Otherwise, read the entire file.
"""
function readbytes(filename::AbstractString, numbytes::Integer=0)
    numbytes = numbytes > 0 ? numbytes : filesize(filename)
    open(filename) do io
        read(io,numbytes) # something broken. read doc does not agree with methods
    end
end

"""
    readlinebytes(stream::IO=STDIN; chomp::Bool=false)

Warning: If `chomp` is true, an error may be raised. For instance
on reading the first (but not subsequent) lines from a file.
Read an array of bytes terminated by '\n' from `stream`.  Lines in the
input end with '\n' or "\r\n" or the end of an input stream. When
`chomp` is true, these trailing newline characters are removed from the
line before it is returned. When `chomp` is false, they are returned as
part of the line.
"""
function readlinebytes(s::IO=STDIN; chomp::Bool=false)
    line = readuntil(s, 0x0a)
    i = length(line)
    if !chomp || i == 0 || line[i] != 0x0a
        return line
    elseif i < 2 || line[i-1] != 0x0d
        return resize!(line,i-1)
    else
        return resize!(line,i-2)
    end
end

"""
    filter_non_ascii(s)

return array of non-ascii characters in string `s::String`, or non-ascii bytes in `s::Array`.
"""
filter_non_ascii(s::AbstractString) = convert(Array{Char,1}, filter(c -> ! isascii(c), s))
filter_non_ascii(a::AbstractArray) = filter(c -> ! isascii(Char(c)), a)

"""
    count_non_ascii(s)

count number of non-ascii characters in string `s::String`, or non-ascii bytes in `s::Array`.
"""
count_non_ascii(str::AbstractString) =  count(c -> ! isascii(c), str)
count_non_ascii(a::AbstractArray) =  count(c -> ! isascii(Char(c)), a)

"""
    countsubstrings(where::String, what::String)

Count the number of occurrences of `what` in `where`. This
is copied from discourse or stackexchange.
"""
function countsubstrings(where::String, what::String)
    numfinds = 0
    starting = 1
    while true
        location = search(where, what, starting)
        isempty(location) && return numfinds
        numfinds += 1
        starting = location.stop + 1
    end
end

"""
    countsubstrings(where::String, what::String)

Count the number of occurrences of `what` in `where`. This
is copied from discourse or stackexchange.
"""
function allstringmatch(where::String, what)
    numfinds = 0
    starting = 1
    res = String[]
    while true
        location = search(where, what, starting)
        isempty(location) && return break
        push!(res,where[location])
        numfinds += 1
        starting = location.stop + 1
    end
    return res
end

# end # module


text:
using Test

A = Int64[1, 2, 3, 4]
B = Complex{Int64}[5+6im, 7+8im, 9+10im]
# getindex
@test reinterpret(Complex{Int64}, A) == [1 + 2im, 3 + 4im]
@test reinterpret(Float64, A) == reinterpret.(Float64, A)

@test reinterpret(NTuple{3, Int64}, B) == [(5,6,7),(8,9,10)]

# setindex
let Ac = copy(A), Bc = copy(B)
    reinterpret(Complex{Int64}, Ac)[2] = -1 - 2im
    @test Ac == [1, 2, -1, -2]
    reinterpret(NTuple{3, Int64}, Bc)[2] = (4,5,6)
    @test Bc == Complex{Int64}[5+6im, 7+4im, 5+6im]
    reinterpret(NTuple{3, Int64}, Bc)[1] = (1,2,3)
    @test Bc == Complex{Int64}[1+2im, 3+4im, 5+6im]

    A1 = reinterpret(Float64, A)
    A2 = reinterpret(Complex{Float64}, A)
    A1[1] = 1.0
    @test real(A2[1]) == 1.0
end

# same-size reinterpret where one of the types is non-primitive
let a = NTuple{4,UInt8}[(0x01,0x02,0x03,0x04)]
    @test reinterpret(Float32, a)[1] == reinterpret(Float32, 0x04030201)
    reinterpret(Float32, a)[1] = 2.0
    @test reinterpret(Float32, a)[1] == 2.0
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
# Functionality to download, clean and parse the text
# files from NCBI into a Julia native datastructure

"""
    clean_genomeset(outio::IO, inio::IO)

Reads in data from `inio`, representing the decompressed data of "genomeset.dat.gz"
and write a cleaned copy to `outio`. The cleaning process renames mistyped or invalid data.

Because it is implemented by manually correcting each observed mistake,
this function may fail, or incompletely clean the data in future versions of
the NCBI Influenza data. In that case, you can assume that parsing incompletely
cleaned data will fail, and not parse invalid data.
"""
function clean_genomeset(outio::IO, inio::IO)
    fields = Vector{SubString{String}}(undef, 11)

    for line in eachline(inio) |> imap(strip) |> ifilter(!isempty)
        # Check the correct number of fields
        if unwrap_or(try_split!(fields, line, '\t'), 0) < 11
            error("Error: Should have 11 fields \"$line\"")
        end

        gi, host, segment, subtype, country, year, len, isolate, age, gender, group = fields

        # Filter isolate
        maybe_alphaname = parse_genomeset_isolate(isolate)
        (isalpha, isolate) = if is_error(maybe_alphaname)
            nothing, ""
        else
            unwrap(maybe_alphaname)
        end

        # Filter subtype - set it to "" if not an influenza A
        subtype_upper = uppercase(subtype)
        subtype = if isalpha === nothing
            ""
        elseif !isalpha
            ""
        elseif startswith(subtype_upper, "MIXED")
            ""
        elseif subtype in ["H1", "H11N9/N2"]
            ""
        elseif subtype in ["H1N2v", "H3N2v"]
            subtype[1:end-1]
        elseif subtype == "H3N6,H3"
            "H3N6"
        elseif subtype == "H6N1,H6"
            "H6N1"
        else
            subtype_upper
        end

        # Filter year
        year = if year in ["NON", "Unknown", "unknown"]
            ""
        else
            year
        end

        println(outio, join([gi, host, segment, subtype, country, year, len, isolate, age, gender, group], '\t'))
    end
end

# Parses a isolate in genomeset.dat during cleaning, returning none if it's malformed.
# If it's unexpectedly malformed, throw an error"
# Returns Option(is_alphavirus, name)
function parse_genomeset_isolate(s::Union{String, SubString})::Option{Tuple{Bool, String}}
    # They look like this: "Influenza A virus (A/swine/Minnesota/66960/2006(H3N2))"
    isempty(s) && return none

    # Strip "Influenza A Virus" beginning off
    occursin(r"^Influenza [AB] [vV]irus", s) || error(s)
    ncodeunits(s) == 17 && return none
    rest = SubString(s, 18 + (codeunit(s, 18) == UInt(' ')):ncodeunits(s))

    # Strip parenthesis off
    isempty(rest) && return none
    rest2 = if codeunit(rest, 1) == UInt8('(') && codeunit(rest, ncodeunits(rest)) == UInt8(')')
        SubString(rest, 2:ncodeunits(rest)-1)
    else
        rest
    end

    # Some of them ends with "(mixed)" - these are useless and not real isolates.
    if endswith(rest2, "(mixed)")
        return none
    end

    # Now it should look like e.g. "A/mallard/Interior Alaska/8BM3586/2008(H3N8)"
    # strip off the serotype designation, if present.
    m = match(r"\([Hh]\d+[nN]\d+\)$", rest2)
    rest3 = if m === nothing
        rest2
    else
        SubString(rest2, 1:ncodeunits(rest2) - ncodeunits(m.match))
    end

    # Now it should look like "A/northern shoveler/Mississippi/09OS168/2009"
    # verify it actually does so!
    validate_isolate_format(rest3) || return none
    isalpha = first(rest3) == 'A'
    return some((isalpha, String(rest3)))
end

function validate_isolate_format(s::AbstractString)
    m = match(r"^([AB])/[^/]+/[^/]+(?:/[^/]+)?/\d{4}$", s)
    return m !== nothing
end

###
function try_parse_from_integer(::Type{Segment}, s::AbstractString)::Option{Segment}
    y = tryparse(UInt8, s)
    y === nothing && return none
    (iszero(y) | (y > 0x08)) && return none
    some(reinterpret(Segment, y - 0x01))
end

function try_parse(::Type{SeroType}, s::AbstractString)::Option{SeroType}
    st = tryparse(SeroType, s)
    st === nothing ? none : some(st)
end

function try_parse(
    ::Type{IncompleteSegmentData},
    fields::Vector{SubString{String}},
    line::Union{String, SubString{String}}
)::Option{IncompleteSegmentData}
    nfields = try_split!(fields, line, '\t')
    unwrap_or(nfields, 0) == 11 || error("Expected 11 fields in \"$line\"")
    serotype = if isempty(strip(fields[4]))
        none(SeroType)
    else
        some(@? try_parse(SeroType, fields[4]))
    end
    year = @? try_parse_year(fields[6])
    host = parse(Host, fields[2])
    gi = String(fields[1])
    isolate = String(fields[8])
    segment = @? try_parse_from_integer(Segment, fields[3])
    some(IncompleteSegmentData(gi, host, segment, serotype, year, isolate, ReferenceProtein[], none(LongDNASeq)))
end

function try_parse(
    ::Type{IncompleteSegmentData},
    line::Union{String, SubString{String}}
)::Option{IncompleteSegmentData}
    try_parse(IncompleteSegmentData, Vector{SubString{String}}(undef, 11), line)
end

function try_parse_year(s::Union{String, SubString})::Option{Int16}
    isempty(s) && return none
    pos_slash_found = findfirst(isequal('/'), s)
    last_byte = pos_slash_found === nothing ? ncodeunits(s) : pos_slash_found - 1
    return some(parse(Int16, SubString(s, 1:last_byte)))
end

function parse_cleaned_genomeset(io::IO)::Dict{String, IncompleteSegmentData}
    result = Dict{String, IncompleteSegmentData}()
    buffer = Vector{SubString{String}}(undef, 11)
    for line in eachline(io) |> imap(strip) |> ifilter(!isempty)
        data = @unwrap_or try_parse(IncompleteSegmentData, buffer, line) continue
        gi = data.id
        if haskey(result, gi)
            error("GB identifier $(gi) not unique")
        end
        result[gi] = data
    end
    result
end

"""
    parse_cleaned_genomeset(::IO) -> Dict{String, IncompleteSegmentData}

Parse the file generated by `clean_genomeset`, and return a `Dict` where the
keys are the GI accession number, and the value is a `IncompleteSegmentData` object.

This function attempts to validate the input data, but SKIPS any row with critical
fields missing or invalid. Erroring on the first invalid record would make the
dataset cleaning too big a task. Therefore, when updating the dataset, you should check
how many individual records are returned, or if an alarming number are skipped.
"""
parse_cleaned_genomeset(s::AbstractString) = open(parse_cleaned_genomeset, s)


text:
export AbstractField, AbstractMetaData, BasicMetaData
export Field

abstract type AbstractField <: AbstractExpression end
abstract type AbstractMetaData <: AbstractExpression end

struct Field{𝒯, 𝒮} <: AbstractField
    data::𝒯
    metadata::𝒮
end

struct BasicMetaData{𝒩} <: AbstractMetaData
    name::𝒩
end

Field() = Field(nothing, nothing)
Field(a; metadata = nothing) = Field(a, metadata)
Field(a::Number) = Field(a, BasicMetaData(string(a)))

# Interpret Numbers as special Fields
for binary_operator in [binary_operators..., ["Negative", "-"]]
    b_name, b_symbol = Meta.parse.(binary_operator)
    @eval $b_symbol(a::Number, b::AbstractExpression) = $b_symbol(Field(a), b)
    @eval $b_symbol(a::AbstractExpression, b::Number) = $b_symbol(a, Field(b))
end

function Base.show(io::IO, ϕ::Field{S, T}) where {S, T <: AbstractMetaData}
    printstyled(io, ϕ.metadata.name, color = 199  )
 end

compute(n::Number) = n
compute(a::Field) = compute(a.data)

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
# ---
# title: 183. Customers Who Never Order
# id: problem183
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/customers-who-never-order/description/>
# hidden: true
# ---
#
# Suppose that a website contains two tables, the `Customers` table and the
# `Orders` table. Write a SQL query to find all customers who never order
# anything.
#
# Table: `Customers`.
#
#
#
#     +----+-------+
#     | Id | Name  |
#     +----+-------+
#     | 1  | Joe   |
#     | 2  | Henry |
#     | 3  | Sam   |
#     | 4  | Max   |
#     +----+-------+
#
#
# Table: `Orders`.
#
#
#
#     +----+------------+
#     | Id | CustomerId |
#     +----+------------+
#     | 1  | 3          |
#     | 2  | 1          |
#     +----+------------+
#
#
# Using the above tables as example, return the following:
#
#
#
#     +-----------+
#     | Customers |
#     +-----------+
#     | Henry     |
#     | Max       |
#     +-----------+
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end