text:
import LinearAlgebra: norm

struct TrackHydrodynamicErrors{T} <: ProcessingMethod
    problem::FluidFlowProblem
    should_process::Bool
    n_steps::Int64
    stop_criteria::StopCriteria
    df::T
end
TrackHydrodynamicErrors(
    problem,
    should_process,
    n_steps,
    stop_criteria = StopCriteria(problem),
) = TrackHydrodynamicErrors(
    problem,
    should_process,
    n_steps,
    stop_criteria,
    Vector{
        NamedTuple{
            (
                :timestep,
                :error_ρ,
                :error_u,
                :error_p,
                :error_σ_xx,
                :error_σ_xy,
                :error_σ_yy,
                :error_σ_yx,
                :mass,
                :momentum,
                :energy,
            ),
            Tuple{
                Int64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
            },
        },
    }(),
)

function next!(process_method::TrackHydrodynamicErrors, q, f_in, t::Int64)
    should_stop = false

    if mod(t, 100) == 0
        if (should_stop!(process_method.stop_criteria, q, f_in))
            @info "Stopping after $t iterations"
            should_stop = true
        end
    end

    if (should_stop == false && t != process_method.n_steps)
        if (!process_method.should_process)
            return false
        end
    end

    problem = process_method.problem
    nx, ny, nf = size(f_in)
    x_range, y_range = range(problem)

    f = Array{Float64}(undef, size(f_in, 3))
    u = zeros(dimension(q))
    expected_u = zeros(dimension(q))

    error_ρ = 0.0
    error_u = 0.0
    error_p = 0.0
    error_σ_xx = 0.0
    error_σ_xy = 0.0
    error_σ_yy = 0.0
    error_σ_yx = 0.0

    time = t * delta_t(problem)
    Δ = Float64(y_range.step) * Float64(x_range.step)
    if (nx == 1)
        Δ = Float64(y_range.step)
        if (ny == 1)
            Δ = 1.0
        end
    elseif ny == 1
        Δ = Float64(x_range.step)
    end
    Δ_ = Δ
    Δ = 1.0

    τ = q.speed_of_sound_squared * LatticeBoltzmann.lattice_viscosity(problem) + 0.5

    D = dimension(q)
    N = div(LatticeBoltzmann.order(q), 2)
    N = 2
    Hs = [[hermite(Val{n}, q.abscissae[:, i], q) for i in 1:length(q.weights)] for n in 1:N]

    total_expected_u_squared = 0.0
    total_expected_σ_xx_squared = 0.0
    total_expected_σ_yy_squared = 0.0
    total_expected_σ_xy_squared = 0.0
    total_expected_σ_yx_squared = 0.0
    total_expected_p_squared = 0.0

    total_mass = 0.0
    total_momentum = 0.0
    total_energy = 0.0
    @inbounds for x_idx in 1:nx, y_idx in 1:ny
        # Analytical
        x = x_range[x_idx]
        y = y_range[y_idx]

        # Compute expected marcroscopic variables
        expected_ρ = LatticeBoltzmann.density(q, problem, x, y, time)
        expected_u = LatticeBoltzmann.velocity(problem, x, y, time)
        expected_p = LatticeBoltzmann.pressure(q, problem, x, y, time)
        expected_ϵ = (dimension(q) / 2) * expected_p / expected_ρ
        expected_T = expected_p / expected_ρ
        expected_σ = deviatoric_tensor(q, problem, x, y, time)

        total_expected_u_squared += sum(expected_u .^ 2)# norm(expected_u, 2)^2
        total_expected_σ_xx_squared += expected_σ[1, 1]^2
        total_expected_σ_yy_squared += expected_σ[2, 2]^2
        total_expected_σ_xy_squared += expected_σ[1, 2]^2
        total_expected_σ_yx_squared += expected_σ[2, 1]^2
        total_expected_p_squared += expected_p^2

        # Compute macroscopic variables
        @inbounds for f_idx in 1:size(f_in, 3)
            f[f_idx] = f_in[x_idx, y_idx, f_idx]
        end

        ρ = density(q, f)
        velocity!(q, f, ρ, u)
        p = pressure(q, f, ρ, u)

        # Adding the forcing term moves the optimal tau for poiseuille flows
        # F = cm.force(x_idx, y_idx, 0.0)
        # u += cm.τ * F

        # Hermite coefficients of \bar{f}
        τ = q.speed_of_sound_squared * LatticeBoltzmann.lattice_viscosity(problem)
        a_bar_2 = sum([f[idx] * Hs[2][idx] for idx in 1:length(q.weights)])
        a_eq_2 = equilibrium_coefficient(Val{2}, q, ρ, u, 1.0)

        # Determin a^2 of f based on \bar{f} and f^eq
        a_2 = (a_bar_2 + (1 / (2 * τ)) * a_eq_2) / (1 + 1 / (2 * τ))

        # Second order convergence?
        # P = (1 - 1 / (2 * τ))a_f[2] - (a_f[1] * a_f[1]') / ρ + ρ * I
        p = 0.0
        cs = 1 / q.speed_of_sound_squared
        for x_idx in 1:D
            # NOTE the (1 - 1/2τ) term probably comes from
            # σ = - (1 - 1/2τ) ∑ f^(1)_i c_i c_i
            # p += (1 - 1 / (2 * τ)) * a_bar_2[x_idx, x_idx] - ρ * (u[x_idx] * u[x_idx] - I)

            # Should be identical / better ?
            p += a_bar_2[x_idx, x_idx] - ρ * (u[x_idx] * u[x_idx] - I)
            # p += a_bar_2[x_idx, x_idx] - ρ * u[x_idx] * u[x_idx] + D * ρ * (1 + (1 - cs) / (2 * τ))
        end
        p /= D
        p1 = p
        # p = tr(a_bar_2 - ρ * (u * u' - I) ) / D

        # Huidige poging
        P = a_2 - ρ * (u * u' - I)
        σ_lb = P - I * tr(P) / D

        τ = q.speed_of_sound_squared * LatticeBoltzmann.lattice_viscosity(problem)
        σ_lb = deviatoric_tensor(q, τ, f, ρ, u)
        # σ_lb = deviatoric_tensor(q, problem.τ * q.speed_of_sound_squared, f, ρ, u)

        # Gives 4th order convergence?
        p = tr(P) / D

        # Determine errors by first scaling from lattice variables to dimensionless
        ρ = dimensionless_density(problem, ρ)
        u = dimensionless_velocity(problem, u)
        p = dimensionless_pressure(q, problem, p)
        σ_lb = dimensionless_stress(problem, σ_lb)

        # - \mu \rho T \Lambda
        σ_err = (expected_σ .- σ_lb)

        error_p += Δ * (p - expected_p)^2
        error_ρ += Δ * (ρ - expected_ρ)^2
        error_u += Δ * ((u[1] - expected_u[1])^2 + (u[2] - expected_u[2])^2)
        error_σ_xx += Δ * (expected_σ[1, 1] .- σ_lb[1, 1])^2
        error_σ_xy += Δ * (expected_σ[1, 2] .- σ_lb[1, 2])^2
        error_σ_yx += Δ * (expected_σ[2, 1] .- σ_lb[2, 1])^2
        error_σ_yy += Δ * (expected_σ[2, 2] .- σ_lb[2, 2])^2

        total_mass += ρ
        total_momentum += ρ * sum(u)
        total_energy += ρ * sum(u .^ 2)
    end

    push!(
        process_method.df,
        (
            timestep = t,
            error_ρ = sqrt(error_ρ),
            error_u = sqrt(error_u / total_expected_u_squared),
            # error_u = sqrt(error_u),# / total_expected_u_squared),
            error_p = sqrt(error_p / total_expected_p_squared),
            error_σ_xx = sqrt(error_σ_xx / total_expected_σ_xx_squared),
            error_σ_xy = sqrt(error_σ_xy / total_expected_σ_xy_squared),
            error_σ_yy = sqrt(error_σ_yy / total_expected_σ_yy_squared),
            error_σ_yx = sqrt(error_σ_yx / total_expected_σ_yx_squared),
            mass = Δ_ * total_mass,
            momentum = Δ_ * total_momentum,
            energy = Δ_ * total_energy,
        ),
    )

    if mod(t, 1) == 0
        should_visualize = false
        if (process_method.should_process)
            if t == process_method.n_steps
                should_visualize = true
            end

            if mod(t, max(10, round(Int, process_method.n_steps / 25))) == 0
                should_visualize = true
            end
        end

        if (should_visualize)
            Δt = delta_t(process_method.problem)
            visualize(process_method.problem, q, f_in, time, process_method.df)
        end
    end

    return should_stop
end


text:
# input: nz, nx, dz, dx, nSteps, nPoints_pml, nPad, dt, f0, survey_fname, data_dir_name, scratch_dir_name, isAc
using JSON
using DataStructures
using Dierckx
using Statistics
using PyCall
using LinearAlgebra
np = pyimport("numpy")

function paraGen(nz, nx, dz, dx, nSteps, dt, f0, nPml, nPad, para_fname, survey_fname, data_dir_name; if_win=false, filter_para=nothing, if_src_update=false, scratch_dir_name::String="", if_cross_misfit=false)

  para = OrderedDict()
  para["nz"] = nz
	para["nx"] = nx
	para["dz"] = dz
	para["dx"] = dx
	para["nSteps"] = nSteps
	para["dt"] = dt
	para["f0"] = f0
	para["nPoints_pml"] = nPml
	para["nPad"] = nPad
  if if_win != false
    para["if_win"] = true
  end
  if filter_para != nothing
    para["filter"] = filter_para
  end
  if if_src_update != false
    para["if_src_update"] = true
  end
	para["survey_fname"] = survey_fname
  para["data_dir_name"] = data_dir_name
  if !isdir(data_dir_name)
    mkdir(data_dir_name)
  end
  # if nStepsWrap != nothing
  #   para["nStepsWrap"] = nStepsWrap
  # end
  if if_cross_misfit != false
    para["if_cross_misfit"] = true
  end

	if(scratch_dir_name != "")
      para["scratch_dir_name"] = scratch_dir_name
      if !isdir(scratch_dir_name)
        mkdir(scratch_dir_name)
      end
  end
  para_string = JSON.json(para)

  open(para_fname,"w") do f
    write(f, para_string)
  end
end

# all shots share the same number of receivers
function surveyGen(z_src, x_src, z_rec, x_rec, survey_fname; Windows=nothing,
  Weights=nothing, Src_Weights=nothing, Src_rxz=nothing, Rec_rxz=nothing)

  nsrc = length(x_src)
  nrec = length(x_rec)
  survey = OrderedDict()
  survey["nShots"] = nsrc
  for i = 1:nsrc
    shot = OrderedDict()
    shot["z_src"] = z_src[i]
    shot["x_src"] = x_src[i]
    shot["nrec"] = nrec
    shot["z_rec"] = z_rec
    shot["x_rec"] = x_rec
    if Windows != nothing
      shot["win_start"] = Windows["shot$(i-1)"][:start]
      shot["win_end"] = Windows["shot$(i-1)"][:end]
    end
    if Weights != nothing
      # shot["weights"] = Int64.(Weights["shot$(i-1)"][:weights])
      shot["weights"] = Weights["shot$(i-1)"][:weights]
    end
    if Src_Weights != nothing
      # shot["weights"] = Int64.(Weights["shot$(i-1)"][:weights])
      shot["src_weight"] = Src_Weights[i]
    end
    if Src_rxz != nothing
      shot["src_rxz"] = Src_rxz[i]
    end
    if Rec_rxz != nothing
      shot["rec_rxz"] = Rec_rxz
    end
    survey["shot$(i-1)"] = shot
  end

  survey_string = JSON.json(survey)
  open(survey_fname,"w") do f
    write(f, survey_string)
  end

end

function sourceGene(f, nStep, delta_t)
#  Ricker wavelet generation and integration for source
#  Dongzhuo Li @ Stanford
#  May, 2015

  e = pi*pi*f*f;
  t_delay = 1.2/f;
  source = Matrix{Float64}(undef, 1, nStep)
  for it = 1:nStep
      source[it] = (1-2*e*(delta_t*(it-1)-t_delay)^2)*exp(-e*(delta_t*(it-1)-t_delay)^2);
  end
  # return source

  for it = 2:nStep
      source[it] = source[it] + source[it-1];
  end
  source = source * delta_t;
end

# get vs high and low bounds from log point cloud
# 1st row of Bounds: vp ref line
# 2nd row of Bounds: vs high ref line
# 3rd row of Bounds: vs low ref line
function cs_bounds_cloud(cpImg, Bounds)
  cs_high_itp = Spline1D(Bounds[1,:], Bounds[2,:]; k=1)
  cs_low_itp = Spline1D(Bounds[1,:], Bounds[3,:]; k=1)
  csHigh = zeros(size(cpImg))
  csLow = zeros(size(cpImg))
  for i = 1:size(cpImg, 1)
    for j = 1:size(cpImg, 2)
      csHigh[i,j] = min(cs_high_itp(cpImg[i,j]), cpImg[i,j]/sqrt(2)-1.0)
      csLow[i,j] = cs_low_itp(cpImg[i,j])
    end
  end
  return csHigh, csLow
end

function klauderWave(fmin, fmax, t_sweep, nStepTotal, nStepDelay, delta_t)
#  Klauder wavelet
#  Dongzhuo Li @ Stanford
#  August, 2019
  nStep = nStepTotal - nStepDelay
  source = Matrix{Float64}(undef, 1, nStep+nStep-1)
  source_half = Matrix{Float64}(undef, 1, nStep-1)
  K = (fmax - fmin) / t_sweep
  f0 = (fmin + fmax) / 2.0
  t_axis = delta_t:delta_t:(nStep-1)*delta_t
  source_half = sin.(pi * K .* t_axis .* (t_sweep .- t_axis)) .* cos.(2.0 * pi * f0 .* t_axis) ./ (pi*K.*t_axis*t_sweep)
  for i = 1:nStep-1
    source[i] = source_half[end-i+1]
  end
  for i = nStep+1:2*nStep-1
    source[i] = source_half[i-nStep]
  end
  source[nStep] = 1.0
  source_crop = source[:,nStep-nStepDelay:end]
  return source_crop
end

# function klauderWave(fmin, fmax, t_sweep, nStep, delta_t)
# #  Klauder wavelet
# #  Dongzhuo Li @ Stanford
# #  August, 2019
#   source = Matrix{Float64}(undef, 1, nStep)
#   K = (fmax - fmin) / t_sweep
#   f0 = (fmin + fmax) / 2.0
#   t_axis = delta_t:delta_t:(nStep-1)*delta_t
#   source_part = sin.(pi * K .* t_axis .* (t_sweep .- t_axis)) .* cos.(2.0 * pi * f0 .* t_axis) ./ (pi*K.*t_axis*t_sweep)
#   for i = 2:nStep
#     source[i] = source_part[i-1]
#   end
#   source[1] = 1.0
#   return source
# end

function computeRsxxzz(Vp, Vs, ind_z, ind_x)
  # since for c++ the index starts from 0
  ind_z .+= 1
  ind_x .+= 1
  Vp_pad = np.pad(Vp, ((4,4),(4,4)), "edge")
  Vs_pad = np.pad(Vs, ((4,4),(4,4)), "edge")
  Vp_pad
  ind_z .+= 4
  ind_x .+= 4
  rxz = zeros(Float64,size(ind_z))
  array_length = length(ind_z)
  Mask = ones(9,9)
  Mask[5,5] = 0.0
  for i = 1:array_length
    vp_ave = mean(Vp_pad[ind_z[i]-4:ind_z[i]+4, ind_x[i]-4:ind_x[i]+4].*Mask)
    vs_ave = mean(Vs_pad[ind_z[i]-4:ind_z[i]+4, ind_x[i]-4:ind_x[i]+4].*Mask)
    # rxz[i] = (vp_ave^2 - vs_ave^2)/(vp_ave^2 - 2*vs_ave^2) #for 2
    rxz[i] = (vp_ave^2)/(vp_ave^2 - 2*vs_ave^2) # for 3
  end
  return rxz
end

function weightObsTraces(obsPath, scratch, scaledPath, nShots, nJump, nt, nx)
  # obsPath: the path to observed data
  # scaledPath: the path to which we save scaled data
  # scratch: where we can find preconditioned observed and synthetic data

  if !isdir(scaledPath)
    mkdir(scaledPath)
  end

  for i=1:nJump:nShots
    S_obs = read(obsPath * "/Shot$(i-1).bin")
    S_obs = reshape(reinterpret(Float32, S_obs), (nt, nx))

    S_condobs = read(scratch * "/CondObs_Shot$(i-1).bin")
    S_condobs = reshape(reinterpret(Float32, S_condobs), (nt, nx))

    S_syn = read(scratch * "/Syn_Shot$(i-1).bin")
    S_syn = reshape(reinterpret(Float32, S_syn), (nt, nx))

    S_fact_obs = zeros(Float32, size(S_obs))

    for iTrace = 1:nx
      fact_nomi = norm(S_syn[:, iTrace])
      fact_denomi = norm(S_condobs[:, iTrace])
      if fact_denomi != 0
        println("ratio = $(fact_nomi/(fact_denomi))")
        S_fact_obs[:, iTrace] = S_obs[:, iTrace] * fact_nomi/(fact_denomi+1e-12)
        # S_fact_obs[:, iTrace] = S_obs[:, iTrace]
      end
    end

    open(scaledPath*"/Shot$(i-1).bin", "w") do f
      write(f, S_fact_obs)
    end
  end

end

ext: jl

max_stars_repo_name: lidongzh/TorchFWI

max_issues_repo_name: lidongzh/TorchFWI

max_forks_repo_name: lidongzh/TorchFWI


text:
using Dates

# mean equatorial circumference, https://en.wikipedia.org/wiki/Earth
const EARTH_CIRCUMFERENCE_METERS = 40_075_017

function parse_gtfsdate(gtfs_date::String)::Date
    return Date(gtfs_date, "yyyymmdd")
end

# TODO automated test
function parse_gtfsdate(gtfs_date::Integer)::Date
    year = gtfs_date ÷ 10_000
    month = (gtfs_date - year * 10_000) ÷ 100
    day = gtfs_date % 100

    return Date(year, month, day)
end

# convert a gtfs time to seconds since midnight
function parse_gtfstime(gtfs_time::String)::Int32
    h, m, s = split(gtfs_time, ':')
    return parse(Int32, h) * 3600 + parse(Int32, m) * 60 + parse(Int32, s)
end

function parse_gtfstime(gtfs_time::Missing)::Int32
    return INT_MISSING
end

# TODO UNIT TESTS
function meters_to_degrees_lat(meters::Real)::Float64
    return meters / EARTH_CIRCUMFERENCE_METERS * 360
end

function meters_to_degrees_lon(meters::Real, lat::Real)::Float64
    return meters_to_degrees_lat(meters) / cosd(lat)
end

function distance_meters(lat1::Real, lon1::Real, lat2::Real, lon2::Real)::Float64
    lat_diff = abs(lat1 - lat2)
    lat_diff_m = lat_diff * EARTH_CIRCUMFERENCE_METERS / 360

    lon_diff = abs(lon1 - lon2)
    lon_diff_m = lon_diff * cosd((lon1 + lon2) / 2) * EARTH_CIRCUMFERENCE_METERS / 360

    return sqrt(lat_diff_m^2 + lon_diff_m^2)
end

# return the indices into destinations that are within the bbox continaing radius_meters around the origin
function bbox_filter(origin::Coordinate, destinations::Vector{Coordinate}, max_dist_meters::Real)::Vector{Int64}
    max_lat_diff = meters_to_degrees_lat(max_dist_meters)
    max_lon_diff = meters_to_degrees_lon(max_dist_meters, origin.latitude)

    candidate_dests = filter(collect(enumerate(destinations))) do t
        idx, dest = t
        return (
            abs(origin.latitude - dest.latitude) < max_lat_diff &&
            abs(origin.longitude - dest.longitude) < max_lon_diff
        )
    end

    return collect(map(c -> c[1], candidate_dests))
end

# primes are used in hashing
# these are from https://primes.utm.edu/lists/small/1000.txt
const primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013,
1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291,
1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583,
1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987,
1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213,
2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,
2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617,
2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819,
2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,
2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181,
3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257,
3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511,
3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571,
3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907,
3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139,
4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231,
4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493,
4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583,
4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831,
4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937,
4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003,
5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179,
5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387,
5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521,
5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639,
5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693,
5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,
5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939,
5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053,
6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,
6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301,
6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,
6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571,
6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673,
6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761,
6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917,
6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997,
7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103,
7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297,
7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411,
7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499,
7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643,
7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723,
7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829,
7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,]


