text:
"""
A ProjMPO computes and stores the projection of an
MPO into a basis defined by an MPS, leaving a
certain number of site indices of the MPO unprojected.
Which sites are unprojected can be shifted by calling
the `position!` method.

Drawing of the network represented by a ProjMPO `P(H)`,
showing the case of `nsite(P)==2` and `position!(P,psi,4)`
for an MPS `psi`:

```
o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o-      -o--o--o--o--o--o |psi>
```
"""
mutable struct ProjMPO
  lpos::Int
  rpos::Int
  nsite::Int
  H::MPO
  LR::Vector{ITensor}
  ProjMPO(H::MPO) = new(0, length(H) + 1, 2, H, Vector{ITensor}(undef, length(H)))
end

"""
    nsite(P::ProjMPO)

Retrieve the number of unprojected (open)
site indices of the ProjMPO object `P`
"""
nsite(P::ProjMPO) = P.nsite

"""
    length(P::ProjMPO)

The length of a ProjMPO is the same as
the length of the MPO used to construct it
"""
Base.length(P::ProjMPO) = length(P.H)

function lproj(P::ProjMPO)
  (P.lpos <= 0) && return nothing
  return P.LR[P.lpos]
end

function rproj(P::ProjMPO)
  (P.rpos >= length(P) + 1) && return nothing
  return P.LR[P.rpos]
end

"""
    product(P::ProjMPO,v::ITensor)::ITensor

    (P::ProjMPO)(v::ITensor)

Efficiently multiply the ProjMPO `P`
by an ITensor `v` in the sense that the
ProjMPO is a generalized square matrix
or linear operator and `v` is a generalized
vector in the space where it acts. The
returned ITensor will have the same indices
as `v`. The operator overload `P(v)` is
shorthand for `product(P,v)`.
"""
function product(P::ProjMPO, v::ITensor)::ITensor
  Hv = v
  if isnothing(lproj(P))
    if !isnothing(rproj(P))
      Hv *= rproj(P)
    end
    for j in (P.rpos - 1):-1:(P.lpos + 1)
      Hv *= P.H[j]
    end
  else #if lproj exists
    Hv *= lproj(P)
    for j in (P.lpos + 1):(P.rpos - 1)
      Hv *= P.H[j]
    end
    if !isnothing(rproj(P))
      Hv *= rproj(P)
    end
  end
  if order(Hv) != order(v)
    error(
      string(
        "The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, ",
        "this is probably due to an index mismatch.\nCommon reasons for this error: \n",
        "(1) You are trying to multiply the ProjMPO with the $(nsite(P))-site wave-function at the wrong position.\n",
        "(2) `orthognalize!` was called, changing the MPS without updating the ProjMPO.\n\n",
        "P*v inds: $(inds(Hv)) \n\n",
        "v inds: $(inds(v))",
      ),
    )
  end
  return noprime(Hv)
end

(P::ProjMPO)(v::ITensor) = product(P, v)

"""
    eltype(P::ProjMPO)

Deduce the element type (such as Float64
or ComplexF64) of the tensors in the ProjMPO
`P`.
"""
function Base.eltype(P::ProjMPO)
  elT = eltype(P.H[P.lpos + 1])
  for j in (P.lpos + 2):(P.rpos - 1)
    elT = promote_type(elT, eltype(P.H[j]))
  end
  if !isnothing(lproj(P))
    elT = promote_type(elT, eltype(lproj(P)))
  end
  if !isnothing(rproj(P))
    elT = promote_type(elT, eltype(rproj(P)))
  end
  return elT
end

"""
    size(P::ProjMPO)

The size of a ProjMPO are its dimensions
`(d,d)` when viewed as a matrix or linear operator
acting on a space of dimension `d`.

For example, if a ProjMPO maps from a space with
indices `(a,s1,s2,b)` to the space `(a',s1',s2',b')`
then the size is `(d,d)` where
`d = dim(a)*dim(s1)*dim(s1)*dim(b)`
"""
function Base.size(P::ProjMPO)::Tuple{Int,Int}
  d = 1
  if !isnothing(lproj(P))
    for i in inds(lproj(P))
      plev(i) > 0 && (d *= dim(i))
    end
  end
  for j in (P.lpos + 1):(P.rpos - 1)
    for i in inds(P.H[j])
      plev(i) > 0 && (d *= dim(i))
    end
  end
  if !isnothing(rproj(P))
    for i in inds(rproj(P))
      plev(i) > 0 && (d *= dim(i))
    end
  end
  return (d, d)
end

function makeL!(P::ProjMPO, psi::MPS, k::Int)
  while P.lpos < k
    ll = P.lpos
    if ll <= 0
      P.LR[1] = psi[1] * P.H[1] * dag(prime(psi[1]))
      P.lpos = 1
    else
      P.LR[ll + 1] = P.LR[ll] * psi[ll + 1] * P.H[ll + 1] * dag(prime(psi[ll + 1]))
      P.lpos += 1
    end
  end
end

function makeR!(P::ProjMPO, psi::MPS, k::Int)
  N = length(P.H)
  while P.rpos > k
    rl = P.rpos
    if rl >= N + 1
      P.LR[N] = psi[N] * P.H[N] * dag(prime(psi[N]))
      P.rpos = N
    else
      P.LR[rl - 1] = P.LR[rl] * psi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))
      P.rpos -= 1
    end
  end
end

"""
    position!(P::ProjMPO, psi::MPS, pos::Int)

Given an MPS `psi`, shift the projection of the
MPO represented by the ProjMPO `P` such that
the set of unprojected sites begins with site `pos`.
This operation efficiently reuses previous projections
of the MPO on sites that have already been projected.
The MPS `psi` must have compatible bond indices with
the previous projected MPO tensors for this
operation to succeed.
"""
function position!(P::ProjMPO, psi::MPS, pos::Int)
  makeL!(P, psi, pos - 1)
  makeR!(P, psi, pos + nsite(P))

  #These next two lines are needed
  #when moving lproj and rproj backward
  P.lpos = pos - 1
  return P.rpos = pos + nsite(P)
end

"""
    noiseterm(P::ProjMPO,
              phi::ITensor,
              ortho::String)

Return a "noise term" or density matrix perturbation
ITensor as proposed in Phys. Rev. B 72, 180403 for aiding
convergence of DMRG calculations. The ITensor `phi`
is the contracted product of MPS tensors acted on by the
ProjMPO `P`, and `ortho` is a String which can take
the values `"left"` or `"right"` depending on the
sweeping direction of the DMRG calculation.
"""
function noiseterm(P::ProjMPO, phi::ITensor, ortho::String)
  if nsite(P) != 2
    error("noise term only defined for 2-site ProjMPO")
  end

  if ortho == "left"
    AL = P.H[P.lpos + 1]
    if !isnothing(lproj(P))
      AL = lproj(P) * AL
    end
    nt = AL * phi
  elseif ortho == "right"
    AR = P.H[P.rpos - 1]
    if !isnothing(rproj(P))
      AR = AR * rproj(P)
    end
    nt = phi * AR
  else
    error("In noiseterm, got ortho = $ortho, only supports `left` and `right`")
  end
  nt = nt * dag(noprime(nt))

  return nt
end

function checkflux(P::ProjMPO)
  checkflux(P.H)
  for n in length(P.LR)
    if isassigned(P.LR, n)
      checkflux(P.LR[n])
    end
  end
end


text:
# # Operations problems with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Sourabh Dalvi

# ## Introduction

# PowerSimulations.jl supports the construction of Operations problems in power system
# with three part cost bids for each time step. MarketBidCost allows the user to pass a
# time-series of variable cost for energy and ancillary services jointly.
# This example shows how to build a Operations problem with MarketBidCost and how to add
# the time-series data to the devices.

# ## Dependencies
using SIIPExamples

# ### Modeling Packages
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using D3TypeTrees

# ### Data management packages
using Dates
using DataFrames
using TimeSeries

# ### Optimization packages
using Cbc #solver

# ### Data
# This data depends upon the [RTS-GMLC](https://github.com/gridmod/rts-gmlc) dataset. Let's
# download and extract the data.

rts_dir = SIIPExamples.download("https://github.com/GridMod/RTS-GMLC")
rts_src_dir = joinpath(rts_dir, "RTS_Data", "SourceData")
rts_siip_dir = joinpath(rts_dir, "RTS_Data", "FormattedData", "SIIP");

# ### Create a `System` from RTS-GMLC data just like we did in the [parsing tabular data example.](../../notebook/2_PowerSystems_examples/parse_tabulardata.jl)
rawsys = PowerSystems.PowerSystemTableData(
    rts_src_dir,
    100.0,
    joinpath(rts_siip_dir, "user_descriptors.yaml"),
    timeseries_metadata_file = joinpath(rts_siip_dir, "timeseries_pointers.json"),
    generator_mapping_file = joinpath(rts_siip_dir, "generator_mapping.yaml"),
);
sys = System(rawsys; time_series_resolution = Dates.Hour(1));

# ### Creating the Time Series data for Energy bid
MultiDay = collect(
    DateTime("2020-01-01T00:00:00"):Hour(1):(DateTime("2020-01-01T00:00:00") + Hour(8783)),
);

# ### Replacing existing ThreePartCost with MarketBidCost
# Here we replace the existing ThreePartCost with MarketBidCost, and add the energy bid
# time series to the system. The TimeSeriesData that holds the energy bid data can be of any
# type (i.e. `SingleTimeSeries` or `Deterministic`) and bid data should be of type
# `Array{Float64}`,`Array{Tuple{Float64, Float64}}` or `Array{Array{Tuple{Float64,Float64}}`.

for gen in get_components(ThermalGen, sys)
    varcost = get_operation_cost(gen)
    market_bid_cost = MarketBidCost(;
        variable = nothing,
        no_load = get_fixed(varcost),
        start_up = (hot = get_start_up(varcost), warm = 0.0, cold = 0.0),
        shut_down = get_shut_down(varcost),
        ancillary_services = Vector{Service}(),
    )
    set_operation_cost!(gen, market_bid_cost)

    data = TimeArray(MultiDay, repeat([get_cost(get_variable(varcost))], 8784))
    _time_series = SingleTimeSeries("variable_cost", data)
    set_variable_cost!(sys, gen, _time_series)
end

# ### Transforming SingleTimeSeries into Deterministic
horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval)

# In the [OperationsProblem example](../../notebook/3_PowerSimulations_examples/1_operations_problems.ipynb)
# we defined a unit-commitment problem with a copper plate representation of the network.
# Here, we want do define unit-commitment problem  with ThermalMultiStartUnitCommitment
# formulation for thermal device representation.

# For now, let's just choose a standard ACOPF formulation.
uc_template = template_unit_commitment(network = CopperPlatePowerModel)

# Currently  energy budget data isn't stored in the RTS-GMLC dataset.
uc_template.devices[:Generators] =
    DeviceModel(ThermalStandard, ThermalMultiStartUnitCommitment)

solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

# Now we can build a 4-hour economic dispatch / ACOPF problem with the RTS data.
problem = OperationsProblem(
    EconomicDispatchProblem,
    uc_template,
    sys,
    horizon = 4,
    optimizer = solver,
    balance_slack_variables = true,
)

# And solve it ...
solve!(problem)

ext: jl


text:
struct Deplol <: DT.Dependency end

@testset "Tokens" begin

    @testset "Untyped Dependencies" begin
        r = root(UntypedDependency)
        @test isroot(r)
        sent = [
            ("The", 2),
            ("cat", 3),
            ("slept", 0),
            (".", 3)
        ]
        for (i, token) in enumerate(sent)
            dep = UntypedDependency(i, token...)
            @test deprel(dep) == ()
            @test form(dep) == token[1]
            @test id(dep) == i
            @test !isroot(dep)
            @test head(dep) == token[2]
            @test untyped(dep) == typed(dep) == ()
        end
        @test !hashead(noval(UntypedDependency))
    end

    @testset "Typed Dependencies" begin
        r = root(TypedDependency{String})
        @test isroot(r)
        @test isroot(root(TypedDependency))
        sent = [
            ("The", "DT", 2),
            ("cat", "NN", 3),
            ("slept", "VBD", 0),
            (".", ".", 3)
        ]
        for (i, token) in enumerate(sent)
            dep = TypedDependency(i, token...)
            @test form(dep) == token[1]
            @test deprel(dep) == token[2]
            @test head(dep) == token[3]
            @test id(dep) == i
            @test !isroot(dep)
            @test untyped(dep) == ()
            @test typed(dep) == (token[2],)
        end
        @test !hashead(noval(TypedDependency))
    end

    @testset "Errors" begin
        d = Deplol()
        function check_error(fn, d)
            try
                fn(d)
                @test false
            catch err
                @test occursin("not implemented", err.msg)
            end
        end
        check_error(DT.dep, d)
        check_error(DT.deprel, d)
        check_error(DT.form, d)
        check_error(DT.hashead, d)
        check_error(DT.head, d)
        check_error(DT.isroot, d)
        check_error(DT.noval, d)
        check_error(DT.root, Deplol)
        check_error(DT.unk, Deplol)
    end