text:
using TermWin

tshow( deleteat! )
tshow( methodswith( Set ) )


text:
using ContextExprRules
using Test

include("test_constraints.jl")
include("operators.jl")
include("test_utils.jl")


text:
immutable Rgb{T}
    r::T
    g::T
    b::T
    Rgb() = new(zero(T), zero(T), zero(T))
    Rgb(r, g, b) = new(r, g, b)
end
Base.zero{T}(::Rgb{T}) = Rgb{T}()
Base.zero{T}(::Type{Rgb{T}}) = Rgb{T}()
Base.convert{T}(::Type{Rgb{T}}, p::Rgb) = Rgb{T}(convert(T, p.r), convert(T, p.g), convert(T, p.b))
Base.show{T <: Integer}(io::IO, p::Rgb{T}) = show(io, "Rgb{$T}($(p.r + 0), $(p.g + 0), $(p.b + 0))")

Base.sum(p::Rgb) = p.r + p.g + p.b

Base.start(p::Rgb) = 1
Base.next{T}(p::Rgb{T}, state::Int) =
    state == 1 ? (p.r, 2) :
    state == 2 ? (p.g, 3) :
    state == 3 ? (p.b, 4) : throw("Iteration state error.")
Base.done(p::Rgb, state::Int) = state == 4
Base.getindex{T}(p::Rgb{T}, n::Int) = n == 1 ? p.r : n == 2 ? p.g : n == 3 ? p.b : throw("Color indexing error.")
Base.size(p::Rgb) = (3,)
Base.size(p::Rgb, n) = n == 1 ? 3 : throw("Invalid dimension.")

Base.rand{T}(::Type{Rgb{T}}) = Rgb{T}(rand(T), rand(T), rand(T))

function sqdiff(x::Rgb, y::Rgb)
    (x.r - y.r) * (x.r - y.r) +
    (x.g - y.g) * (x.g - y.g) +
    (x.b - y.b) * (x.b - y.b)
end

typealias Pixel Rgb{Uint8}
typealias GLPixel Rgb{Uint8} # BUG: Should be Rgb{GLubyte}, but this causes odd dependency chains

-{X, Y}(x::Rgb{X}, y::Rgb{Y}) = Rgb{promote_type(X, Y)}(x.r - y.r, x.g - y.g, x.b - y.b)
+{X, Y}(x::Rgb{X}, y::Rgb{Y}) = Rgb{promote_type(X, Y)}(x.r + y.r, x.g + y.g, x.b + y.b)
/{X, Y}(x::Rgb{X}, y::Rgb{Y}) = Rgb{promote_type(X, Y)}(x.r / y.r, x.g / y.g, x.b / y.b)
*{X, Y}(x::Rgb{X}, y::Rgb{Y}) = Rgb{promote_type(X, Y)}(x.r * y.r, x.g * y.g, x.b * y.b)

# TODO: Could these be sped up by specializing on the type of number?
-{T}(x::Rgb{T}, y::Number) = Rgb{T}(x.r - y, x.g - y, x.b - y)
+{T}(x::Rgb{T}, y::Number) = Rgb{T}(x.r + y, x.g + y, x.b + y)
/{T}(x::Rgb{T}, y::Number) = Rgb{T}(x.r / y, x.g / y, x.b / y)
*{T}(x::Rgb{T}, y::Number) = Rgb{T}(x.r * y, x.g * y, x.b * y)

-{T}(x::Number, y::Rgb{T}) = Rgb{T}(x - y.r, x - y.g, x - y.b)
+{T}(x::Number, y::Rgb{T}) = Rgb{T}(x + y.r, x + y.g, x + y.b)
/{T}(x::Number, y::Rgb{T}) = Rgb{T}(x / y.r, x / y.g, x / y.b)
*{T}(x::Number, y::Rgb{T}) = Rgb{T}(x * y.r, x * y.g, x * y.b)

function hue(c::Pixel)
    c_min = min(c.r, c.g, c.b)
    c_max = max(c.r, c.g, c.b)
    l = (c_max - c_min) / 2

    if c_max == c_min
        return 0.0
    end

    if l < 0.5; s = (c_max - c_min) / (c_max + c_min)
    else;       s = (c_max - c_min) / (2.0 - c_max - c_min)
    end

    if c_max == c.r
        h = (c.g - c.b) / (c_max - c_min)
    elseif c_max == c.g
        h = 2.0 + (c.b - c.r) / (c_max - c_min)
    else
        h = 4.0 + (c.r - c.g) / (c_max - c_min)
    end

    h *= 60
    if h < 0
        h += 360
    elseif h > 360
        h -= 360
    end

    return h
end

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
module Multivectors

import Base:
    (*), (+), (-),
    show, iterate, getindex, length
import StaticArrays:
    FieldVector

export Vector3D, norm²,
    Scalar, ŝ,
    VectorBasis3D, x̂, ŷ, ẑ,
    PseudoVectorBasis3D, x̂ŷ, x̂ẑ, ŷẑ,
    PseudoScalar3D, x̂ŷẑ, ı̂,
    MultiVector3D

include("bases.jl")
include("multivector3d.jl")

end # module

max_stars_repo_name: Samuele-Colombo/Multivectors

max_issues_repo_name: Samuele-Colombo/Multivectors

max_forks_repo_name: Samuele-Colombo/Multivectors


text:
using Distributions, LaTeXStrings, Random, Plots
Random.seed!(1234)

# Simulate n coin flips with a bias of 0.6
n = 1000
flips = rand(n) .< 0.6
p = sum(flips)
q = n - p

# Create plot of Beta(α,β) density
α = 10
β = 5

x_values = 0:0.001:1

# Get prior and posterior values
prior = [pdf(Beta(α,β),x) for x in x_values]
posterior = [pdf(Beta(α+p,β+q),x) for x in x_values]

# Create plots
plot(x_values, posterior, seriestype = :line, label = "posterior",
 yaxis = L"f(\theta|y)", xaxis = L"\theta", grid = false)

plot!(x_values, prior, seriestype = :line, label = "prior", grid = false)

posterior_quantile = cdf(Beta(α+p,β+q),0.65) - cdf(Beta(α+p,β+q),0.55)
confidence_interval = [quantile(Beta(α+p,β+q),0.025), quantile(Beta(α+p,β+q),0.975)]
bayes_mean = mean(Beta(α+p,β+q))