end


text:
"""
Scale Transform
```
    l = 2.0
    tr = ScaleTransform(l)
```
Multiply every element of the input by `l`
"""
struct ScaleTransform{T<:Real} <: Transform
    s::Vector{T}
end

function ScaleTransform(s::T=1.0) where {T<:Real}
    @check_args(ScaleTransform, s, s > zero(T), "s > 0")
    ScaleTransform{T}([s])
end

set!(t::ScaleTransform,ρ::Real) = t.s .= [ρ]
params(t::ScaleTransform) = first(t.s)
dim(str::ScaleTransform) = 1

transform(t::ScaleTransform,x::AbstractVecOrMat,obsdim::Int=defaultobs) = first(t.s) * x

Base.isequal(t::ScaleTransform,t2::ScaleTransform) = isequal(first(t.s),first(t2.s))


text:
@testset "CDS_scaled_add_tests" begin

#test addition structured matrix
TI=Int64
TF=Float64;
n1=30;
n2=20;
N=n1*n2;
x=randn(TF,N);
y=zeros(TF,N);

#test structured matrix
#comp_grid = compgrid( (25, 25),( n1,n2 ) );
comp_grid = compgrid( (TF(25), TF(25)),( n1,n2 ) );
(TV_OP, AtA_diag, dense, TD_n, banded)=get_TD_operator(comp_grid,"TV",TF);
A=TV_OP'*TV_OP;
A=convert(SparseMatrixCSC{TF,TI},A);
(D_z, AtA_diag, dense, TD_n, banded)=get_TD_operator(comp_grid,"D_z",TF);
B=D_z'*D_z;
B=convert(SparseMatrixCSC{TF,TI},B);

C=A+B;

(R_A,offset_A)=mat2CDS(A)
(R_B,offset_B)=mat2CDS(B)
(R_C,offset_C)=mat2CDS(C)

rho=1.0;
CDS_scaled_add!(R_A,R_B,offset_A,offset_B,rho )

#test if CDS storage of the sum is the same as the sum of the CDS storage matrices
@test R_C==R_A
@test offset_C==offset_A

#test if it still works with a matrix-vector product
Cx_native = C*x;
fill!(y,0)
Cx_CDS = CDS_MVp(size(R_A,1),size(R_A,2),R_A,offset_A,x,y)
@test isapprox(Cx_native,Cx_CDS,rtol=10*eps())

#test random matrix
A=sprandn(100,100,0.01)
B=sprandn(100,100,0.01)
x=randn(100)
y=zeros(100)

A=A+B #make sure A has got all offsets of A and B

(R_A,offset_A)=mat2CDS(A)
(R_B,offset_B)=mat2CDS(B)

C=A.+B;
Cx_native = C*x;

rho=1.0;
CDS_scaled_add!(R_A,R_B,offset_A,offset_B,rho )

fill!(y,0)
Cx_CDS = CDS_MVp(size(R_A,1),size(R_A,2),R_A,offset_A,x,y)
@test isapprox(Cx_native,Cx_CDS,rtol=10*eps())

(R_C,offset_C)=mat2CDS(C)
@test R_C==R_A
@test offset_C==offset_A

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
import ClimateMachine.BalanceLaws:
    # declaration
    vars_state,
    # initialization
    init_state_prognostic!,
    init_state_auxiliary!,
    nodal_init_state_auxiliary!,
    # rhs computation
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
    # boundary conditions
    boundary_conditions,
    boundary_state!

import ClimateMachine.DGMethods.NumericalFluxes:
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder,
    NumericalFluxFirstOrder,
    NumericalFluxSecondOrder,
    RusanovNumericalFlux,
    numerical_boundary_flux_second_order!,
    numerical_flux_second_order!,
    numerical_boundary_flux_first_order!,
    numerical_flux_first_order!,
    normal_boundary_flux_second_order!

using Unitful
"""
    Declaration of state variables

    vars_state returns a NamedTuple of data types.
"""
function vars_state(model::ModelSetup, aux::Auxiliary, T)
    # orientation = model.physics.orientation

    @vars begin
        x::T
        y::T
        z::T
        Φ::T
        ∇Φ::SVector{3, T} # TODO: only needed for the linear model
        # orientation::vars_state(orientation, aux, T)
        ρθ_secondary::T  # stores opposite face for primary (atmospheric import)
        F_ρθ_prescribed::T # stores prescribed flux for secondary (ocean import)
    end
end

function vars_state(::ModelSetup, ::Prognostic, T)
    @vars begin
        ρ::T
        ρu::SVector{3, T}
        ρθ::T
        F_ρθ_accum::T
    end
end

function vars_state(::ModelSetup, ::Gradient, T)
    @vars begin
        ∇ρ::T
        ∇u::SVector{3, T}
        ∇θ::T
    end
end

function vars_state(::ModelSetup, ::GradientFlux, T)
    @vars begin
        μ∇ρ::SVector{3, T}
        ν∇u::SMatrix{3, 3, T, 9}
        κ∇θ::SVector{3, T}
    end
end

"""
    Initialization of state variables

    init_state_xyz! sets up the initial fields within our state variables
    (e.g., prognostic, auxiliary, etc.), however it seems to not initialized
    the gradient flux variables by default.
"""

function nodal_init_state_auxiliary!(m::ModelSetup, state_auxiliary, tmp, geom)
    init_state_auxiliary!(m, m.physics.orientation, state_auxiliary, geom)
    # init_state_auxiliary!(m, m.physics.ref_state, state_auxiliary, geom)

    state_auxiliary.ρθ_secondary = 0
    state_auxiliary.F_ρθ_prescribed = 0
end

function init_state_auxiliary!(::ModelSetup, ::SphericalOrientation, state_auxiliary, geom)
    FT = eltype(state_auxiliary)
    _grav = FT(grav(param_set))
    r = norm(geom.coord)
    state_auxiliary.x = geom.coord[1]
    state_auxiliary.y = geom.coord[2]
    state_auxiliary.z = geom.coord[3]
    state_auxiliary.Φ = _grav * r
    state_auxiliary.∇Φ = _grav * geom.coord / r
end

function init_state_auxiliary!(::ModelSetup, ::FlatOrientation, state_auxiliary, geom)
    FT = eltype(state_auxiliary)
    _grav = FT(grav(param_set))
    @inbounds r = geom.coord[3]
    state_auxiliary.x = geom.coord[1]
    state_auxiliary.y = geom.coord[2]
    state_auxiliary.z = geom.coord[3]
    state_auxiliary.Φ = _grav * r
    state_auxiliary.∇Φ = SVector{3, FT}(0, 0, _grav)
end

function init_state_prognostic!(model::ModelSetup, state::Vars, aux::Vars, localgeo, t)
    x = aux.x
    y = aux.y
    z = aux.z

    parameters = model.parameters
    ic = model.initial_conditions

    state.ρ = ic.ρ(parameters, x, y, z)
    state.ρu = ic.ρu(parameters, x, y, z)
    state.ρθ = ic.ρθ(parameters, x, y, z)

    state.F_ρθ_accum = 0

    return nothing
end

"""
    LHS computations
"""
@inline function flux_first_order!(model::ModelSetup, flux::Grad, state::Vars, aux::Vars, t::Real, direction)
    flux.ρu += calc_pressure(model.physics.eos, state) * I

    calc_advective_flux!(flux, model.physics.advection, state, aux, t)

    return nothing
end

@inline function compute_gradient_argument!(model::ModelSetup, grad::Vars, state::Vars, aux::Vars, t::Real)
    calc_diffusive_flux_argument!(grad, model.physics.diffusion, state, aux, t)

    return nothing
end

@inline function compute_gradient_flux!(model::ModelSetup, gradflux::Vars, grad::Grad, state::Vars, aux::Vars, t::Real)
    calc_diffusive_flux!(gradflux, model.physics.diffusion, grad, state, aux, t)

    return nothing
end