plot(x_values, normal_values)
savefig("beta_prior_coin_ex_2.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


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

using Markdown
using InteractiveUtils

# ╔═╡ e913d16a-8ddd-11eb-094f-3705b4bd1222
begin
using ColorTypes, FileIO, ImageDraw, BenchmarkTools
img = load("test2.png")
fill_color = RGB(0);

function floodfillcolor(img, x, y, current_color, fill_color)
    if (img[y, x] != current_color || img[y, x] == fill_color) return end
    if checkbounds(Bool, img, y, x) img[y, x] = fill_color end
    floodfillcolor(img , x + 1, y, current_color, fill_color)
    floodfillcolor(img , x - 1, y, current_color, fill_color)
    floodfillcolor( img, x, y + 1, current_color, fill_color)
    floodfillcolor( img, x, y - 1, current_color, fill_color)
end
function floodfill(img, x, y, fill_color)
    current_color = img[y,x]
    floodfillcolor(img, x, y, current_color, fill_color)
	img
end
img_new = floodfill(img, 10,100, fill_color); # rectangle boundary
img_case2 = floodfill(img, 51,200, fill_color); # circle boundary
img_case3 = floodfill(img, 350,200, fill_color); # circle fill
img_case4 = floodfill(img, 550,200, fill_color) ;# square fill
end

# ╔═╡ fc39f9d6-8de5-11eb-01f5-ad97913d8b7e

# ╔═╡ Cell order:
# ╠═e913d16a-8ddd-11eb-094f-3705b4bd1222
# ╠═fc39f9d6-8de5-11eb-01f5-ad97913d8b7e

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
export SerialSCPIInstrument

mutable struct SerialSCPIInstrument <: SCPIInstrument
    serial_port::SerialPort
    is_connected::Bool
    delimiter::String

    function SerialSCPIInstrument(port::String,
        baudrate::Integer;
        mode::SPMode=SP_MODE_READ_WRITE,
        ndatabits::Integer=8,
        parity::SPParity=SP_PARITY_NONE,
        nstopbits::Integer=1,
        flow_control::Bool=false)
        inst =  new(open(port, baudrate, mode=mode), true, "\n")

        set_frame(inst.serial_port, ndatabits=ndatabits, parity=parity, nstopbits=nstopbits)

        if flow_control
            set_flow_control(inst.serial_port, rts= SP_RTS_FLOW_CONTROL, cts=SP_CTS_FLOW_CONTROL)
        else
            set_flow_control(inst.serial_port, rts= SP_RTS_OFF, cts=SP_CTS_IGNORE)
        end

        return inst
    end
end

function Base.close(inst::SerialSCPIInstrument)
    close(inst.serial_port)
    inst.is_connected = false
end

"""
Send a command to the instrument
"""
function command(inst::SerialSCPIInstrument, cmd::String)
    write(inst.serial_port, cmd*inst.delimiter)
    return nothing
end

"""
Perform a query to the instrument. Return String.
"""
function query(inst::SerialSCPIInstrument, cmd::String)
    command(inst, cmd)
    return readline(inst.serial_port)[1:end]
end

"""
Perform a query to the instrument. Parse result as type T.
"""
function query(inst::SerialSCPIInstrument, cmd::String, T::Type)
    a = query(inst, cmd)
    return parse(T, a)
end


text:
# Test @assume and @predict macros on a model with conditioning.
# We may want to add comparison of p to the test. REVIEW: what does this comment mean? (Kai)

using Turing
using Distributions
using Test

prior = Beta(2,2)
obs = [0,1,0,1,1,1,1,1,1,1]
exact = Beta(prior.α + sum(obs), prior.β + length(obs) - sum(obs))
meanp = exact.α / (exact.α + exact.β)

@model testbb(obs) = begin
  p ~ Beta(2,2)
  x ~ Bernoulli(p)
  for i = 1:length(obs)
    obs[i] ~ Bernoulli(p)
  end
  p, x
end

s = SMC(10000)
p = PG(100,1000)
g = Gibbs(1500, HMC(1, 0.2, 3, :p), PG(100, 1, :x))

check_numerical(
  sample(testbb(obs), s),
  [:p], [meanp], eps=0.05
)

check_numerical(
  sample(testbb(obs), p),
  [:x], [meanp], eps=0.1
)

check_numerical(
  sample(testbb(obs), g),
  [:x], [meanp], eps=0.1
)


text:
using GalerkinSparseGrids
using Test
using HCubature
using StaticArrays
using SparseArrays

#--------------------------------------
# Constructing the Orthogonal DG Basis
#--------------------------------------

@testset "basis_constructor.jl" begin
    import GalerkinSparseGrids.dg_basis
    @info "Testing that the DG basis is created orthogonally"
    @testset "dg_basis does not throw error" begin
        for k in 1:5
            Q = dg_basis(k)
            @test typeof(Q) == Array{Array{Float64,1},1}
        end
    end
end


text:
using Surrogates
using Test

#1D
lb = 0.0
ub = 5.0
obj_1D = x -> log(x)*exp(x)
x = sample(5,lb,ub,SobolSample())
y = obj_1D.(x)
my_second_order_poly = SecondOrderPolynomialSurrogate(x,y,lb,ub)
val = my_second_order_poly(5.0)
add_point!(my_second_order_poly,5.0,238.86)
add_point!(my_second_order_poly,[6.0,7.0],[722.84,2133.94])

#ND
lb = [0.0,0.0]
ub = [10.0,10.0]
obj_ND = x -> log(x[1])*exp(x[2])
x = sample(10,lb,ub,UniformSample())
y = obj_ND.(x)
my_second_order_poly = SecondOrderPolynomialSurrogate(x,y,lb,ub)
val = my_second_order_poly((5.0,7.0))
add_point!(my_second_order_poly,(5.0,7.0),1764.96)
add_point!(my_second_order_poly,[(1.5,1.5),(3.4,5.4)],[1.817,270.95])

# Multi-output #98
f  = x -> [x^2, x]
lb = 1.0
ub = 10.0
x  = sample(5, lb, ub, SobolSample())
push!(x, 2.0)
y  = f.(x)
surrogate = SecondOrderPolynomialSurrogate(x, y, lb, ub)
# should be exact
@test surrogate.β ≈ [0 0; 0 1; 1 0]

@test surrogate(2.0) ≈ [4, 2]
@test surrogate(1.0) ≈ [1, 1]

f  = x -> [x[1], x[2]^2]
lb = [1.0, 2.0]
ub = [10.0, 8.5]
x  = sample(20, lb, ub, SobolSample())
push!(x, (1.0, 2.0))
y  = f.(x)
surrogate = SecondOrderPolynomialSurrogate(x, y, lb, ub)
@test surrogate.β ≈ [0 0; 1 0; 0 0; 0 0; 0 0; 0 1]
@test surrogate((1.0, 2.0)) ≈ [1, 4]
x_new = (2.0, 2.0)
y_new = f(x_new)
@test surrogate(x_new) ≈ y_new
add_point!(surrogate, x_new, y_new)
@test surrogate(x_new) ≈ y_new


text:
module DynamicPipe

export @>>
export @>

pipe_symbols = [:|>, :.|>]

function get_first_part_of_pipe(pipe)
    global pipe_subsection = pipe
    while pipe_subsection isa Expr
        if pipe_subsection.head == :call && pipe_subsection.args[1] ∈ pipe_symbols
            # previous_pipe_subsection for when the first part is just _ as we can't re write this in-place
            global previous_pipe_subsection = pipe_subsection
            global pipe_subsection = pipe_subsection.args[2]
        else
            break
    end end
    return pipe_subsection, previous_pipe_subsection
end

# function to impute the first part of @>> if it contains underscores
# this allows you to use pipes in pipes @>> 2 |> sum([3, @>> _ |> sqrt(_)]) == sum([3, sqrt(2)])
function rewrite_first_part_of_pipe_macro(pipe::Expr, new_arg::Symbol)
    pipe_subsection = pipe.args[3]
    # if of the form @>> arg begin ... end
    if (pipe.args |> length) == 4 && pipe.args[end].head == :block
        pipe.args[3], contains_underscore =create_internal_function_code(pipe_subsection
                                                                      ,new_arg
                                                                      ,impute_as_first_arg = false)
        return pipe, contains_underscore
    else
        # if of the form @>> block .. end
        if pipe_subsection.head == :block
            pipe_subsection = pipe_subsection.args[2]
            pipe.args[3].args[2], contains_underscore = create_internal_function_code(pipe_subsection
                                                                                ,new_arg
                                                                                ,impute_as_first_arg = false)
            return pipe, contains_underscore
        else
            pipe_subsection, previous_pipe_subsection = get_first_part_of_pipe(pipe_subsection)
            contains_underscore = false
            # rewrite the the first part of the pipe
            if pipe_subsection isa Expr
                new_code, contains_underscore = create_internal_function_code(pipe_subsection,
                                                                              new_arg,
                                                                              impute_as_first_arg = false)
                # may not need to do this now we are not copying the code in rewrite_function_internal
                pipe_subsection.args = new_code.args
            elseif pipe_subsection == :_
                contains_underscore = true
                previous_pipe_subsection.args[2] = new_arg
            end
            return pipe, contains_underscore
        end
    end
end

# function to convert the expressions between the pipes |> ... |>
# replace underscores with the new argument and possible impute the first argument
# if there is no underscore
function create_internal_function_code(code::Expr
                                      ,new_arg::Symbol
                                      ;impute_as_first_arg = true)
    excluded_macros = (Symbol("@>>"), Symbol("@>"))

    contains_underscore = false

    loop_start_point = code.head ∈ [:call, :kw] ? min(2, length(code.args)) :
                       code.head == :macrocall ? min(3, length(code.args)) :
                       1

    new_code = copy(code)
    if new_code.head == :macrocall
        if new_code.args[1] ∈ excluded_macros
            if new_code.args[1] == Symbol("@>>")
                return  rewrite_first_part_of_pipe_macro(new_code, new_arg)
            else # it is a @> macro
                return new_code, false
            end
        end
    end
    # loop through replacing underscores with the new argument
    for (location, arg) in [enumerate(new_code.args)...][loop_start_point:end]
        if arg == :_
            new_code.args[location] =  new_arg
            contains_underscore = true
        elseif arg isa Expr
            if arg.head == :call && arg.args[1] ∈ pipe_symbols
                # new_arg = gensym()
                if arg.args[end] isa Expr
                    if arg.args[end].head == :macrocall
                        if arg.args[end].args[1] == Symbol("@>") # we don't want to rewrite this
                            arg.args[2], had_underscore = create_internal_function_code(arg.args[2], new_arg, impute_as_first_arg = false)
                            contains_underscore = contains_underscore || had_underscore
                            continue
                end end end

                new_code.args[location], had_underscore = form_pipe(arg.args[1]
                                                                ,arg.args[2]
                                                                ,:($new_arg -> $(create_internal_function_code(arg.args[end], new_arg)[1]))
                                                                ,new_arg
                                                                ,rewrite_first_part = true
                                                                ,outer_pipe_arg = new_arg)
                contains_underscore = contains_underscore || had_underscore

            elseif arg.head ≠ :macrocall || (arg.head == :macrocall && arg.args[1] ∉ excluded_macros)
                new_code.args[location], had_underscore = create_internal_function_code(arg, new_arg, impute_as_first_arg = false)
                contains_underscore = had_underscore || contains_underscore
            elseif arg.head == :macrocall && arg.args[1] == Symbol("@>>")
                new_code.args[location], had_underscore = rewrite_first_part_of_pipe_macro(arg, new_arg)
                contains_underscore = contains_underscore || had_underscore
            end
        end
    end

    # if does not contain and we are imputing the first argument
    if impute_as_first_arg && !contains_underscore
        if  new_code.head == :call
            if length(new_code.args) > 1
                new_code.args = [new_code.args[1], new_arg, new_code.args[2:end]...]
            else
                new_code.args = [new_code.args[1], new_arg]
            end
        elseif new_code.head == :macrocall
            if new_code.args[1] ∉ excluded_macros
                if length(new_code.args) > 2
                    new_code.args = [new_code.args[1], new_code.args[2], new_arg, new_code.args[3:end]...]
                else
                    new_code.args = [new_code.args[1], new_code.args[2], new_arg]
                end
            end
        #vectorised function handling
        elseif new_code.head == :.
            if !(new_code.args[2] isa QuoteNode)
                new_code.args[2].args = [new_arg, new_code.args[2].args...]
            end
        end
    end
    return new_code, contains_underscore
end

function create_internal_function_code(code::Symbol, new_arg::Symbol; impute_as_first_arg = true)
    if code == :_
        return new_arg, true
    elseif impute_as_first_arg
        return :($code($new_arg)), false
    else
        return code, false
    end
end

function create_internal_function_code(code, new_arg::Symbol; impute_as_first_arg = true)
    return code, false
end

function create_anonymous_function(code
                                   ,new_arg::Symbol
                                   ;impute_as_first_arg = true)

    new_internal_code, had_underscore = create_internal_function_code(code, new_arg, impute_as_first_arg = impute_as_first_arg)
    return :($new_arg -> $new_internal_code), had_underscore
end

function create_anonymous_function(code::Expr
                                   ,new_arg::Symbol
                                   ;impute_as_first_arg = true)
    # don't do anything if @> macrocall
    if code.head == :macrocall
        if code.args[1]  == Symbol("@>")
            return code, false
        end
    end
    new_internal_code, had_underscore = create_internal_function_code(code, new_arg, impute_as_first_arg = impute_as_first_arg)
    return :($new_arg -> $new_internal_code), had_underscore
end

# function to recursivly convert the internal parts of the pipe, |> ... |>
function form_pipe(pipe_func, first_part, second_part, new_arg; rewrite_first_part = false, outer_pipe_arg = nothing)
    # if we have not yet got to the start of the pipe
    if first_part isa Expr
        if first_part.head == :call
            if first_part.args[1] ∈ pipe_symbols
                return form_pipe(first_part.args[1],
                                first_part.args[2],
                                #Expr(:call, pipe_func, create_anonymous_function(first_part.args[end], new_arg)[1], second_part),
                                :($new_arg -> $(Expr(:call, pipe_func, create_internal_function_code(first_part.args[end], new_arg)[1], second_part))),
                                new_arg,
                                rewrite_first_part = rewrite_first_part,
                                outer_pipe_arg = outer_pipe_arg)
            end
        end
    end
    # we are at the start of the pipe
    if rewrite_first_part
         if outer_pipe_arg |> isnothing
            #Expr(:call, pipe_func, create_anonymous_function(first_part, new_arg)[1], second_part)
            return :($new_arg -> $(Expr(:call, pipe_func, create_internal_function_code(first_part, new_arg)[1], second_part)))
         else
            # we are in forming a pipe in a pipe so we need to know if there is an underscore in the first part
            new_first_part, has_underscore = create_internal_function_code(first_part, outer_pipe_arg, impute_as_first_arg = false)
            return Expr(:call, pipe_func, new_first_part, second_part), has_underscore
         end
    else
        return Expr(:call, pipe_func, first_part, second_part)
    end
end

# function that writes the supplied expression
function rewrite_code(code::Expr; rewrite_first_part = false)
    new_arg = gensym()
    if code.head == :call
        if code.args[1] ∈ pipe_symbols
            return form_pipe(code.args[1],
                            code.args[2],
                            #create_anonymous_function(code.args[end], new_arg)[1],
                            :($new_arg -> $(create_internal_function_code(code.args[end], new_arg)[1])),
                            new_arg,
                            rewrite_first_part = rewrite_first_part,
                            outer_pipe_arg = nothing)
        end
    end
    if !rewrite_first_part
        error("If not acting on block @>> must contain at least one pipe, |>.")
    else
        return create_anonymous_function(code, new_arg)[1]
        # return :($new_arg ->  $create_internal_function_code(code, new_arg)[1]))
    end
end

function rewrite_code(code::Symbol)
    new_arg = gensym()
    return create_anonymous_function(code, new_arg)[1]
    # return :($new_arg -> $create_internal_function_code(code, new_arg)[1]))
end

# convert block of code to a function
function rewrite_block(code; rewrite_first_part = false)
    new_arg = gensym()
    # just get the functions remove the line line nodes
    functions = code.args[map((!(x -> isa(x, LineNumberNode))), code.args)]

    # transform into a pipe so we can use form_pipe
    if length(functions) > 1
        pipe = Expr(:call, :|>, functions[1], functions[2])
        if length(functions) > 2
            for func in functions[3:end]
                pipe = Expr(:call, :|>, pipe, func)
            end
        end

        return form_pipe(pipe.args[1],
                        pipe.args[2],
                        :($new_arg -> $(create_internal_function_code(pipe.args[end], new_arg)[1])),
                        #create_anonymous_function(pipe.args[end], new_arg)[1],
                        new_arg,
                        rewrite_first_part = rewrite_first_part)
    else
        if !rewrite_first_part
            error("If @>> is acting on a block then the block must contain multiple lines.")
        else
            return create_anonymous_function(functions[1], new_arg)[1]
        end
    end
end

"""
    @>(code)

Rewrites code to create an anonymous function that takes one argument. Designed to pipe one object
through multiple functions.

The functions are created using the following rules:
1. Underscores are treated as the function argument. `@> sum(_)` is equivalent to `x -> sum(x)`
2. If the is no underscore in the expression then then the the first argument is imputed,
`@> +(3)` is equivalent to `x -> +(x, 3)`
3. If expression is symbol then it is treated as function so `@> print` is interpreted as
`x -> print(x)`
4. The above 2 rules are applied to expressions separate by the pipe operator, `|>`. Hence
`@> [_, 1, 2] |> sum()` is equivalent to `x -> [x, 1, 2] |> x -> sum(x)`
5. These rules also apply to macros, with the exception of @> and @>>.`
6. The input can also be a `begin end block`. A separate function is created for each line in the
block with result of previous function passed into the next using the re-writing rules.
7. If there is a sequence of pipe characters within a sequence of pipe characters then a new pipe is
is created if the first part of the pipe does not contain an underscore so @> [_, 2, sqrt(36) |> _/2]
is equivalent to x -> [x, 3, sqrt(36) |> y -> y/2].

The macro can also be used in itself. Although will often require the use of brackets to get the
desired effect.

Examples:
```julia-repl
julia> 1 |>
           @>  [_ |> @>(_ + 2), 1, 1] |>
               sum
5
```

```julia-repl
julia> 1 |>
           @>  begin
                [1, 1, _ |> @> +(1, 2)]
                sum
            end
6
```
"""
macro >(code::Symbol)
    return :($(esc(rewrite_code(code))))
end

"""
    @>(code)

Rewrites code to create an anonymous function that takes one argument. Designed to pipe one object
through multiple functions.

The functions are created using the following rules:
1. Underscores are treated as the function argument. `@> sum(_)` is equivalent to `x -> sum(x)`
2. If the is no underscore in the expression then then the the first argument is imputed,
`@> +(3)` is equivalent to `x -> +(x, 3)`
3. If expression is symbol then it is treated as function so `@> print` is interpreted as
`x -> print(x)`
4. The above 2 rules are applied to expressions separate by the pipe operator, `|>`. Hence
`@> [_, 1, 2] |> sum()` is equivalent to `x -> [x, 1, 2] |> x -> sum(x)`
5. These rules also apply to macros, with the exception of @> and @>>.`
6. The input can also be a `begin end block`. A separate function is created for each line in the
block with result of previous function passed into the next using the re-writing rules.
7. If there is a sequence of pipe characters within a sequence of pipe characters then a new pipe is
is created if the first part of the pipe does not contain an underscore so @> [_, 2, sqrt(36) |> _/2]
is equivalent to x -> [x, 3, sqrt(36) |> y -> y/2].

The macro can also be used in itself. Although will often require the use of brackets to get the
desired effect.

Examples:
```julia-repl
julia> 1 |>
           @>  [_ |> @>(_ + 2), 1, 1] |>
               sum
5
```

```julia-repl
julia> 1 |>
           @>  begin
                [1, 1, _ |> @> +(1, 2)]
                sum
            end
6
```
"""
macro >(code::Expr)
    if code.head == :block
        return :($(esc(rewrite_block(code, rewrite_first_part = true))))
    else
        return :($(esc(rewrite_code(code, rewrite_first_part = true))))
    end
end

"""
    @>(code)

Rewrites code to create an anonymous function that takes one argument. Designed to pipe one object
through multiple functions.

The functions are created using the following rules:
1. Underscores are treated as the function argument. `@> sum(_)` is equivalent to `x -> sum(x)`
2. If the is no underscore in the expression then then the the first argument is imputed,
`@> +(3)` is equivalent to `x -> +(x, 3)`
3. If expression is symbol then it is treated as function so `@> print` is interpreted as
`x -> print(x)`
4. The above 2 rules are applied to expressions separate by the pipe operator, `|>`. Hence
`@> [_, 1, 2] |> sum()` is equivalent to `x -> [x, 1, 2] |> x -> sum(x)`
5. These rules also apply to macros, with the exception of @> and @>>.`
6. The input can also be a `begin end block`. A separate function is created for each line in the
block with result of previous function passed into the next using the re-writing rules.
7. If there is a sequence of pipe characters within a sequence of pipe characters then a new pipe is
is created if the first part of the pipe does not contain an underscore so @> [_, 2, sqrt(36) |> _/2]
is equivalent to x -> [x, 3, sqrt(36) |> y -> y/2].

The macro can also be used in itself. Although will often require the use of brackets to get the
desired effect.

Examples:
```julia-repl
julia> 1 |>
           @>  [_ |> @>(_ + 2), 1, 1] |>
               sum
5
```

```julia-repl
julia> 1 |>
           @>  begin
                [1, 1, _ |> @> +(1, 2)]
                sum
            end
6
```
"""
macro >(code)
    error("@> can only act on Symbols and Expresions not $(typeof(code))")
end

"""
    @>>(code)

Rewrites code using the same rewriting rules as @> only the first element of the pipe is not
rewritten and is instead passed through the pipe.

The functions are created using the following rules:
1. Underscores are treated as the function argument. `@>> [1,2,3] |> sum(_)` is equivalent to
`[1,2,3] |> x -> sum(x)`
2. If the is no underscore in the expression then then the the first argument is imputed,
`@>> [1,2,3] |> .+(3)` is equivalent to `[1,2,3] |> x -> .+(x, 3)`
3. If expression is symbol then it is treated as function so `@> "hello world" |> print` is interpreted as
`"hello world" |> x -> print(x)`
4. The above 2 rules are applied to expressions separated by the pipe operator, `|>`. Hence
`@>> 1 |> [_, 1, 2] |> sum()` is equivalent to `1 |> x -> [x, 1, 2] |> x -> sum(x)`
5. These rules also applied to macros, with the exception of @> and @>>. Hence `@>> "hello world" |> @show`
is equivalent to `@>> "hello world" |> @show(_)`
6. The macro can also be used in itself. If the @>> appears within itself or @> and _ is in the first
section of the pipe then it comes from the surrounding context. See the example below.
7. Instead of pipe characters, |>, separating the different functions that to be created a `begin ... end`
block can be used with and separate function is created for each line in the block. See the example below.
7. If there is a sequence of pipe characters within a sequence of pipe characters then a new pipe is
is created if the first part of the pipe does not contain an underscore so @>> 1 |> [_, 2, sqrt(36) |> _/2]
is equivalent to 1 |> x -> [x, 3, sqrt(36) |> y -> y/2].

Examples:
```julia-repl
julia> @>> 1 |>
            [(@>> _ |> +(1, 2)), 1, 1] |>
            sum
6
```

```julia-repl
julia> @>> begin
            1
            [(@>> _ |> _ + 2), 1, 1]
            sum
        end
5
```
"""
macro >>(code::Expr)
    if code.head == :block
        return :($(esc(rewrite_block(code, rewrite_first_part = false))))
    else
        return :($(esc(rewrite_code(code, rewrite_first_part = false))))
    end
end

# below should have the same doc string as >>(code::Expr)
"""
    @>>(code)

Rewrites code using the same rewriting rules as @> only the first element of the pipe is not
rewritten and is instead passed through the pipe.

The functions are created using the following rules:
1. Underscores are treated as the function argument. `@>> [1,2,3] |> sum(_)` is equivalent to
`[1,2,3] |> x -> sum(x)`
2. If the is no underscore in the expression then then the the first argument is imputed,
`@>> [1,2,3] |> .+(3)` is equivalent to `[1,2,3] |> x -> .+(x, 3)`
3. If expression is symbol then it is treated as function so `@> "hello world" |> print` is interpreted as
`"hello world" |> x -> print(x)`
4. The above 2 rules are applied to expressions separated by the pipe operator, `|>`. Hence
`@>> 1 |> [_, 1, 2] |> sum()` is equivalent to `1 |> x -> [x, 1, 2] |> x -> sum(x)`
5. These rules also applied to macros, with the exception of @> and @>>. Hence `@>> "hello world" |> @show`
is equivalent to `@>> "hello world" |> @show(_)`
6. The macro can also be used in itself. If the @>> appears within itself or @> and _ is in the first
section of the pipe then it comes from the surrounding context. See the example below.
7. Instead of pipe characters, |>, separating the different functions that to be created a `begin ... end`
block can be used with and separate function is created for each line in the block. See the example below.
7. If there is a sequence of pipe characters within a sequence of pipe characters then a new pipe is
is created if the first part of the pipe does not contain an underscore so @>> 1 |> [_, 2, sqrt(36) |> _/2]
is equivalent to 1 |> x -> [x, 3, sqrt(36) |> y -> y/2].

Examples:
```julia-repl
julia> @>> 1 |>
            [(@>> _ |> +(1, 2)), 1, 1] |>
            sum
6
```

```julia-repl
julia> @>> begin
            1
            [(@>> _ |> _ + 2), 1, 1]
            sum
        end
5
```
"""
macro >>(code)
    error("@>> can only act on expressions containing a pipe, |>, or a `begin ... end` block not $(typeof(code)).")
end

"""
    @>>(first_arg, pipe::Expr)

Rewrites `pipe` into an anonymous function and passes `first_arg` into it. `pipe` must be a `begin ... end` block.
A function is created for each line in the `begin ... end` block. The result of the previous function is passed
into the next.

The pipe is created using the following rules:
1. Underscores are treated as the function argument. `@>> [1,2,3] begin sum(_) end` is equivalent to
`[1,2,3] |> x -> sum(x)`
2. If the is no underscore in the expression then then the the first argument is imputed,
`@>> [1,2,3] begin .+(3) end` is equivalent to `[1,2,3] |> x -> .+(x, 3)`
3. If expression is a symbol then it is treated as function so `@>> "hello world" begin print end` is interpreted as
`"hello world" |> x -> print(x)`
4. The above 2 rules are applied to each line of a begin end block. Hence
```
    @>> 1 begin
        [_, 1, 2]
        sum()
    end
```
is equivalent to `1 |> x -> [x, 1, 2] |> x -> sum(x)`
5. These rules also applied to macros, with the exception of @> and @>>. Hence `@>> "hello world" begin @show end`
is equivalent to `"hello world" |> x -> @show(x )`
6. The macro can also be used in itself. If the @>> appears within itself or @> and _ is in the first
section of the pipe then it comes from the surrounding context. See the example below.
7. If there is a sequence of pipe characters within a sequence of pipe characters then a new pipe is
is created if the first part of the pipe does not contain an underscore so
```
    @>> 1 begin
        [_, 2, sqrt(36) |> _/2]
    end

````
is equivalent to 1 |> x -> [x, 3, sqrt(36) |> y -> y/2].

Examples:
```julia-repl
julia> @>> 1 begin
            [(@>> _ |> +(2)), 1, 1]
            sum
        end
5
```
"""
macro >>(first_arg, pipe::Expr)
    if pipe.head == :block
        pipe.args = [first_arg, pipe.args...]
        return esc(:(@>> $pipe))
    end
    error("When @>> is acting on 2 arguments the second argument must be a \"begin ... end\" block expression.")
end

# below should have the same docstring as @>>(first_arg, pipe::Expr)
"""
    @>>(first_arg, pipe::Expr)

Rewrites `pipe` into an anonymous function and passes `first_arg` into it. `pipe` must be a `begin ... end` block.
A function is created for each line in the `begin ... end` block. The result of the previous function is passed
into the next.

The pipe is created using the following rules:
1. Underscores are treated as the function argument. `@>> [1,2,3] begin sum(_) end` is equivalent to
`[1,2,3] |> x -> sum(x)`
2. If the is no underscore in the expression then then the the first argument is imputed,
`@>> [1,2,3] begin .+(3) end` is equivalent to `[1,2,3] |> x -> .+(x, 3)`
3. If expression is a symbol then it is treated as function so `@>> "hello world" begin print end` is interpreted as
`"hello world" |> x -> print(x)`
4. The above 2 rules are applied to each line of a begin end block. Hence
```
    @>> 1 begin
        [_, 1, 2]
        sum()
    end
```
is equivalent to `1 |> x -> [x, 1, 2] |> x -> sum(x)`
5. These rules also applied to macros, with the exception of @> and @>>. Hence `@>> "hello world" begin @show end`
is equivalent to `"hello world" |> x -> @show(x )`
6. The macro can also be used in itself. If the @>> appears within itself or @> and _ is in the first
section of the pipe then it comes from the surrounding context. See the example below.
7. If there is a sequence of pipe characters within a sequence of pipe characters then a new pipe is
is created if the first part of the pipe does not contain an underscore so
```
    @>> 1 begin
        [_, 2, sqrt(36) |> _/2]
    end

````
is equivalent to 1 |> x -> [x, 3, sqrt(36) |> y -> y/2].

Examples:
```julia-repl
julia> @>> 1 begin
            [(@>> _ |> +(2)), 1, 1]
            sum
        end
5
```
"""
macro >>(first_arg, second_arg)
    error("When @>> is acting on 2 arguments the second argument must be a \"begin ... end\" block expression not $(typeof(second_arg)). ")
end

end #module

max_stars_repo_name: ConallTimoney/DynamicPipe

max_issues_repo_name: ConallTimoney/DynamicPipes

max_forks_repo_name: ConallTimoney/DynamicPipes


text:
function updatediscretization!(system::CVSystem)
    # determine time shift
    system.solverparams.tshift = system.t[end]-
        sum(system.heart.activation.th[1:end-1]);

    # update total number of time steps
    ntoadd = ceil((system.heart.activation.th[end]-
        system.solverparams.tshift)/system.solverparams.h);
    system.solverparams.numsteps+=ntoadd;

    # update discrete times
    ttoadd = [0:1:ntoadd;]*system.solverparams.h + system.t[end];
    append!(system.t,ttoadd[2:end]);

    # allocate additional space for solution variables
    for i = 1:length(system.branches.ID)
        system.branches.A[i] = [system.branches.A[i] zeros(system.solverparams.JL,length(ttoadd[2:end]))];
        system.branches.Q[i] = [system.branches.Q[i] zeros(system.solverparams.JL,length(ttoadd[2:end]))];
        system.branches.P[i] = [system.branches.P[i] zeros(system.solverparams.JL,length(ttoadd[2:end]))];
        if isempty(system.branches.children[i])
            system.branches.term[i].P = [system.branches.term[i].P;zeros(length(ttoadd[2:end]),5)];
            system.branches.term[i].V = [system.branches.term[i].V;zeros(length(ttoadd[2:end]),5)];
            system.branches.term[i].Q = [system.branches.term[i].Q;zeros(length(ttoadd[2:end]),5)];
        end
    end

    append!(system.svc.P,zeros(length(ttoadd[2:end])));
    append!(system.svc.V,zeros(length(ttoadd[2:end])));
    append!(system.svc.Q,zeros(length(ttoadd[2:end])));
    append!(system.ivc.P,zeros(length(ttoadd[2:end])));
    append!(system.ivc.V,zeros(length(ttoadd[2:end])));
    append!(system.ivc.Q,zeros(length(ttoadd[2:end])));

    append!(system.heart.lv.V,zeros(length(ttoadd[2:end])));
    append!(system.heart.lv.P,zeros(length(ttoadd[2:end])));
    append!(system.heart.lv.E,zeros(length(ttoadd[2:end])));

    append!(system.heart.av.zeta,zeros(length(ttoadd[2:end])));

    append!(system.heart.rv.V,zeros(length(ttoadd[2:end])));
    append!(system.heart.rv.P,zeros(length(ttoadd[2:end])));
    append!(system.heart.rv.E,zeros(length(ttoadd[2:end])));
    append!(system.heart.rv.Q,zeros(length(ttoadd[2:end])));

    append!(system.heart.la.V,zeros(length(ttoadd[2:end])));
    append!(system.heart.la.P,zeros(length(ttoadd[2:end])));
    append!(system.heart.la.Q,zeros(length(ttoadd[2:end])));

    append!(system.heart.ra.V,zeros(length(ttoadd[2:end])));
    append!(system.heart.ra.P,zeros(length(ttoadd[2:end])));
    append!(system.heart.ra.Q,zeros(length(ttoadd[2:end])));

    append!(system.lungs.Pp,zeros(length(ttoadd[2:end])));
    system.lungs.Pa = [system.lungs.Pa;zeros(length(ttoadd[2:end]),3)];
    system.lungs.Va = [system.lungs.Va;zeros(length(ttoadd[2:end]),3)];
    system.lungs.Qa = [system.lungs.Qa;zeros(length(ttoadd[2:end]),3)];
    system.lungs.Pv = [system.lungs.Pv;zeros(length(ttoadd[2:end]),2)];
    system.lungs.Vv = [system.lungs.Vv;zeros(length(ttoadd[2:end]),2)];
    system.lungs.Qv = [system.lungs.Qv;zeros(length(ttoadd[2:end]),2)];

    system.liver.P = [system.liver.P;zeros(length(ttoadd[2:end]),4)];
    system.liver.V = [system.liver.V;zeros(length(ttoadd[2:end]),4)];
    system.liver.Q = [system.liver.Q;zeros(length(ttoadd[2:end]),4)];

    append!(system.cns.H,zeros(length(ttoadd[2:end])));
    append!(system.cns.Emaxlv,zeros(length(ttoadd[2:end])));
    append!(system.cns.Emaxrv,zeros(length(ttoadd[2:end])));
    append!(system.cns.R2L,zeros(length(ttoadd[2:end])));
    append!(system.cns.R3L,zeros(length(ttoadd[2:end])));
    append!(system.cns.R2U,zeros(length(ttoadd[2:end])));
    append!(system.cns.R3U,zeros(length(ttoadd[2:end])));
    append!(system.cns.C4L,zeros(length(ttoadd[2:end])));
    append!(system.cns.C5L,zeros(length(ttoadd[2:end])));
    append!(system.cns.C4U,zeros(length(ttoadd[2:end])));
    append!(system.cns.C5U,zeros(length(ttoadd[2:end])));
    append!(system.cns.V4L,zeros(length(ttoadd[2:end])));
    append!(system.cns.V5L,zeros(length(ttoadd[2:end])));
    append!(system.cns.V4U,zeros(length(ttoadd[2:end])));
    append!(system.cns.V5U,zeros(length(ttoadd[2:end])));

end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Plots;
pyplot();
using LinearAlgebra

function angular_position(points)
    cpoints = points[1] + im * points[2]
    return rad2deg(angle(cpoints))
end

function perturbed_circle(point, meanradius, amplitude, frequency)
    angle = angular_position(point)
    interfaceposition = meanradius + amplitude * sind(angle * frequency)
    distance = interfaceposition - norm(point)
    return distance
end

xrange = -1:1e-1:1
meanradius = 0.5
amplitude = 0.1
frequency = 5

plot(
    xrange,
    xrange,
    (x,y) -> -perturbed_circle([x,y], meanradius, 0., frequency),
    st = :surface,
    #camera = (-30, 30),
)
plot!(xrange,xrange,(x,y)->0.2,st=:surface)


text:
function momintegrals!(op::MWSingleLayer3D, g::RTRefSpace, f::RTRefSpace, t, s, z, strat::BogaertStrategy)

    T, GG  = GetIntegrals(t, s, op.gamma, strat)

    # Get the (u,v,w) → (x,y,z) tf matrix for tcell == bcell
    P = @SMatrix [
        t[1][1] t[2][1] t[3][1]
        t[1][2] t[2][2] t[3][2]
        t[1][3] t[2][3] t[3][3]
    ]

    R = @SMatrix [
        s[1][1] s[2][1] s[3][1]
        s[1][2] s[2][2] s[3][2]
        s[1][3] s[2][3] s[3][3]
    ]

    ∫G = sum(T)

    Q = P*T
    ∫xG = @SVector [
        Q[1,1]+Q[1,2]+Q[1,3],
        Q[2,1]+Q[2,2]+Q[2,3],
        Q[3,1]+Q[3,2]+Q[3,3],
    ]

    Q = R*T.'
    ∫Gy = @SVector [
        Q[1,1]+Q[1,2]+Q[1,3],
        Q[2,1]+Q[2,2]+Q[2,3],
        Q[3,1]+Q[3,2]+Q[3,3]
    ]

    Q = P*T*R.'
    ∫xGy = Q[1,1] + Q[2,2] + Q[3,3]

    c₁ = op.α
    c₂ = op.β

    # Now build the shape-shape interactions from these
    α = 1 / volume(t) / volume(s) / 4
    γ = op.gamma
    for i in 1:3
        a = t[i]
        for j in 1:3
            b = s[j]
            z[i,j] = α*c₁*(∫xGy - dot(a,∫Gy) - dot(b,∫xG) + dot(a,b)*∫G) + 4α*c₂*∫G
        end
    end

end

function momintegrals!(op::MWDoubleLayer3D, g::RTRefSpace, f::RTRefSpace,
    τ, σ, z, strat::BogaertStrategy)

    # Get the primitives
    r = τ.vertices
    s = σ.vertices

    G, GG = GetIntegrals(τ, σ, op.gamma, strat)

    # representation of RT elements in terms of
    # the tangents tu, tv and coordinates u,v,w
    # f_1 = [ (-v-w) t_u + v t_v ] / j
    # f_2 = [ u t_u + (-u-w) t_v ] / j
    # f_3 = [ u t_u + v t_v ] / j

    α = @SMatrix [
        0 -1 -1
        1  0  0
        1  0  0]

    β = @SMatrix [
         0  1  0
        -1  0 -1
         0  1  0]

    duu = τ.tangents[1] × σ.tangents[1]
    duv = τ.tangents[1] × σ.tangents[2]
    dvu = τ.tangents[2] × σ.tangents[1]
    dvv = τ.tangents[2] × σ.tangents[2]

    J = 4 * volume(σ) * volume(τ)
    for i in 1:3
        for j in 1:3
            Iuu = sum([α[i,p] * GG[p,q] * α[j,q] for p in 1:3, q in 1:3])
            Iuv = sum([α[i,p] * GG[p,q] * β[j,q] for p in 1:3, q in 1:3])
            Ivu = sum([β[i,p] * GG[p,q] * α[j,q] for p in 1:3, q in 1:3])
            Ivv = sum([β[i,p] * GG[p,q] * β[j,q] for p in 1:3, q in 1:3])
            z[i,j] = (dot(duu, Iuu) + dot(duv, Iuv) + dot(dvu, Ivu) + dot(dvv, Ivv)) / J
        end
    end

end


text:
# v0.6.0

function magicsquareodd(base::Int)
    if base & 1 == 0 || base < 3; error("base must be odd and >3") end

    square = fill(0, base, base)
    r, number = 1, 1
    size = base * base

    c = div(base, 2) + 1
    while number ≤ size
        square[r, c] = number
        fr = r == 1 ? base : r - 1
        fc = c == base ? 1 : c + 1
        if square[fr, fc] != 0
            fr = r == base ? 1 : r + 1
            fc = c
        end
        r, c = fr, fc
        number += 1
    end

    return square
end

for n in 3:2:7
    println("Magic square with size $n - magic constant = ", div(n ^ 3 + n, 2))
    println("----------------------------------------------------")
    square = magicsquareodd(n)
    for i in 1:n
        println(square[i, :])
    end
    println()
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using StaticArrays

@testset "StaticArrays Size" begin
    @test Size(static(2:4)) === Size(3)
    @test Size(static(3)) === Size(3)
    @test Size(static(3), static(4)) === Size(3, 4)
    @test static(Size(3)) === (static(3),)
    @test typeof(Size()) === Size{()} # Issue #8
    @test typeof(Size(())) === Size{()} # Issue #8
end

@testset "static length indexing" begin
    @test SOneTo(StaticOneTo(3)) === SOneTo(3)
    @test StaticOneTo(SOneTo(3)) === StaticOneTo(3)
    @test StaticArrays.SUnitRange(LengthUnitRange(3:4)) === StaticArrays.SUnitRange(3,4)
    @test LengthUnitRange(StaticArrays.SUnitRange(3,4)) === static(3):static(4)
    @test staticlength(StaticArrays.SUnitRange(3,4)) === static(3):static(4)

    for x in ([1, 2, 3, 4], @MVector([1, 2, 3, 4]), [1 3; 2 4], [1 2; 3 4]', @MMatrix([1 3; 2 4]), @MMatrix([1 2; 3 4])')
        xi = x[static(1):static(2)]
        v = @MVector([1,2])
        @test xi == v
        @test typeof(xi) == typeof(v)
        xi = @stat x[1:2]
        @test xi == v
        @test typeof(xi) == typeof(v)

        # Test of trystatic(lastindex, x)
        xi = @stat x[1:end-2]
        @test xi == v
        if x isa MArray
            @test xi isa MArray
        else
            @test !(xi isa MArray)
        end
    end

    for x in (@SVector([1, 2, 3, 4]), )
        xi = x[static(1):static(2)]
        v = @SVector([1,2])
        @test xi == v
        @test typeof(xi) == typeof(v)

        f(t) = @stat t[2:end-1]
        @test f(x) === SVector((2,3))
        Test.@inferred f(x)
        Test.@inferred f((1,2,3,4))

        @stat g(t,k) = t[k .+ (0:1)]
        @test g(x, 2) === SVector((2,3))
        Test.@inferred g((1,2,3,4), 2)
    end

    @test SVector(@stat 1:2) === SVector{2}([1,2])
    @test SVector(i^2 for i in @stat 1:2) === SVector{2}([1,4])

    f(k) = SVector(@stat 1:k)
    @test f(static(2)) === SVector{2}([1,2])
    Test.@inferred(f(static(2)))

    g(k) = SVector(i^2 for i in @stat 1:k)
    @test g(static(2)) === SVector{2}([1,4])
    Test.@inferred(g(static(2)))

    @test SVector(i^2 for i in static(1):static(3)) === SVector{3}([1, 4, 9])

    @test CleanNamespace.get_one_to(SVector(1, 2, 3), static(2)) === SVector(1, 2)
    @test CleanNamespace.get_one_to(MVector(1, 2, 3), static(2)) == MVector(1, 2)
    @test CleanNamespace.get_one_to(MVector(1, 2, 3), static(2)) isa MVector
end

@inline function testmul(A::SMatrix, B::SMatrix)
    SArray((Base.@_inline_meta; +(Tuple(A[i,k]*B[k,j] for k in static(axes(A,2)))...))
        for i in static(axes(A,1)), j in static(axes(B,2)))
end

@testset "SMatrix generators" begin
    A = SArray(i + j for i = static(2:2:4), j = static(10:10:30))
    @test A === @SMatrix [12 22 32; 14 24 34]
    B = SArray(10i + j for i = static(1:3), j = static(4:5))
    @test B === @SMatrix [14 15; 24 25; 34 35]
    @test testmul(A, B) == [1784  1850; 1928  2000]
    Test.@inferred testmul(A, B)
end

@testset "SMatrix from view" begin
    A = rand(5,5)
    v = view(A, range(2, length=static(3)), range(3, length=static(2)))
    S = Test.@inferred SMatrix(v)
    @test S isa SMatrix
    @test S == A[2:4,3:4]
end


text:
snaps = ["ca-AstroPhA","condmat2005A","loc-Brightkite","EmailEnronA","roadNet-CA","roadNet-TX","web-Google","web-BerkStan","amazon-communities","com-Youtube"]

SizeMat = zeros(6,10)
EdMat = zeros(6,10)
AvgMat = zeros(6,10)
Avg2Mat = zeros(6,10)
MaxMat = zeros(6,10)
RunMat = zeros(6,10)
next = 1
for i = 1:length(snaps)
    global next
    graph = snaps[i]
    F = matread("../data/snap/standard_$graph.mat")
    A = F["A"]
    n = size(A,1)
    M = matread("Output/$(graph)_output.mat")
    C = matread("Output/$(graph)_maxcore_stats.mat")
    mc_stats = C["stats"]  # [length(S) EdgeDen MinDeg MeanDeg Mean2 MaxDeg]

    Sizes = M["Sizes"]
    Times = M["Times"]
    EdgeDen = M["EdgeDens"]
    MeanDeg = M["MeanDegree"]
    Mean2 = M["Mean2Degree"]
    MaxDeg = M["MaxDegree"]
    # MinDeg = M["EdgeDens"]

    # sizes
    s = [mc_stats[1] Sizes[[1 2 5 3 4]]]
    SizeMat[:,next] = s

    # edge densities
    s = [mc_stats[2] EdgeDen[[1 2 5 3 4]]]
    EdMat[:,next] = s

    # mindeg
    # s = [mc_stats[3]; MinDeg[1; 2; 5; 3; 4]]
    # SizeMat[:,next] = s

    # mean degree
    s = [mc_stats[4] MeanDeg[[1 2 5 3 4]]]
    AvgMat[:,next] = s

    # mean 2- degree
    s = [mc_stats[5] Mean2[[1 2 5 3 4]]]
    Avg2Mat[:,next] = s

    # max degree
    s = [mc_stats[6] MaxDeg[[1 2 5 3 4]]]
    MaxMat[:,next] = s

    # time
    s = [Times[2] Times[[1 2 5 3 4]]]
    RunMat[:,next] = s

    next += 1
end

## Try printing this to a Latex Table

print("Size")
G = round.(Int64,SizeMat)
methods = ["\\emph{maxcore}", "\$p = 0.5\$", "\$p = 1.0\$", "\$p = 1.05\$" , "\$p = 1.5\$", "\$p = 2.0\$"]
for i = 1:6
        print(" & $(methods[i])")
    for j = 1:10
        print("& $(G[i,j])")
    end
    print(" \\\\")
    println("")
end

## EDensity
print("Edge Density")
G = round.(EdMat,digits = 3)
methods = ["\\emph{maxcore}", "\$p = 0.5\$", "\$p = 1.0\$", "\$p = 1.05\$" , "\$p = 1.5\$", "\$p = 2.0\$"]
for i = 1:6
        print(" & $(methods[i])")
    for j = 1:10
        print("& $(G[i,j])")
    end
    print(" \\\\")
    println("")
end

## Avg Deg
print("Avg Degree")
G = round.(AvgMat,digits = 2)
methods = ["\\emph{maxcore}", "\$p = 0.5\$", "\$p = 1.0\$", "\$p = 1.05\$" , "\$p = 1.5\$", "\$p = 2.0\$"]
for i = 1:6
        print(" & $(methods[i])")
    for j = 1:10
        print("& $(G[i,j])")
    end
    print(" \\\\")
    println("")
end

## Avg Deg 2
print("Avg")
G = round.(Avg2Mat,digits = 1)
methods = ["\\emph{maxcore}", "\$p = 0.5\$", "\$p = 1.0\$", "\$p = 1.05\$" , "\$p = 1.5\$", "\$p = 2.0\$"]
for i = 1:6
        print(" & $(methods[i])")
    for j = 1:10
        print("& $(G[i,j])")
    end
    print(" \\\\")
    println("")
end

## Max Degree
print("Max")
G = round.(Int64,MaxMat)
methods = ["\\emph{maxcore}", "\$p = 0.5\$", "\$p = 1.0\$", "\$p = 1.05\$" , "\$p = 1.5\$", "\$p = 2.0\$"]
for i = 1:6
        print(" & $(methods[i])")
    for j = 1:10
        print("& $(G[i,j])")
    end
    print(" \\\\")
    println("")
end

## Max Degree
print("Time")
G = round.(RunMat,digits = 2)
methods = ["\\emph{maxcore}", "\$p = 0.5\$", "\$p = 1.0\$", "\$p = 1.05\$" , "\$p = 1.5\$", "\$p = 2.0\$"]
for i = 1:6
        print(" & $(methods[i])")
    for j = 1:10
        print("& $(G[i,j])")
    end
    print(" \\\\")
    println("")
end

## Graph stats

snaps = ["ca-AstroPhA","condmat2005A","loc-Brightkite","EmailEnronA","roadNet-CA","roadNet-TX","web-Google","web-BerkStan","amazon-communities","com-Youtube"]
for i = 1:length(snaps)
    graph = snaps[i]
    F = matread("../data/snap/standard_$graph.mat")
    A = F["A"]
    n = size(A,1)
    m = round(Int64,sum(A/2))
    print("$n $m")
end

max_stars_repo_name: nveldt/GenMeanDSG

max_issues_repo_name: nveldt/GenMeanDSG

max_forks_repo_name: nveldt/GenMeanDSG