@inline function flux_second_order!(
    model::ModelSetup,
    flux::Grad,
    state::Vars,
    gradflux::Vars,
    ::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρ += gradflux.μ∇ρ
    flux.ρu += gradflux.ν∇u
    flux.ρθ += gradflux.κ∇θ

    return nothing
end

"""
    RHS computations
"""
@inline function source!(model::ModelSetup, source::Vars, state::Vars, gradflux::Vars, aux::Vars, t::Real, direction)
    coriolis = model.physics.coriolis
    gravity = model.physics.gravity
    orientation = model.physics.orientation

    calc_force!(source, coriolis, state, aux, orientation, t)
    calc_force!(source, gravity, state, aux, orientation, t)

    source.F_ρθ_accum = (state.ρθ - aux.ρθ_secondary) * model.parameters.λ_coupler
    return nothing
end

"""
    Boundary conditions
"""
@inline boundary_conditions(model::ModelSetup) = model.boundary_conditions

@inline function boundary_state!(numerical_flux, bc::FluidBC, model::ModelSetup, args...)
    # We need to apply boundary conditions for state variables.
    # This depends on first, second, and high-order
    # operations, hence the diffusion model dependence.
    # TODO!: make work for higher-order diffusion
    diffusion = model.physics.diffusion # defaults to `nothing`
    calc_boundary_state!(numerical_flux, bc.ρu, model, diffusion, args...)
    calc_boundary_state!(numerical_flux, bc.ρθ, model, diffusion, args...)
end

calc_boundary_state!(numerical_flux, bc, model, diffusion, args...) = nothing

@inline vertical_unit_vector(::Orientation, aux) = aux.∇Φ / grav(param_set)
@inline vertical_unit_vector(::NoOrientation, aux) = @SVector [0, 0, 1]

function get_components(csolver)
    mA = csolver.component_list.domainA.component_model
    mB = csolver.component_list.domainB.component_model
    return mA, mB
end

"""
function preB(csolver)
    - saves and regrids the EnergyA couplerfield (i.e. regridded mA.state.ρθ[mA.boundary]) on coupler grid to mB.state.ρθ_secondary[mB.boundary] on the domainB grid
    csolver::CplSolver
"""
function preB(csolver)
    mA, mB = get_components(csolver)
    # Set boundary ρθ used in B to ρθ of A surface at start of coupling cycle.
    mB.rhs.state_auxiliary.ρθ_secondary[mB.boundary] .=
        coupler_get(csolver.coupler, :EnergyA, mB.grid.numerical, DateTime(0), u"J")
    # Set boundary flux accumulator to 0.
    mB.state.F_ρθ_accum .= 0

    idx = varsindex(vars(mA.state), :ρθ)[1]

    @info(
        "preB",
        time = string(csolver.t) * " / " * string(mA.time.finish),
        total_ρθA = weightedsum(mA.state, idx),
        total_ρθB = weightedsum(mB.state, idx),
        total_ρθ = weightedsum(mA.state, idx) + weightedsum(mB.state, idx),
        ρθ_surface_maxA = maximum(mA.state.ρθ[mA.boundary]),
        ρθ_surface_maxB = maximum(mB.state.ρθ[mB.boundary]),
    )

    sphere_sfc_area = 4π * mA.grid.domain.radius^2
    isnothing(csolver.fluxlog) ? nothing :
    csolver.fluxlog.A[csolver.steps] = weightedsum(mA.state, idx) ./ csolver.dt ./ sphere_sfc_area # tracer units / s / m^2
    isnothing(csolver.fluxlog) ? nothing :
    csolver.fluxlog.B[csolver.steps] = weightedsum(mB.state, idx) ./ csolver.dt ./ sphere_sfc_area # tracer units / s / m^2
end

"""
function postB(csolver)
    - updates couplerfield EnergyFluxB with mB.state.F_ρθ_accum[mB.boundary] regridded to the coupler grid, and updates the coupler time
    csolver::CplSolver
"""
function postB(csolver)
    mA, mB = get_components(csolver)
    # Pass exports to "coupler" namespace
    # 1. Save mean θ flux at B boundary during the coupling period
    coupler_put!(
        csolver.coupler,
        :EnergyFluxB,
        mB.state.F_ρθ_accum[mB.boundary] ./ csolver.dt,
        mB.grid.numerical,
        DateTime(0),
        u"J",
    )

    # idx = varsindex(vars(mA.state), :ρθ)[1]

    # @info(
    #     "postatmos",
    #     time = time = csolver.t + csolver.dt,
    #     total_ρθ_atmos = weightedsum(mB.state, idx),
    #     total_ρθ_ocean = weightedsum(mA.state, idx),
    #     total_F_accum = mean(mB.state.F_accum[mB.boundary]) * 1e6 * 1e6,
    #     total_ρθ =
    #         weightedsum(mB.state, idx) +
    #         weightedsum(mA.state, idx) +
    #         mean(mB.state.F_ρθ_accum[mB.boundary]) * 1e6 * 1e6,
    #     F_accum_max = maximum(mB.state.F_accum[mB.boundary]),
    #     F_avg_max = maximum(mB.state.F_accum[mB.boundary] ./ csolver.dt),
    #     ρθ_surface_maxB = maximum(mB.state.ρθ[mB.boundary]),
    #     ρθ_surface_maxA = maximum(mA.state.ρθ[mA.boundary]),
    # )
end

"""
function preA(csolver)
    - saves and regrids the EnergyFluxB couplerfield (i.e. regridded mB.state.F_ρθ_accum[mB.boundary]) to ρθ_prescribed[mA.boundary] on the domainA grid
    csolver::CplSolver
"""
function preA(csolver)
    mA, mB = get_components(csolver)

    mA.rhs.state_auxiliary.F_ρθ_prescribed[mA.boundary] .=
        coupler_get(csolver.coupler, :EnergyFluxB, mA.grid, DateTime(0), u"J")
    # Set boundary flux accumulator to 0. (this isn't used)
    mA.state.F_ρθ_accum .= 0

    # @info(
    #     "preA",
    #     time = csolver.t,
    #     F_prescribed_max =
    #         maximum(mA.discretization.state_auxiliary.F_prescribed[mA.boundary]),
    #     F_prescribed_min =
    #         maximum(mA.discretization.state_auxiliary.F_prescribed[mA.boundary]),
    #     ρθ_surface_maxA = maximum(mA.state.ρθ[mA.boundary]),
    #     ρθ_surface_minA = maximum(mA.state.ρθ[mA.boundary]),
    # )
end

"""
function postA(csolver)
    - updates couplerfield EnergyA with mA.state.ρθ[mA.boundary] regridded to the coupler grid, and updates the coupler time
    - updates EnergyA with mA.state.ρθ[mA.boundary], and the coupler time
    csolver::CplSolver
"""
function postA(csolver)
    mA, mB = get_components(csolver)

    # @info(
    #     "postA",
    #     time = csolver.t + csolver.dt,
    #     ρθ_surface_maxA = maximum(mA.state.ρθ[mA.boundary]),
    #     ρθ_surface_minA = maximum(mA.state.ρθ[mA.boundary]),
    # )

    # Pass exports to "coupler" namespace
    coupler_put!(csolver.coupler, :EnergyA, mA.state.ρθ[mA.boundary], mA.grid.numerical, DateTime(0), u"J")
end

abstract type AbstractCouplerBoundary end

export PenaltyNumFluxDiffusive, CoupledPrimaryBoundary, CoupledSecondaryBoundary

"""
struct PenaltyNumFluxDiffusive <: NumericalFluxSecondOrder end
- PenaltyNumFluxDiffusive: additional NumericalFluxSecondOrder to account for an additional penalty term
"""
struct PenaltyNumFluxDiffusive <: NumericalFluxSecondOrder end

# ## 2. CoupledPrimaryBoundary
"""
struct CoupledPrimaryBoundary <: AbstractCouplerBoundary end
# # compute flux based on opposite face
# # also need to accumulate net flux across boundary
"""
struct CoupledPrimaryBoundary <: AbstractCouplerBoundary end

## 3. CoupledSecondaryBoundary
"""
struct CoupledSecondaryBoundary <: AbstractCouplerBoundary end
# # use prescribed flux computed in primary
"""
struct CoupledSecondaryBoundary <: AbstractCouplerBoundary end

function numerical_boundary_flux_second_order!(
    numerical_flux::Union{PenaltyNumFluxDiffusive},
    bctype::FluidBC,
    balance_law::ModelSetup,
    args...,
) where {S, D, A, HD}
    normal_boundary_flux_second_order!( # for ρu BCs (to be implemented in boundary_state!)
        numerical_flux,
        bctype,
        balance_law,
        args...,
    )

    numerical_boundary_flux_second_order!( # this overwrites the total flux for ρθ (using methods below)
        numerical_flux,
        bctype.ρθ,
        balance_law,
        args...,
    )

end

# customized flux for CoupledPrimaryBoundary
function numerical_boundary_flux_second_order!(
    numerical_flux::Union{PenaltyNumFluxDiffusive},
    bctype::CoupledPrimaryBoundary,
    balance_law::ModelSetup,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}

    fluxᵀn.ρθ = (state_prognostic⁻.ρθ - state_auxiliary⁺.ρθ_secondary) * balance_law.parameters.λ_coupler # W/m^2

end

# customized flux for CoupledSecondaryBoundary
function numerical_boundary_flux_second_order!(
    numerical_flux::Union{PenaltyNumFluxDiffusive},
    bctype::CoupledSecondaryBoundary,
    balance_law::ModelSetup,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, A, HD}

    fluxᵀn.ρθ = -state_auxiliary⁺.F_ρθ_prescribed

end

# customized flux for CoupledSecondaryBoundary
function numerical_boundary_flux_second_order!(
    numerical_flux::PenaltyNumFluxDiffusive,
    bctype::Insulating,
    balance_law,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, HD, A}
    fluxᵀn.ρθ = -Float64(0)
    return nothing
end

"""
function numerical_flux_second_order!(::PenaltyNumFluxDiffusive,
  Penalty flux formulation of second order numerical flux. This formulation
  computes the CentralNumericalFluxSecondOrder term first (which is just the average
  of the + and - fluxes and an edge), and then adds a "penalty" flux that relaxes
  the edge state + and - toward each other.
"""
function numerical_flux_second_order!(
    ::PenaltyNumFluxDiffusive,
    bl::ModelSetup,
    fluxᵀn::Vars{S},
    n::SVector,
    state⁻::Vars{S},
    diff⁻::Vars{D},
    hyperdiff⁻::Vars{HD},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    diff⁺::Vars{D},
    hyperdiff⁺::Vars{HD},
    aux⁺::Vars{A},
    t,
) where {S, HD, D, A}

    numerical_flux_second_order!(
        CentralNumericalFluxSecondOrder(),
        bl,
        fluxᵀn,
        n,
        state⁻,
        diff⁻,
        hyperdiff⁻,
        aux⁻,
        state⁺,
        diff⁺,
        hyperdiff⁺,
        aux⁺,
        t,
    )

    Fᵀn = parent(fluxᵀn)
    FT = eltype(Fᵀn)
    tau = FT(0.0) # decide if want to use penalty
    Fᵀn .+= tau * (parent(state⁻) - parent(state⁺))
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
__precompile__(true)

module Simplices
using LinearAlgebra
using Printf

using Distributions

include("simplexsplitting/SimplexSplitting.jl")
include("barycentric-coordinates.jl")
include("Binary.jl")
include("Circumsphere.jl")
include("complementary.jl")
include("embed.jl")
include("geometry.jl")
include("heaviside.jl")
include("intersection-of-boundaries.jl")
include("intersection-of-boundaries-loop.jl")
include("init_functions.jl")
include("NullSpace.jl")
include("polytope-generating-vertices.jl")
include("QR.jl")
include("shared-vertices.jl")
include("sharing-a-face.jl")
include("some-vertex-in-circumsphere.jl")
include("simplexoperations.jl")
include("simplex_sampling.jl")
include("TriangulationPolytopeFaces.jl")
include("TriangulationNonSimplicialFaces.jl")
include("Update.jl")
include("volume-computation.jl")

include("simplexintersect.jl")
include("intersecting_vertices.jl")

export BarycentricCoordinates,
        centroid,
        Circumsphere,
        childsimplex,
        childpoint,
        childpoints,
        Delaunay,
        even_sampling, even_sampling_rules, evenly_sample,
        heaviside,
        heaviside0,
	insidepoints,
        intersecting_simplices,
        intersectingvertices,
        IntersectionOfBoundaries,
        IntersectionOfBoundaries_NoStorage,
        issingular,
        nontrivially_intersecting_simplices,
        nontrivial_intersection,
        outsidepoint,
        orientation,
        PolytopeGeneratingVertices,
        radius,
        SharedVertices,
        SharedFaceVolume,
        SharedFaceVertices,
        shared_vertex_intersection,
	SimplexSplitting,
        simplices_sharing_vertices,
        simplexintersection,
        SomeVertexInCircumsphere,
        tensordecomp,
        volume,
        VolumeComputation,
	outsidepoint,
	outsidepoints

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1KeyToPath <: SwaggerModel
    key::Any # spec type: Union{ Nothing, String } # spec name: key
    mode::Any # spec type: Union{ Nothing, Int32 } # spec name: mode
    path::Any # spec type: Union{ Nothing, String } # spec name: path

    function IoK8sApiCoreV1KeyToPath(;key=nothing, mode=nothing, path=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1KeyToPath, Symbol("key"), key)
        setfield!(o, Symbol("key"), key)
        validate_property(IoK8sApiCoreV1KeyToPath, Symbol("mode"), mode)
        setfield!(o, Symbol("mode"), mode)
        validate_property(IoK8sApiCoreV1KeyToPath, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        o
    end
end # type IoK8sApiCoreV1KeyToPath

const _property_map_IoK8sApiCoreV1KeyToPath = Dict{Symbol,Symbol}(Symbol("key")=>Symbol("key"), Symbol("mode")=>Symbol("mode"), Symbol("path")=>Symbol("path"))
const _property_types_IoK8sApiCoreV1KeyToPath = Dict{Symbol,String}(Symbol("key")=>"String", Symbol("mode")=>"Int32", Symbol("path")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1KeyToPath }) = collect(keys(_property_map_IoK8sApiCoreV1KeyToPath))
Swagger.property_type(::Type{ IoK8sApiCoreV1KeyToPath }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1KeyToPath[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1KeyToPath }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1KeyToPath[property_name]

function check_required(o::IoK8sApiCoreV1KeyToPath)
    (getproperty(o, Symbol("key")) === nothing) && (return false)
    (getproperty(o, Symbol("path")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1KeyToPath }, name::Symbol, val)
end


text:
"""
`Basics.extractNoOpenShells(conf::Configuration)`
    ... determine the number of open (nonrelativistic) shells in conf; a singleton of type AbstractOpenShell is returned.
"""
function Basics.extractNoOpenShells(conf::Configuration)
    ##x println("conf = $conf")
    ns = 0
    for  (shell, occ)  in conf.shells
        if      occ == 0  ||  occ == 2*(2*shell.l + 1)
        else    ns = ns + 1
        end
    end
    ##x println("conf = $conf  ns = $ns")

    if      ns == 0   return ( LSjj.ZeroOpenShell() )
    elseif  ns == 1   return ( LSjj.OneOpenShell() )
    elseif  ns == 2   return ( LSjj.TwoOpenShells() )
    elseif  ns == 3   return ( LSjj.ThreeOpenShells() )
    else    error("stop a")
    end

end

"""
`Basics.extractNonrelativisticShellList(subshells::Array{Subshell,1})
    ... extract a nonrelativistic shell list from the given subshell list; the shells appear in the same order
        as the corresponding subshells. A shellList::Array{Subshell,1} is returned.
"""
function Basics.extractNonrelativisticShellList(subshells::Array{Subshell,1})
    shellList = Shell[]

    for  subsh  in subshells
        l = Basics.subshell_l(subsh);   shell = Shell( subsh.n, l)
        if  shell in shellList      else  push!(shellList, shell)   end
    end
    ##x println("Basics.extractNonrelativisticShellList()::  shellList = $shellList")

    return( shellList )
end

"""
`Basics.extractNonrelativisticConfigurations(basis::Basis)
    ... extract all nonrelativistic configurations that contribute to the given set of CSF in basis.csfs.
        A confList::Array{Configuration,1} is returned.
"""
function Basics.extractNonrelativisticConfigurations(basis::Basis)
    confList  = Configuration[]
    subshells = basis.subshells

    for  csf  in basis.csfs
        newShells = Dict{Shell,Int64}();    NoElectrons = 0
        for  i = 1:length(subshells)
            n = subshells[i].n;    l = Basics.subshell_l(subshells[i]);    occ = csf.occupation[i]
            shell = Shell(n,l)
            if    haskey(newShells, shell)    newShells[shell] = newShells[shell] + occ;   NoElectrons = NoElectrons + occ
            else
                if   occ > 0  newShells = Base.merge( newShells, Dict( shell => occ));     NoElectrons = NoElectrons + occ   end
            end
        end
        if  basis.NoElectrons != NoElectrons    error("stop a")
        else
            conf = Configuration(newShells, NoElectrons)
            if  conf in confList    continue;    else    push!( confList,  conf)    end
        end
    end
    ##x println("Basics.extractNonrelativisticConfigurations()::  confList = $confList")

    return( confList )
end

"""
`Basics.extractNonrelativisticConfigurationFromCsfR(csf::CsfR,  basis::Basis)`
    ... extract the nonrelativistic configuration from the given CSF that is defined in basis.csfs.
        A conf::Configuration is returned.
"""
function  Basics.extractNonrelativisticConfigurationFromCsfR(csf::CsfR,  basis::Basis)
    subshells = basis.subshells

    newShells = Dict{Shell,Int64}();    NoElectrons = 0
    for  i = 1:length(subshells)
        n = subshells[i].n;    l = Basics.subshell_l(subshells[i]);    occ = csf.occupation[i]
        shell = Shell(n,l)
        if    haskey(newShells, shell)    newShells[shell] = newShells[shell] + occ;   NoElectrons = NoElectrons + occ
        else
            if   occ > 0  newShells = Base.merge( newShells, Dict( shell => occ));     NoElectrons = NoElectrons + occ   end
        end
    end

    if      basis.NoElectrons != NoElectrons    error("stop a")
    else    conf = Configuration(newShells, NoElectrons)
    end
    ##x println("Basics.extractNonrelativisticConfigurationFromCsfR()::  conf = $conf")

    return( conf )
end

"""
`Basics.extractOpenShellQNfromCsfR(csfR::CsfR, basis::Basis)`
    ... extracts the quantum numbers of the nonrelativistic open shells from the given csfr. Here, we assume that csfR is
        defined in the given basis, and it is checked that the subshells of each nonrelativistic shell occur subsequently
        in basis and in the right sequence j = l - 1/2, j = l + 1/2, respectively.
        An rQN::Dict{Shell, Tuple{Tuple{Int64,Int64,Int64,Int64,Int64}, Tuple{Int64,Int64,Int64,Int64,Int64}}}
                = Dict(shell => (rQNm, rQNp) ) is returned with rQNm = (idxm, Nm, Qm, Jm, Jmx) and
        rQNp = (idxp, Np, Qp, Jp, Jpx). The two tuples rQNm, rQNp contains all the subshell quantum numbers and their
        coupling in csfR.
"""
function Basics.extractOpenShellQNfromCsfR(csfR::CsfR, basis::Basis)
    wa = Dict{Shell,Tuple{NTuple{5,Int64},NTuple{5,Int64}}}()
    for  s = 1:length(basis.subshells)
        sx = s + 1;     subsh = basis.subshells[s];     n = subsh.n;        occ   = csfR.occupation[s]
        l  = Basics.subshell_l(subsh);    j2 = Basics.subshell_2j(subsh);   shell = Shell(n,l)
        # Get l,j from next subshell ... if still defined
        if  sx <= length(basis.subshells)       subshx = basis.subshells[sx]
            nx = subshx.n;    lx = Basics.subshell_l(subshx);   j2x = Basics.subshell_2j(subshx);   occx   = csfR.occupation[sx]
        end
        # Now collect the quantum numbers into tuples
        if  l == 0   &&   0 < occ < 2
            wa = Base.merge( wa, Dict( shell => ( (s, csfR.occupation[s], csfR.seniority[s],
                                                   Basics.twice(csfR.subshellJ[s]), Basics.twice(csfR.subshellX[s]) ),
                                                  (-9, -9, -9, -9, -9) ) ) )
        elseif  2l + 1 == j2    continue
        elseif  2l - 1 == j2    &&    (n != nx  ||  l != lx  ||  j2x != j2 + 2)
                error("stop a")
        elseif  0 < occ + occx < 2* (2l + 1)
            wa = Base.merge( wa, Dict(shell => ( (s, csfR.occupation[s], csfR.seniority[s],
                                                  Basics.twice(csfR.subshellJ[s]), Basics.twice(csfR.subshellX[s]) ),
                                                 (sx, csfR.occupation[sx], csfR.seniority[sx],
                                                  Basics.twice(csfR.subshellJ[sx]), Basics.twice(csfR.subshellX[sx]) ) ) ) )
        end
    end

    ##x println("Basics.extractOpenShellQNfromCsfR()::  csfR = $csfR  \n  wa = $wa")

    return( wa )
end

"""
`Basics.extractShellOccupationFromCsfR(csfR::CsfR, basis::Basis)`
    ... extracts the nonrelativistic shell list as well as the (nonrelativistic) occupation of the given CsfR.
        Here, we assume that csfR is defined in the given basis and that all subshells occur in standard sequence
        for a jjJ-LSJ transformation. A tuple tp::Tuple{Array{Shell,1},Array{Int64,1}}  of two arrays of equal
        length is returned that contains the nonrelativistic shells and their occupation in the specification of
        csfR.
"""
function Basics.extractShellOccupationFromCsfR(csfR::CsfR, basis::Basis)
    shellList = Basics.extractNonrelativisticShellList(basis.subshells)
    occList   = Int64[]

    for  s = 1:length(basis.subshells)
        l = Basics.subshell_l(basis.subshells[s]);   j2 = Basics.subshell_2j(basis.subshells[s])
        if      l == 0          push!( occList, csfR.occupation[s])
        elseif  2*l + 1 == j2   continue    # contributions from j = l + 1/2 are collected together with j = l - 1/2 below
        else    push!( occList, csfR.occupation[s] + csfR.occupation[s+1])
        end
    end
    ##x println("Basics.extractShellOccupationFromCsfR()::  shellList = $shellList,  occList = $occList")

    return( (shellList, occList) )
end


text:
# A Tables.jl interface for AbstractModel

Tables.istable(::AbstractModel) = true
Tables.columnnames(m::AbstractModel) = keys(m)
Tables.schema(m::AbstractModel) = Tables.Schema(Tables.columnnames(m), _columntypes(m))

_columntypes(m) = map(k -> Union{map(typeof, getindex(m, k))...}, keys(m))

# As Columns
Tables.columnaccess(::Type{<:AbstractModel}) = true
Tables.columns(m::AbstractModel) = m
Tables.getcolumn(m::AbstractModel, nm::Symbol) = collect(getindex(m, nm))
Tables.getcolumn(m::AbstractModel, i::Int) = collect(getindex(m, i))
Tables.getcolumn(m::AbstractModel, ::Type{T}, col::Int, nm::Symbol) where T =
    collect(getindex(m, nm))

# As rows
Tables.rowaccess(::Type{<:AbstractModel}) = true
# Vector of NamedTuple already has a row interface defined,
# so we take a shortcut and return that.
Tables.rows(m::AbstractModel) = [nt for nt in map((p,typ,name) -> merge((component=typ,fieldname=name), parent(p)), params(m), paramparenttypes(m), paramfieldnames(m))]


text:
# Estimates the basic Nerlove Cobb-Douglas model
using DelimitedFiles
function main()
data = readdlm("nerlove.data")
data = data[:,2:6]
data = log.(data)
n = size(data,1)
y = data[:,1]
x = data[:,2:end]
x = [ones(n,1) x]
names = ["constant", "output", "labor", "fuel", "capital"]
b, junk, junk, junk = ols(y, x, names=names)

# First HOD1
R = [0 0 1 1 1]
r = 1
println()
println("Imposing and testing HOD1")
ols(y, x, R=R, r=r, names=names)
TestStatistics(y, x, R, r)

# Now CRTS
R = [0 1 0 0 0]
r = 1
println()
println("Imposing and testing CRTS")
ols(y, x, R=R, r=r, names=names)
TestStatistics(y, x, R, r)
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using Distributions
using Plots

function calc_midpoints(edges::AbstractVector)::Vector{Float64}
    return Float64[0.5 * (edges[i] + edges[i+1]) for i ∈ 1:length(edges)-1]
end

function bin_data1d(data::AbstractVector, bins::Integer)
    low, high = extrema(data)
    edges = collect(range(low, stop = high, length = bins + 1))
    midpoints = calc_midpoints(edges)
    buckets = Int[max(2, min(searchsortedfirst(edges, x), length(edges))) - 1 for x ∈ data]
    counts = zeros(Int, length(midpoints))
    for b in buckets
        counts[b] += 1
    end
    return edges, midpoints, buckets, counts
end

function bin_data2d(data::AbstractMatrix, xbins::Integer, ybins::Integer)
    xmin, ymin = minimum(data, dims = 2)
    xmax, ymax = maximum(data, dims = 2)
    xedges = collect(range(xmin, stop = xmax, length = xbins + 1))
    yedges = collect(range(ymin, stop = ymax, length = ybins + 1))
    xmidpts = calc_midpoints(xedges)
    ymidpts = calc_midpoints(yedges)
    buckets = []
    d, c = size(data)
    for i = 1:c
        x = data[1, i]
        y = data[2, i]
        xb = max(2, min(searchsortedfirst(xedges, x), length(xedges))) - 1
        yb = max(2, min(searchsortedfirst(yedges, y), length(yedges))) - 1
        push!(buckets, (yb, xb))
    end
    counts = zeros(Int, length(ymidpts), length(xmidpts))
    for (yb, xb) in buckets
        counts[yb, xb] += 1
    end
    xedges, yedges, xmidpts, ymidpts, buckets, counts
end

μ = [
    -1
    -2
]

Σ = [
    1.0 -0.7
    -0.7 1.0
]

using Distributions
d = MvNormal(μ, Σ)
xydata = rand(d, 10000)

p1 = plot(xlabel = "x", ylabel = "y")
xedges, yedges, xmidpts, ymidpts, buckets, counts = bin_data2d(xydata, 30, 30)
r, c = size(counts)
scatter!(p1, xydata[1, :], xydata[2, :], alpha = 0.2)
heatmap!(p1, xmidpts, ymidpts, counts, alpha = 0.85)
xind = 18
plot!(
    p1,
    [xmidpts[xind], xmidpts[xind]],
    [minimum(ymidpts), maximum(ymidpts)],
    color = :red,
    legend = :false,
)
p2 = plot(
    ymidpts,
    counts[:, xind] / sum(counts[:, xind]),
    label = "p(y|x=$(xmidpts[xind])",
)
plot(p1, p2, layout = (2, 1))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function add_drag(w, range, point_id, slider_length, slideridx_s)
    m2id = mouse2id(w)
    # interaction
    @materialize mouse_buttons_pressed, mouseposition = w.inputs
    isoverpoint = const_lift(is_same_id, m2id, point_id)
    # single left mousekey pressed (while no other mouse key is pressed)
    key_pressed = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    # dragg while key_pressed. Drag only starts if isoverpoint is true
    mousedragg = GLAbstraction.dragged(mouseposition, key_pressed, isoverpoint)
    preserve(foldp(value(slideridx_s), mousedragg) do v0, dragg
        if dragg != Vec2f0(0)
            idx_steps = round(Int, (dragg[1]/slider_length)*length(range))
            new_idx = clamp(v0 + idx_steps, 1, length(range))
            push!(slideridx_s, new_idx)
            return v0
        else # dragging started
            return value(slideridx_s)
        end
    end)
    slideridx_s
end
function add_play(slideridx_s, play_signal, range, rate=30.0)
    play_s = fpswhen(play_signal, rate)
    preserve(map(play_s, init=nothing) do t
        push!(slideridx_s, mod1(value(slideridx_s)+1, length(range)))
        nothing
    end)
end

function maxdigits(range)
    if eltype(range) <: AbstractFloat
        return 7
    else
        return max(ndigits(first(range)), ndigits(last(range)))
    end
end
function play_widget(
        range, window = GLPlot.widget_screen!();
        startidx::Int = 1
    )
    numberbox = map(icon_size()) do is
        AABB(Vec3f0(0, -1mm, 0), Vec3f0(2is, is, 1))
    end
    target_s = map(x-> x * 0.4, icon_size())
    sliderlen = 70mm - (3 * value(icon_size())) - 2mm
    play_button, play_stop_signal = GLVisualize.toggle_button(
        rot180(GLPlot.imload("play.png")), GLPlot.imload("break.png"), window
    )
    play_s = map(!, play_stop_signal)
    slider_w, slider_s = GLVisualize.slider(range, window,
        startidx = startidx, play_signal = play_s,
        slider_length = sliderlen
    )
    number = visualize(
        map(GLVisualize.printforslider, slider_s),
        color = RGBA{Float32}(0.6, 0.6, 0.6,1),
        boundingbox = numberbox,
        relative_scale = target_s
    )
    add_widget!(play_button, number, slider_w, window=window)

    slider_s
end
export play_widget
export add_widget!
function add_widget!(widgets...;
        delete = Signal(false),
        window = widget_screen!(delete = delete),
        height = value(icon_size())
    )
    scalings = map(widgets) do widget
        bb = value(boundingbox(widget))
        w = widths(bb)
        if w[2] > height # only scale when to big
            s = height/w[2]
            w.*s, minimum(bb)
        else
            w, minimum(bb)
        end
    end
    last_x = 0f0
    foreach(zip(scalings, widgets)) do scaling_widget
        scaling, widget = scaling_widget
        scale, offset = scaling
        xwidth = scale[1] - 1mm
        place = SimpleRectangle{Float32}(
            last_x - offset[1] + 0.5mm, -offset[2] + 0.5mm, xwidth, scale[2] - 1mm
        )
        layout!(place, widget)
        _view(widget, window, camera = :fixed_pixel)
        last_x += xwidth
    end
end

function item_area(la, deleted, item_height)
    y = la.y - item_height - 2
    deleted && return IRect(la.x, la.y, la.w, 0)
    return IRect(la.x, y, la.w, item_height)
end

function widget_screen!(parentscreen = edit_screen();
        left_gap = 1.5mm, delete = Signal(false), height = icon_size()
    )
    scroll = parentscreen.inputs[:menu_scroll]
    if isempty(parentscreen.children)
        last_area = map(parentscreen.area, height, scroll) do a, ih, s
            return IRect(left_gap, a.h-ih+s, a.w-left_gap, ih)
        end
    else
        last_area = last(parentscreen.children).area
    end
    itemarea = map(GLPlot.item_area, last_area, delete, height)
    Screen(parentscreen, area = itemarea)
end

function choices(x::Vector)
    signal, vis = GLVisualize.choice_widget(x, edit_screen(), area=(48mm, 8mm))
    w = widget_screen!()
    _view(vis, w, camera = :fixed_pixel)
    signal
end

function edit_toggle()
    edit_button, no_edit_signal = toggle_button(
        imload("play.png"), rotr90(imload("play.png")), edit_screen()
    )
end

function visible_toggle()
    toggle_button(
        imload("showing.png"), imload("notshowing.png"), edit_screen()
    )
end
function delete_toggle()
    delete_button, del_signal = button(
        imload("delete.png"), edit_screen()
    )
end

function async_map2(f, init, inputs...; typ = typeof(init))
    node = Signal(typ, init, inputs)
    worker_task = @async init
        map(inputs...) do args...
            outer_task = current_task()
            hasworked = istaskdone(worker_task) #
            if istaskdone(worker_task) #
                worker_task = @async begin
                try
                    inner_worker = @async begin
                        x = f(args...)
                        push!(node, x)
                    end
                    wait(inner_worker)
                catch err
                    Base.throwto(outer_task, CapturedException(err, catch_backtrace()))
                end
            end
        end
         worker_task
    end, node
end


text:
# converters from selected RSEXPREC to Hash
# They are used to translate SEXPREC attributes into Hash

using Dates

function Base.convert(::Type{Hash}, pl::RPairList)
    res = Hash()
    for i in eachindex(pl.items)
        @inbounds setindex!(res, pl.items[i], pl.tags[i])
    end
    res
end

##############################################################################
##
## Conversion of intermediate R objects into Vector{T} and DataFrame objects
##
##############################################################################

isna(x::Int32) = x == R_NA_INT32
isna(x::Float64) = isna_float64(reinterpret(UInt64, x))
# if re or im is NA, the whole complex number is NA
isna(x::ComplexF64) = isna(real(x)) || isna(imag(x))

# convert R vector into Vector holding elements of type T
# if force_missing is true, the result is always Vector{Union{T,Missing}},
# otherwise it's Vector{T} if `rv` doesn't contain NAs
function jlvec(::Type{T}, rv::RVEC, force_missing::Bool=true) where T
    anyna = any(isna, rv.data)
    if force_missing || anyna
        res = convert(Vector{Union{T,Missing}}, rv.data)
        if anyna
            @inbounds for (i,x) in enumerate(rv.data)
                isna(x) && (res[i] = missing)
            end
        end
        return res
    else
        return convert(Vector{T}, rv.data)
    end
end

# convert R nullable vector (has an explicit NA mask) into Vector{T[?]}
function jlvec(::Type{T}, rv::RNullableVector{R}, force_missing::Bool=true) where {T, R}
    anyna = any(rv.na)
    if force_missing || anyna
        res = convert(Vector{Union{T,Missing}}, rv.data)
        anyna && @inbounds res[rv.na] .= missing
        return res
    else
        return convert(Vector{T}, rv.data)
    end
end

# convert R vector into Vector of appropriate type
function jlvec(rv::RVEC, force_missing::Bool=true)
    if inherits(rv, R_Date_Class)
        return jlvec(Dates.Date, rv, force_missing)
    elseif inherits(rv, R_POSIXct_Class)
        return jlvec(ZonedDateTime, rv, force_missing)
    else
        return jlvec(eltype(rv.data), rv, force_missing)
    end
end

# convert R logical vector (uses Int32 to store values) into Vector{Bool[?]}
function jlvec(rl::RLogicalVector, force_missing::Bool=true)
    anyna = any(isna, rl.data)
    if force_missing || anyna
        return Union{Bool,Missing}[ifelse(isna(x), missing, x != 0) for x in rl.data]
    else
        return Bool[x != 0 for x in rl.data]
    end
end

# kernel method that converts Vector{Int32} into Vector{R} replacing R_NA_INT32 with 0
# it's assumed that v fits into R
na2zero(::Type{R}, v::Vector{Int32}) where R =
    [ifelse(!isna(x), x % R, zero(R)) for x in v]

# convert to CategoricalVector{String[?]} if `ri` is a factor,
# to Dates.Date is it's a vector of dates and
# to Vector{Int32[?]} otherwise
function jlvec(ri::RIntegerVector, force_missing::Bool=true)
    if isfactor(ri)
        return jlvec(CategoricalArray, ri, force_missing)
    elseif inherits(ri, R_Date_Class)
        return jlvec(Dates.Date, ri, force_missing)
    else
        return jlvec(eltype(ri.data), ri, force_missing)
    end
end

# convert R factor to CategoricalArray
function jlvec(::Type{CategoricalArray}, ri::RVEC, force_missing::Bool=true)
    @assert isfactor(ri)

    rlevels = getattr(ri, "levels", emptystrvec)
    sz = length(rlevels)
    REFTYPE = sz <= typemax(UInt8)  ? UInt8 :
              sz <= typemax(UInt16) ? UInt16 :
              sz <= typemax(UInt32) ? UInt32 :
                                      UInt64
    refs = na2zero(REFTYPE, ri.data)
    anyna = any(iszero, refs)
    pool = CategoricalPool{String, REFTYPE}(rlevels, inherits(ri, "ordered"))
    if force_missing || anyna
        return CategoricalArray{Union{String, Missing}, 1}(refs, pool)
    else
        return CategoricalArray{String, 1}(refs, pool)
    end
end

# convert R Date to Dates.Date
function jlvec(::Type{Dates.Date}, rv::RVEC, force_missing::Bool=true)
    @assert inherits(rv, R_Date_Class)
    nas = isnan.(rv.data)
    if force_missing || any(nas)
        dates = Union{Dates.Date, Missing}[isna ? missing : rdays2date(dtfloat)
                 for (isna, dtfloat) in zip(nas, rv.data)]
    else
        dates = rdays2date.(rv.data)
    end
    return dates
end

# convert R POSIXct to ZonedDateTime
function jlvec(::Type{ZonedDateTime}, rv::RVEC, force_missing::Bool=true)
    @assert class(rv) == R_POSIXct_Class
    tz, validtz = getjuliatz(rv)
    nas = isnan.(rv.data)
    if force_missing || any(nas)
        datetimes = Union{ZonedDateTime, Missing}[isna ? missing : _unix2zdt(dtfloat, tz=tz)
                     for (isna, dtfloat) in zip(nas, rv.data)]
    else
        datetimes = _unix2zdt.(rv.data, tz=tz)
    end
    return datetimes
end

function sexp2julia(rex::RSEXPREC)
    @warn "Conversion of $(typeof(rex)) to Julia is not implemented"
    return nothing
end

function sexp2julia(rv::RVEC)
    # TODO dimnames?
    # FIXME add force_missing option to control whether always convert to Union{T, Missing}
    jv = jlvec(rv, false)
    if hasnames(rv)
        # if data has no NA, convert to simple Vector
        return DictoVec(jv, names(rv))
    else
        hasdims = hasdim(rv)
        if !hasdims && length(rv.data)==1
            # scalar
            return jv[1]
        elseif !hasdims
            # vectors
            return jv
        else
            # matrices and so on
            dims = tuple(convert(Vector{Int}, getattr(rv, "dim"))...)
            return reshape(jv, dims)
        end
    end
end

function sexp2julia(rl::RList)
    if isdataframe(rl)
        # FIXME add force_missing option to control whether always convert to Union{T, Missing}
        DataFrame(Any[isa(col, RAltRep) ? sexp2julia(col) : jlvec(col, false) for col in rl.data],
                  identifier.(names(rl)), makeunique=true)
    elseif hasnames(rl)
        DictoVec(Any[sexp2julia(item) for item in rl.data], names(rl))
    else
        # FIXME return DictoVec if forceDictoVec is on
        map(sexp2julia, rl.data)
    end
end

function sexp2julia(ar::RAltRep)
    # detect if it's a standard R 3.5+ wrapper
    if isa(ar.info, RPairList) && length(ar.info) >= 1 &&
        isa(ar.info.items[1], RSymbol) && startswith(string(ar.info.items[1]), "wrap_")
        # the first element of the AltRep state should be the wrapped one
        if isa(ar.state, RPairList) && length(ar.state) >= 1
            return sexp2julia(ar.state.items[1])
        else
            error("Unexpected state of \"$(ar.info.items[1])\" AltRep SEXP")
        end
    else
        # TODO support compact_intseq and compact_realseq AltRep
        @warn "Unsupported AltRep SEXP"
        return nothing
    end
end

function rdays2date(days::Real)
    epoch_conv = 719528 # Dates.date2epochdays(Date("1970-01-01"))
    Dates.epochdays2date(days + epoch_conv)
end

# gets R timezone from the data attribute and converts it to TimeZones.TimeZone
# see r2juliatz()
function getjuliatz(rv::RVEC, deftz=tz"UTC")
    tzattr = getattr(rv, "tzone", [""])[1]
    if tzattr == ""
        return deftz, true # R will store a blank for tzone
    else
        return r2juliatz(tzattr, deftz)
    end
end

# converts R timezone code to TimeZones.TimeZone
# returns a tuple:
#  - timezone (or `deftz` if `rtz` is not recognized as a valid time zone)
#  - boolean flag: true if `rtz` is not recognized, false otherwise
function r2juliatz(rtz::AbstractString, deftz=tz"UTC")
    valid = istimezone(rtz)
    if !valid
        @warn "Could not determine the timezone of '$(rtz)', treating as $deftz."
        return deftz, false
    else
        return TimeZone(rtz), true
    end
end

# version with user-specified tz (unix2zdt(seconds) is fixed to tz"UTZ")
_unix2zdt(seconds::Real; tz::TimeZone=tz"UTC") =
    ZonedDateTime(Dates.unix2datetime(seconds), tz, from_utc=true)


text:
using DrWatson
@quickactivate "KrotovTests"
import Downloads

DOWNLOADS = Dict(
    "https://github.com/JuliaQuantumControl/Krotov.jl/raw/data-dump/DissGateOCT%23J_T%3DJ_T_re%23iter_stop%3D3000%23method%3Dkrotov.jld2" =>
        joinpath(datadir(), "DissGateOCT#J_T=J_T_re#iter_stop=3000#method=krotov.jld2"),
    "https://github.com/JuliaQuantumControl/Krotov.jl/raw/data-dump/parametrization%23opt_result_logisticsq.jld2" =>
        joinpath(datadir(), "parametrization#opt_result_logisticsq.jld2"),
    "https://github.com/JuliaQuantumControl/Krotov.jl/raw/data-dump/parametrization%23opt_result_positive.jld2" =>
        joinpath(datadir(), "parametrization#opt_result_positive.jld2"),
    "https://github.com/JuliaQuantumControl/Krotov.jl/raw/data-dump/parametrization%23opt_result_tanh.jld2" =>
        joinpath(datadir(), "parametrization#opt_result_tanh.jld2"),
    "https://github.com/JuliaQuantumControl/Krotov.jl/raw/data-dump/parametrization%23opt_result_tanhsq.jld2" =>
        joinpath(datadir(), "parametrization#opt_result_tanhsq.jld2"),
    "https://github.com/JuliaQuantumControl/Krotov.jl/raw/data-dump/PE_OCT.jld2" =>
        joinpath(datadir(), "PE_OCT.jld2"),
    "https://github.com/JuliaQuantumControl/Krotov.jl/raw/data-dump/PE_OCT_direct.jld2" =>
        joinpath(datadir(), "PE_OCT_direct.jld2"),
)

function download_dump(url, destination; force=false, verbose=true)
    if !isfile(destination) || force
        verbose && (@info "Downloading $url => $destination")
        Downloads.download(url, destination)
    else
        verbose && (@info "$destination OK")
    end
end

@info "Download Dumps"
for (url, destination) in DOWNLOADS
    download_dump(url, destination)
end


text:
using Plots

"""
    Strategy(decisions, scores, energies, activeMons)

Struct for information about an entire match and the strategy used throughout.
This includes the decisions made, the scores at each point along the battle,
the energies of the mons (to determine fast and charged move applications), and
the active mons used by each team.
"""
mutable struct Strategy
    decisions::Array{Tuple{UInt8,UInt8}}
    scores::Array{Float64}
    hps::Array{Tuple{Tuple{UInt16, UInt16, UInt16}, Tuple{UInt16, UInt16, UInt16}}}
    activeMons::Array{Tuple{UInt8, UInt8}}
end

"""
    plot_strategy(strat, static_s)

Given the Strategy and the StaticState, plot the battle using a PvPoke-like
notation.
"""
function plot_strategy(strat::Strategy, static_s::StaticState)
    gr()
    plt1 = plot(1:length(strat.scores), strat.scores, width = 2,
        label = "possible battle scores", ylims = [0, 1],
        ylabel = "Battle Score", xlabel = "Decisions", size = (950, 400))
    plot!(plt1, 1:length(strat.scores),
        map(i -> strat.hps[i][1][1], 1:length(strat.scores)) ./
        static_s[0x01][0x01].stats.hitpoints, color = :green, label = "")
    plot!(plt1, 1:length(strat.scores),
        map(i -> strat.hps[i][1][2], 1:length(strat.scores)) ./
        static_s[0x01][0x02].stats.hitpoints, color = :green, label = "")
    plot!(plt1, 1:length(strat.scores),
        map(i -> strat.hps[i][1][3], 1:length(strat.scores)) ./
        static_s[0x01][0x03].stats.hitpoints, color = :green, label = "")
    plot!(plt1, 1:length(strat.scores),
        map(i -> strat.hps[i][2][1], 1:length(strat.scores)) ./
        static_s[0x02][0x01].stats.hitpoints, color = :purple, label = "")
    plot!(plt1, 1:length(strat.scores),
        map(i -> strat.hps[i][2][2], 1:length(strat.scores)) ./
        static_s[0x02][0x02].stats.hitpoints, color = :purple, label = "")
    plot!(plt1, 1:length(strat.scores),
        map(i -> strat.hps[i][2][3], 1:length(strat.scores)) ./
        static_s[0x02][0x03].stats.hitpoints, color = :purple, label = "")

    hline!(plt1, [0.5], label = "win/loss")
    plt2 = plot(xlims = [0, length(strat.scores)], ylims = [-3, 0],
        legend = false, size = (950, 200), axis = nothing)
    for i = 1:length(strat.scores), j = 0x01:0x02
        o_agent = j == 0x01 ? 0x02 : 0x01
        if strat.decisions[i][j] == 0x01
            scatter!(plt2, [i], [-Int64(j)], markershape = :hexagon, markersize = 12,
                alpha = 0.5, color = shieldColor)
        elseif strat.decisions[i][j] == 0x05 || strat.decisions[i][j] == 0x06
            color = colors[static_s[j][strat.activeMons[i][j]].primary_type]
            scatter!(plt2, [i - .25], [-Int64(j) + .25], markershape = :utriangle,
                alpha = 0.5, color = color)
            if 1 <= static_s[j][strat.activeMons[i][j]].secondary_type <= 18
                color = colors[static_s[j][strat.activeMons[i][j]].secondary_type]
            end
            scatter!(plt2, [i + .25], [-Int64(j) - .25], markershape = :dtriangle,
                alpha = 0.5, color = color)
        elseif strat.decisions[i][j] == 0x07 || strat.decisions[i][j] == 0x08
            color = colors[strat.decisions[i][j] == 0x07 ? static_s[j][
                strat.activeMons[i][j]].charged_move_1.move_type : static_s[j][
                    strat.activeMons[i][j]].charged_move_2.move_type]
            scatter!(plt2, [i], [-Int64(j)], markershape = :circle, markersize = 10,
                alpha = 0.5, color = color)
        elseif i > 1 && sum(strat.hps[i][o_agent]) < sum(
            strat.hps[i - 1][o_agent])
            color = colors[static_s[j][
                strat.activeMons[i][j]].fast_move.move_type]
            scatter!(plt2, [i], [-Int64(j)], markershape = :square, alpha = 0.5,
                color = color)
        end
        if i > 1
            if strat.decisions[i][j] != 0x07 && strat.decisions[i][j] != 0x08 &&
                sum(strat.hps[i][o_agent]) <
                sum(strat.hps[i - 1][o_agent])
                color = colors[static_s[j][
                    strat.activeMons[i-1][j]].fast_move.move_type]
                scatter!(plt2, [i], [-Int64(j)], markershape = :square, alpha = 0.5,
                    color = color)
            end
            if count(iszero, strat.hps[i][j]) >
                count(iszero, strat.hps[i - 1][j])
                scatter!(plt2, [i], [-Int64(j)], markershape = :xcross, alpha = 0.5,
                    markersize = 10, color = :red)
            end
        end
    end

    l = @layout [a; b{0.2h}]
    plot(plt1, plt2, layout = l, size = (950, 600))
end


text:
#!/usr/bin/julia

# Copyright (C) 2021 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Twelve, Passage Pathing."

#= The recursive solution works on a good time, therefore no need to use
 = dynamic programming, at least now.  Maybe I will return and refactor this
 = solution someday.
 =#

include("Utils.jl")

using DelimitedFiles
using .Utils

"Create adjacency matrix from `paths`."
function create_adjacency_matrix(paths :: Matrix{String})
    vertices = unique(paths)
    vert_len = length(vertices)

    vert_assoc = Dict(zip(vertices, 1:vert_len))
    path_assoc = zip(paths[:,1], paths[:,2])

    adj_mat = BitArray{2}(zeros(Bool, (vert_len, vert_len)))

    for path in path_assoc
        if "start" in path
            to = path[findfirst(path .!= "start")]
            adj_mat[vert_assoc["start"], vert_assoc[to]] = 1
        elseif "end" in path
            from = path[findfirst(path .!= "end")]
            adj_mat[vert_assoc[from], vert_assoc["end"]] = 1
        else
            from, to = path

            adj_mat[vert_assoc[from], vert_assoc[to]] = 1
            adj_mat[vert_assoc[to], vert_assoc[from]] = 1
        end
    end

    return adj_mat
end

"""Find the number of possible `paths` in the cavern system.  If `can_twice`,
allow at most two visits to a single small cave; if not, only allow one visit
per small cave."""
function solve(paths :: Matrix{String}, can_twice = false)
    vertices = unique(paths)

    vert_assoc = Dict(zip(vertices, 1:length(vertices)))
    adj_mat = create_adjacency_matrix(paths)

    small_cave = Array{Int}([])

    for v in vertices
        if all(islowercase.(collect(v))) && v ∉ ["start", "end"]
            push!(small_cave, vert_assoc[v])
        end
    end

    """Recurse through the paths, recording the small caves `visited` and the
    `start` point. For `can_twice`, vide supra."""
    function rec_path(can_twice = false, visited = [], start = vert_assoc["start"])
        if start == vert_assoc["end"]
            return 1
        elseif start in visited
            if !can_twice
                return 0
            else
                can_twice = false
            end
        end

        new_visited = visited

        if start in small_cave
            new_visited = copy(visited)
            push!(new_visited, start)
        end

        starts = findall(adj_mat[start,:])

        return sum([rec_path(can_twice, new_visited, st) for st in starts])
    end

    return rec_path(can_twice)
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = Matrix{String}(readdlm(first(ARGS), '-'))

    println(solve(data))
    println(solve(data, true))
end

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


text:
using DelimitedFiles

input = joinpath(@__DIR__, "input")
crabs = readdlm(input, ',', Int)

function align_crabs(crabs, fuel_calc)
    acc = []
    for i in minimum(crabs):maximum(crabs)
        append!(acc, sum(map(n -> fuel_calc(i, n), crabs)))
    end
    return acc
end

p1 = minimum(align_crabs(crabs, (i, n) -> abs(i - n)))
p2 = minimum(align_crabs(crabs, (i, n) -> div(abs(i - n) * (abs(i - n) + 1), 2)))

println("-----------------------------------------------------------------------")
println("the treachery of whales -- part one :: $p1")
println("the treachery of whales -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 343468)
@assert(p2 == 96086265)


text:
# This file contains the definition for the PowerPlants struct.

mutable struct PowerPlants
    unit
    node
    pmax
    pmin
    rup
    rdn
    timeup
    timedn
    mc
    varcost
    fuel
    efficiency
    emission

    function PowerPlants(pp_df::DataFrame)

        G = pp_df[!, 1]                                    # U = unit
        pp_dict = df_to_dict_with_id(pp_df)

        return new(G,
                   pp_dict[:node],
                   pp_dict[:pmax],
                   pp_dict[:pmin],
                   pp_dict[:rup],
                   pp_dict[:rdn],
                   pp_dict[:timeup],
                   pp_dict[:timedn],
                   pp_dict[:mc],
                   0,
                   0,
                   0,
                   0,
                   )
    end

    function PowerPlants(pp_df::DataFrame, fuelcost_df::DataFrame)

        G = pp_df[!, 1]                                    # U = unit
        pp_dict = df_to_dict_with_id(pp_df)

        # Marginal generation cost from fuel cost, efficiency, variable cost
        # and CO2 price

        fuelcost_dict = df_to_dict(fuelcost_df)
        mc = Dict{Int64, Array{Float64, 1}}()

        for p in keys(pp_dict[:fuel])
            f = Symbol(pp_dict[:fuel][p])
            fuelcost = fuelcost_dict[f]
            η = pp_dict[:efficiency][p]
            co_price = fuelcost_dict[:CO2]
            emission = pp_dict[:emission][p]
            vc = pp_dict[:varcost][p]

            mc[p] = fuelcost / η + co_price * emission .+ vc
        end

        return new(G,
                   pp_dict[:node],
                   pp_dict[:pmax],
                   pp_dict[:pmin],
                   pp_dict[:rup],
                   pp_dict[:rdn],
                   pp_dict[:timeup],
                   pp_dict[:timedn],
                   mc,
                   pp_dict[:varcost],
                   pp_dict[:fuel],
                   pp_dict[:efficiency],
                   pp_dict[:emission]
                   )
    end
end


text:
function setup(backend::GPUBackend, layer::GaussianReconLossLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    aux_ones = make_blob(backend, data_type, size(inputs[2])...) # for summing
    fill!(aux_ones, 1.0)
    tmp_blob = make_blob(backend, data_type, size(inputs[2])...)
    scaled_sq_errs = make_blob(backend, data_type, size(inputs[2])...)
    tmp = @compat(Dict(:aux_ones => aux_ones,
                       :tmp => tmp_blob,
                       :scaled_sq_errs => scaled_sq_errs))
    state = GaussianReconLossLayerState(layer, zero(data_type), zero(data_type), 0, tmp)
    return state
end

function forward(backend::GPUBackend, state::GaussianReconLossLayerState, inputs::Vector{Blob})
    data_type = eltype(inputs[1])
    n = get_num(inputs[1])
    np = length(inputs[1])

    μ = inputs[1]
    σ = inputs[2]
    x = inputs[3]
    aux_ones = state.tmp_blobs[:aux_ones]
    logσ² = state.tmp_blobs[:tmp]
    scaled_sq_errs = state.tmp_blobs[:scaled_sq_errs]

    copy!(logσ², σ)
    CuVec.log!(backend, logσ²)
    CuVec.mul_scal!(backend, logσ², convert(data_type, 2.))
    Σlogσ² = CuBLAS.dot(backend.cublas_ctx, data_type, np, logσ².ptr, 1, aux_ones.ptr, 1)

    copy!(scaled_sq_errs, μ)
    CuVec.sub!(backend, scaled_sq_errs, x)
    CuVec.div!(backend, scaled_sq_errs, σ)

    CuVec.pow!(backend, scaled_sq_errs, convert(data_type, 2.))

    Σscaled_sq_errs = CuBLAS.dot(backend.cublas_ctx, data_type, np,
                              scaled_sq_errs.ptr, 1, aux_ones.ptr, 1)

    state.loss = np * log(2pi) + Σlogσ² + Σscaled_sq_errs
    state.loss *= 0.5 * state.layer.weight / n

    # accumulate statistics
    state.loss_accum *= state.n_accum
    state.loss_accum += state.loss * n
    state.loss_accum /= state.n_accum + n

    state.n_accum += n
end

function backward(backend::GPUBackend, state::GaussianReconLossLayerState,
            inputs::Vector{Blob}, diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    n = get_num(inputs[1])
    np = length(inputs[1])

    μ = inputs[1]
    σ = inputs[2]
    x = inputs[3]

    tmp = state.tmp_blobs[:tmp]
    scaled_sq_errs = state.tmp_blobs[:scaled_sq_errs]

    c = convert(data_type, state.layer.weight / n)

    if isa(diffs[1], CuTensorBlob)
        copy!(diffs[1], μ)
        CuVec.sub!(backend, diffs[1], x)
        CuVec.div!(backend, diffs[1], σ)
        CuVec.div!(backend, diffs[1], σ)

        CuVec.mul_scal!(backend, diffs[1], c)
    end

    if isa(diffs[2], CuTensorBlob)
        copy!(diffs[2], σ)
        CuVec.pow!(backend, diffs[2], convert(data_type, -1.0))

        copy!(tmp, scaled_sq_errs)
        CuVec.div!(backend, tmp, σ)

        CuVec.sub!(backend, diffs[2], tmp)
        CuVec.mul_scal!(backend, diffs[2], c)
    end

    if isa(diffs[3], CuTensorBlob)
        raise("last blob should be a data blob")
    end
end


text:
function CjSjCoefficient(arg0::jdouble, arg1::jdouble)
    return CjSjCoefficient((jdouble, jdouble), arg0, arg1)
end

function get_cj(obj::CjSjCoefficient, arg0::jint)
    return jcall(obj, "getCj", jdouble, (jint,), arg0)
end

function get_dcj_dh(obj::CjSjCoefficient, arg0::jint)
    return jcall(obj, "getDcjDh", jdouble, (jint,), arg0)
end

function get_dcj_dk(obj::CjSjCoefficient, arg0::jint)
    return jcall(obj, "getDcjDk", jdouble, (jint,), arg0)
end

function get_dsj_dh(obj::CjSjCoefficient, arg0::jint)
    return jcall(obj, "getDsjDh", jdouble, (jint,), arg0)
end

function get_dsj_dk(obj::CjSjCoefficient, arg0::jint)
    return jcall(obj, "getDsjDk", jdouble, (jint,), arg0)
end

function get_sj(obj::CjSjCoefficient, arg0::jint)
    return jcall(obj, "getSj", jdouble, (jint,), arg0)
end


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

"""
    soft_string_drag!(dF, F, p, s)

Defines the soft string ordinary differential equation for the axial force in
a drillstring while tripping in/out of a well, i.e.

\$\\frac{dF_t}{ds} = -w_{bp} \\hat{t}_3 + \\mu w_c,\$

where

\$w_{bp} = A g (\\rho_p - \\rho_m),\$

and

\$w_c = \\sqrt{\\left(F_t \\kappa + w_{bp} \\hat{n}_3\\right)^2 + \\left(w_{bp} \\hat{b}_3 \\right)^2}.\$

\$\\hat{t}, \\hat{n}, \\hat{b}\$ are the tangential, normal, and bi-normal vectors
in the Frenet frame that follows the drillstring trajectory.  From the position
vector along trajectory curve defined by the B-Splines \$\\vec{r} = \\vec{r}(s)\$
these are defined as

\$\\hat{t} = \\frac{\\vec{r}^\\prime}{\\left\\Vert \\vec{r}^\\prime \\right\\Vert},\$

\$\\hat{b} = \\frac{\\vec{r}^{\\prime\\prime}}{\\left\\Vert \\vec{r}^{\\prime\\prime} \\right\\Vert},\$

\$\\hat{n} = \\hat{t} \\times \\hat{b},\$

and finally, \$\\kappa\$ is the curvature defined by

\$\\kappa = \\frac{\\left\\Vert \\vec{r}^{\\prime} \\times \\vec{r}^{\\prime\\prime} \\right\\Vert}{\\left\\Vert \\vec{r}^\\prime \\right\\Vert^3}\$.

\$\\mu\$ is the friction factor, \$\\rho_p\$ and \$\\rho_m\$ are the densities of the
drill pipe and drilling mud respectively, \$A\$ is the cross-sectional area of
the pipe, and \$g\$ is acceration due to gravity.

# Arguments
- `p::Vector{<:Number}`: Vector that defines the parameters of the problem.
  Must be in the following order:

  p = [μ, rₒ, rᵢ, g, ρₛ, ρₘ, c]

  where

  μ - friction factor
  rₒ - outer radius of drill pipe
  rᵢ - inner radius of drill pipe
  g - acceleration due to gravity
  ρₛ - density of drill pipe
  ρₘ - density of drilling mud
  c - `BSplineCurve` that defines the wellbore trajectory
"""
function soft_string_drag!(dF,F,p,s)
  Fₜ, = F
  μ, rₒ, rᵢ, g, ρₛ, ρₘ, c = p

  w_dp = π / 4 * (rₒ ^ 2 - rᵢ ^ 2)  * g * (ρₛ- ρₘ)

  curves = c(s)

  r⃗′   = curves[2,:]
  r⃗′′  = curves[3,:]

  t̂ = r⃗′ / norm(r⃗′)
  n̂ = r⃗′′ / norm(r⃗′′)
  b̂ = t̂ × n̂

  κ = norm(r⃗′ × r⃗′′) / norm(r⃗′) ^ 3

  wₛ = sqrt((Fₜ * κ + w_dp * n̂[3])^ 2 + (w_dp * b̂[3]) ^ 2)

  dF[1] =  -w_dp * t̂[3] + μ * wₛ
end

function stiff_string_drag!(dF,F,p,s)
  Fₜ, Fₙ, F_b = F
  μ, E, rₒ , rᵢ, g, ρₛ, ρₘ, c = p

  EI = E * π * (rₒ ^ 4 - rᵢ ^ 4) / 64
  w_dp = π / 4 * (rₒ ^ 2 - rᵢ ^ 2)  * g * (ρₛ- ρₘ)

  curves = c(s)

  r⃗′    = curves[2,:]
  r⃗′′   = curves[3,:]
  r⃗′′′  = curves[4,:]

  t̂ = r⃗′ / norm(r⃗′)
  b̂ = (r⃗′ × r⃗′′) /  norm(r⃗′ × r⃗′′)
  n̂ = b̂ × t̂

  κ = norm(r⃗′ × r⃗′′) / norm(r⃗′) ^ 3
  τ = r⃗′ ⋅ (r⃗′′ × r⃗′′′) / norm(r⃗′ × r⃗′′) ^ 2

  # c⃗′′′ = (r⃗′′′ * norm(r⃗′) ^ 3 /
          # ((r⃗′ ⋅ r⃗′) * (r⃗′ ⋅ r⃗′′′ + r⃗′′ ⋅ r⃗′′) - (r⃗′ ⋅ r⃗′′) ^ 2))
  # κ′ = c⃗′′′ ⋅ n̂

  # a = Fₙ+ EI * κ′
  # b = -F_b - EI * κ * τ
  a = Fₜ * κ + w_dp * n̂[3]
  b = w_dp * b̂[3]
  wₛ = sqrt(a ^ 2 + b ^ 2) # / μ / rₒ
  θ = atan(b, a)

  dF[1] =  κ * Fₙ - w_dp * t̂[3] + μ * wₛ
  dF[2] = -κ * Fₜ - w_dp * n̂[3] - wₛ * cos(θ) - F_b * τ
  dF[3] =  τ * Fₙ - w_dp * b̂[3] + wₛ * sin(θ)
end

export stiff_string_drag!, soft_string_drag!

# %%
#


text:
using HorizonSideRobots
include("BackPath.jl")
include("Functions.jl")

function try_move!(r::Robot,side::HorizonSide)::Int
    ort_side = left(side)
    n=0
    while isborder(r,side)
        if !isborder(r,ort_side)
            move!(r,ort_side)
            n += 1
        else
            break
        end
    end
    if isborder(r,side)
        movements!(r,reverse_side(ort_side),n)
        return -2
    end
    move!(r,side)
    m = 0
    while (isborder(r,reverse_side(ort_side)) && n != 0)
        move!(r,side)
        m += 1
    end
    movements!(r,reverse_side(ort_side),n)
    return m
end

function try_move_to_wall!(r::Robot, side::HorizonSide)
    while (try_move!(r, side) != -2)
    end
end

function move_n_steps_putmarkers!(r::Robot, side::HorizonSide, n::Int)
    putmarker!(r)
    while n > 0
        m = try_move!(r, side)
        if m < 0
            m = 0
        end
        n = n - 1 - m
        if n > -1
            putmarker!(r)
        end
    end
end

function paint_stair(r::Robot)
    back_path = BackPath(r)
    len_cur_step = move_with_path!(r, Ost)
    movements!(r, West, len_cur_step)
    while (!isborder(r, Nord) || !isborder(r, West)) && (len_cur_step > 0)
        move_n_steps_putmarkers!(r, Ost, len_cur_step)
        len_cur_step -= 1
        try_move_to_wall!(r, West)
        move!(r, Nord)
    end
    if isborder(r, Nord)
        move_n_steps_putmarkers!(r, Ost, len_cur_step)
    end

    to_sud_west_corner = BackPath(r)
    back!(r, back_path)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
using JLD2, LinearAlgebra, Oceananigans, Printf
record_interaction = false
include(pwd() * "/scripts/vizinanigans.jl")
include(pwd() * "/scripts/states.jl")
include(pwd() * "/scripts/compare.jl")
include(pwd() * "/scripts/zonalstates.jl")

include(pwd() * "/analysis_scripts/" * "post_analysis.jl")

files = [
    pwd() * "/Channel_1_zonal_averages.jld2",
    pwd() * "/Channel_16_zonal_averages.jld2",
    pwd() * "/Channel_24_zonal_averages.jld2",
]

new_files = [
    pwd() * "/Channel_8_zonal_averages.jld2",
    pwd() * "/Ridge_8_zonal_averages.jld2",
    pwd() * "/Relaxation_Channel_16_zonal_averages.jld2",
    pwd() * "/Abernathy_16_zonal_averages.jld2",
    pwd() * "/Abernathy_16_τ0.05_Q10_zonal_averages.jld2",
    pwd() * "/Abernathy_1_τ0.20_Q10_zonal_averages.jld2",
    pwd() * "/Abernathy_16_τ0.20_Q10_zonal_averages.jld2",
    pwd() * "/Abernathy_16_τ0.10_Q5_zonal_averages.jld2",
    pwd() * "/Abernathy_1_τ0.10_Q5_zonal_averages.jld2",

]
##
file = new_files[end] # new_files[1]
file = pwd() * "/Abernathy_16_τ0.20_Q10_zonal_averages.jld2"
zonalstatistics = jldopen(file)
tkeys = keys(zonalstatistics["timeseries"]["t"])
states, statenames, units, domain = grabzonalstates(file, ghost = 3, startind = 20)
li = 1 # bottom
mval = 0 # positive
ui = length(domain[2])-mval # top
newstates = [state[:, li:end-mval] for state in states]
xlims = (domain[1][1], domain[1][end])
ylims = (domain[2][li], domain[2][ui])
scene = visualize(newstates, statenames = statenames, xlims = xlims, ylims = ylims, units = units)

##
seconds = 20
fps = 10
frames = round(Int, fps * seconds )
if record_interaction
record(scene, pwd() * "/zonalandtemporal.mp4"; framerate = fps) do io
    for i = 1:frames
        sleep(1/fps)
        recordframe!(io)
    end
end
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
# lazy.jl

# A collection of lazy sequences that perform a computation on-the-fly whenever
# getindex is called.
#
# In general, these lazy sequences do not specialize on the type of the underlying
# sequence(s). This avoids nested type parameters, yet it means operations
# are typically slow.

# A lazy sequence can be transformed into a real one using `collect`
function collect(s::Sequence)
    @assert iscompact(s)
    r = nzrange(s)
    data = OffsetArray([s[i] for i in r], r)
    ZeroPaddingSequence(data)
end

"""
A lazy representation of the convolution of two sequences `f` and `g`. Their
convolution `h = f * g` is defined element-wise by:
`h[n] = \sum_k f[k]*g[n-k].`

The lazy convolution performs the above calculation on the fly with every
call to `getindex`.
"""
immutable Convolution{T} <: Sequence{T}
    f   ::  Sequence{T}
    g   ::  Sequence{T}
end

# * for sequences means convolution
(*)(f::Sequence, g::Sequence) = convolve(f, g)

convolve(f::Sequence, g::Sequence) = Convolution(f, g)

getindex(s::Convolution, n) = _getindex(s, n, s.f, s.g)

function _getindex(s::Convolution, n, f, g)
    z = zero(eltype(s))
    if iscompact(f)
        for k in nzrange(f)
            z += f[k] * g[n-k]
        end
    elseif iscompact(g)
        for k in nzrange(g)
            z += g[k] * f[n-k]
        end
    else
        throw(InexactError())
    end
    z
end

iscompact(s::Convolution) = _iscompact(s, s.f, s.g)
_iscompact(s::Convolution, f, g) = iscompact(f) && iscompact(g)

nzrange(s::Convolution) = _nzrange(s, s.f, s.g)

function _nzrange(s::Convolution, f, g)
    @assert _iscompact(s, f, g)
    # We have to determine the range of indices such that f[k] and g[n-k] overlap.
    a = first(nzrange(f))
    b = last(nzrange(f))
    c = first(nzrange(g))
    d = last(nzrange(g))
    # The intervals are [a,b] and [n-d,n-c].
    # Hence, the lower bound is when n-c = a.
    # The upper bound is when n-d = b
    a+c:d+b
end

ztransform(s::Convolution, z) = _ztransform(s, z, s.f, s.g)
_ztransform(s::Convolution, z, f, g) = ztransform(f, z) * ztransform(g, z)

"""
A ShiftedSequence is the lazy representation of a sequence that is shifted by
`k` units:
`g[n] = x[n-k]`.
"""
immutable ShiftedSequence{T} <: Sequence{T}
    f   ::  Sequence{T}
    k   ::  Int
end

sequence_shift(s::ShiftedSequence) = s.k

"Shift a sequence by `k` positions forward."
shift(s::Sequence, k) = ShiftedSequence(s, k)

# Don't shift a shifted sequence twice
shift(s::ShiftedSequence, k) = ShiftedSequence(s.f, s.k+k)

getindex(s::ShiftedSequence, n) = _getindex(s, n, s.f, s.k)
_getindex(s::ShiftedSequence, n, f, k) = f[n-k]

iscompact(s::ShiftedSequence) = _iscompact(s, s.f)
_iscompact(s::ShiftedSequence, f) = iscompact(f)

nzrange(s::ShiftedSequence) = _nzrange(s, s.f, s.k)

function _nzrange(s::ShiftedSequence, f, k)
    @assert _iscompact(s, f)
    r = nzrange(f)
    first(r)+k:last(r)+k
end

ztransform(s::ShiftedSequence, z) = _ztransform(s, z, s.f, s.k)
_ztransform(s::ShiftedSequence, z, f, k) = z^(-k)*ztransform(f, z)

"""
A ReversedSequence is the lazy representation of a sequence that is reversed in
time:
`g[k] = x[-k]`.
"""
immutable ReversedSequence{T} <: Sequence{T}
    f   ::  Sequence{T}
end

reverse(s::Sequence) = ReversedSequence(s)
reverse(s::ReversedSequence) = s.f

getindex(s::ReversedSequence, k) = _getindex(s, k, s.f)
_getindex(s::ReversedSequence, k, f) = f[-k]

iscompact(s::ReversedSequence) = _iscompact(s, s.f)
_iscompact(s::ReversedSequence, f) = iscompact(f)

nzrange(s::ReversedSequence) = _nzrange(s, s.f)

function _nzrange(s::ReversedSequence, f)
    @assert _iscompact(s, f)
    r = nzrange(f)
    -last(r):-first(r)
end

ztransform(s::ReversedSequence, z) = _ztransform(s, z, s.f)
_ztransform(s::ReversedSequence, z, f) = ztransform(f, 1/z)

"""
A DownsampledSequence is the lazy representation of a sequence that is
downsampled by a factor `M`. It is defined by:
`g[k] = f[M*k]`.
"""
immutable DownsampledSequence{T} <: Sequence{T}
    f   ::  Sequence{T}
    M   ::  Int
end

samplefactor(s::DownsampledSequence) = s.M

"Return the sequence downsampled by a factor `M`. By default, `M=2`."
downsample(s::Sequence, M = 2) = DownsampledSequence(s, M)

downsample(s::DownsampledSequence, M = 2) = DownsampledSequence(s.f, s.M*M)

getindex(s::DownsampledSequence, k) = _getindex(s, k, s.f, s.M)
_getindex(s::DownsampledSequence, k, f, M) = f[M*k]

iscompact(s::DownsampledSequence) = _iscompact(s, s.f)
_iscompact(s::DownsampledSequence, f) = iscompact(f)

nzrange(s::DownsampledSequence) = _nzrange(s, s.f, s.M)

function _nzrange(s::DownsampledSequence, f, M)
    @assert _iscompact(s, f)
    r = nzrange(f)
    r1 = first(r)
    r2 = last(r)
    div(r1-1,M)+1 : div(r2-1,M)+1
end

ztransform(s::DownsampledSequence, z) = _ztransform(s, z, s.f, s.M)

function _ztransform(s::DownsampledSequence, z, f, M)
    T = promote_type(eltype(s), eltype(z), eltype(1/z))
    u = zero(T)
    for m = 0:M-1
        u += 1/M * ztransform(f, exp(-im*2*T(pi)*m/M)*z^(1/M))
    end
    u
end

"""
An UpsampledSequence is the lazy representation of a sequence that is
upsampled by a factor `M`, with intermediate zeros. It is defined by:
- `g[k] = f[k/M]`, if `k` is a multiple of `M`, e.g. `k = lM`
- `g[k] = 0`, otherwise.
"""
immutable UpsampledSequence{T} <: Sequence{T}
    f   ::  Sequence{T}
    M   ::  Int
end

samplefactor(s::UpsampledSequence) = s.M

"Return the sequence upsampled by a factor `M`. By default, `M=2`."
upsample(s::Sequence, M = 2) = UpsampledSequence(s, M)

upsample(s::UpsampledSequence, M = 2) = UpsampledSequence(s.f, s.M*M)

getindex(s::UpsampledSequence, k) = _getindex(s, k, s.f, s.M)
_getindex(s::UpsampledSequence, k, f, M) = mod(k, M) == 0 ? f[div(k,M)] : zero(eltype(s))

iscompact(s::UpsampledSequence) = _iscompact(s, s.f)
_iscompact(s::UpsampledSequence, f) = iscompact(f)

nzrange(s::UpsampledSequence) = _nzrange(s, s.f, s.M)

function _nzrange(s::UpsampledSequence, f, M)
    @assert _iscompact(s, f)
    r = nzrange(f)
    r1 = first(r)
    r2 = last(r)
    r1*M : r2*M
end

ztransform(s::UpsampledSequence, z) = _ztransform(s, z, s.f, s.M)
_ztransform(s::UpsampledSequence, z, f, M) = ztransform(f, z^M)

"""
A ModulatedSequence is the lazy representation of a sequence that is modulated:
`g[k] = (-1)^k x[k]`.
"""
immutable ModulatedSequence{T} <: Sequence{T}
    f   ::  Sequence{T}
end

modulate(s::Sequence) = ModulatedSequence(s)
modulate(s::ModulatedSequence) = s.f

getindex(s::ModulatedSequence, k) = _getindex(s, k, s.f)
# We avoid writing (-1)^k, since this throws a DomainError if k is negative
_getindex(s::ModulatedSequence, k, f) = (-1)^isodd(k) * f[k]

iscompact(s::ModulatedSequence) = iscompact(s.f)

nzrange(s::ModulatedSequence) = nzrange(s.f)

ztransform(s::ModulatedSequence, z) = ztransform(s.f, -z)


text:
module SubVectorsTests

using Gridap.Arrays

a = collect(1:10)
pini = 1
pend = 5
b = SubVector(a,pini,pend)
test_array(b,collect(pini:pend))

pini = 7
pend = 9
b = SubVector(a,pini,pend)
test_array(b,collect(pini:pend))

end # module


text:
#! /usr/bin/julia
#
#	csv_update.jl
#
#					Jul/05/2018
# --------------------------------------------------------------------
include("text_manipulate.jl")
#
# --------------------------------------------------------------------
println("*** 開始 ***")
fname_csv = ARGS[1]
key_in = ARGS[2]
population_in = ARGS[3]
println(fname_csv)
println(key_in,'\t',population_in)
#
dict_aa = csv_to_dict_proc(fname_csv)
#
if ( haskey(dict_aa,key_in))
	dict_aa = dict_update_proc(dict_aa,key_in,population_in)
	dict_display_proc(dict_aa)
	dict_to_csv_proc(dict_aa,fname_csv)
end
#
println("*** 終了 ***")
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense