text:
# In this example, we use EKI to tune the closure parameters of a HydrostaticFreeSurfaceModel
# with a CATKEVerticalDiffusivity closure in order to align the predictions of the model
# to those of a high-resolution LES data generated in LESbrary.jl. Here `predictions` refers to the
# 1-D profiles of temperature, velocity, and turbulent kinetic energy horizontally averaged over a
# 3-D physical domain.

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "../.."))

using Oceananigans
using LinearAlgebra, Distributions, JLD2, DataDeps, Random
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity
using OceanBoundaryLayerParameterEstimation
using ParameterEstimocean
using ParameterEstimocean.Parameters: closure_with_parameters
using ParameterEstimocean.PseudoSteppingSchemes
using ParameterEstimocean.EnsembleKalmanInversions: eki_objective

Random.seed!(1234)

Nz = 32
Nensemble = 128
architecture = GPU()
Δt = 5minutes
prior_type = "scaled_logit_normal"
# prior_type = "normal"
description = "Calibrating to days 1-3 of 4-day suite."

directory = "calibrate_catke_to_lesbrary_4day_5minute_take2/"
isdir(directory) || mkpath(directory)

path = joinpath(directory, "calibration_setup.txt")
o = open_output_file(path)
write(o, "$description \n Δt: $Δt \n Nz: $Nz \n Nensemble: $Nensemble \n Prior type: $prior_type \n")

#####
##### Set up ensemble model
#####

begin
    field_names = (:b, :u, :v, :e)
    fields_by_case = Dict("free_convection" => (:b, :e),
                          "weak_wind_strong_cooling" => (:b, :u, :v, :e),
                          "strong_wind_weak_cooling" => (:b, :u, :v, :e),
                          "strong_wind" => (:b, :u, :v, :e),
                          "strong_wind_no_rotation" => (:b, :u, :e)
                         )

    parameter_set = CATKEParametersRiDependent

    parameter_names = (:CᵂwΔ,  :Cᵂu★, :Cᴰ,
                    :Cˢc,   :Cˢu,  :Cˢe,
                    :Cᵇc,   :Cᵇu,  :Cᵇe,
                    :Cᴷc⁻,  :Cᴷu⁻, :Cᴷe⁻,
                    :Cᴷcʳ,  :Cᴷuʳ, :Cᴷeʳ,
                    :CᴷRiᶜ, :CᴷRiʷ)

    parameter_set = ParameterSet{CATKEVerticalDiffusivity}(Set(parameter_names),
                                nullify = Set([:Cᴬu, :Cᴬc, :Cᴬe]))

    transformation = (b = Transformation(normalization=ZScore()),
                    u = Transformation(normalization=ZScore()),
                    v = Transformation(normalization=ZScore()),
                    e = Transformation(normalization=RescaledZScore(1e-1)))

    closure = closure_with_parameters(CATKEVerticalDiffusivity(Float64;), parameter_set.settings)
end

#####
##### Build free parameters
#####

function build_prior(name)
    b = bounds(name, parameter_set)
    prior_type == "scaled_logit_normal" && return ScaledLogitNormal(bounds=b)
    prior_type == "normal" && return Normal(mean(b), (b[2]-b[1])/3)
end

free_parameters = FreeParameters(named_tuple_map(parameter_set.names, build_prior))

#####
##### Build the Inverse Problem
#####

output_map = ConcatenatedOutputMap()

function inverse_problem(path_fn, Nensemble, times)
    observations = SyntheticObservationsBatch(path_fn, times, Nz; architecture, transformation, field_names, fields_by_case)
    simulation = lesbrary_ensemble_simulation(observations; Nensemble, architecture, closure, Δt)
    ip = InverseProblem(observations, simulation, free_parameters; output_map)
    return ip
end

training_times = [1.0day, 1.75days, 2.5days, 3.25days, 4.0days]
validation_times = [0.5days, 1.0days, 1.5days, 2.0days]
testing_times = [1.0days, 3.0days, 6.0days]

training = inverse_problem(four_day_suite_path_2m, Nensemble, training_times)
validation = inverse_problem(two_day_suite_path_2m, Nensemble, validation_times)
testing = inverse_problem(six_day_suite_path_2m, Nensemble, testing_times)

write(o, "Training observations: $(summary(training.observations)) \n")
write(o, "Validation observations: $(summary(validation.observations)) \n")
write(o, "Testing observations: $(summary(testing.observations)) \n")

write(o, "Training inverse problem: $(summary(training)) \n")
write(o, "Validation inverse problem: $(summary(validation)) \n")
write(o, "Testing inverse problem: $(summary(testing)) \n")

# y = observation_map(training);
# θ = named_tuple_map(parameter_set.names, name -> default(name, parameter_set))
# G = forward_map(training, [θ])
# zc = [mapslices(norm, G .- y, dims = 1)...]

###
### Calibrate
###

iterations = 10

function estimate_noise_covariance(data_path_fns, times)
    obsns_various_resolutions = [SyntheticObservationsBatch(dp, times, Nz; architecture, transformation, field_names, fields_by_case) for dp in data_path_fns]
    Nobs = Nz * (length(times) - 1) * sum([length(obs.forward_map_names) for obs in obsns_various_resolutions[1]])
    noise_covariance = estimate_η_covariance(output_map, obsns_various_resolutions)
    noise_covariance = noise_covariance + I(Nobs) * mean(noise_covariance)/10
    return noise_covariance
end

noise_covariance = estimate_noise_covariance([four_day_suite_path_1m, four_day_suite_path_2m, four_day_suite_path_4m], training_times)

resampler = Resampler(acceptable_failure_fraction=0.2, only_failed_particles=true)

pseudo_stepping = Iglesias2021()
eki = EnsembleKalmanInversion(training; noise_covariance, pseudo_stepping, resampler, tikhonov = true)
final_params = iterate!(eki; iterations = 10, show_progress=false, pseudo_stepping)
visualize!(training, final_params;
    field_names = [:u, :v, :b, :e],
    directory,
    filename = "realizations_training_iglesias2021.png"
)
visualize!(validation, final_params;
    field_names = [:u, :v, :b, :e],
    directory,
    filename = "realizations_validation_iglesias2021.png"
)
visualize!(testing, final_params;
    field_names = [:u, :v, :b, :e],
    directory,
    filename = "realizations_testing_iglesias2021.png"
)

write(o, "Final ensemble mean: $(final_params) \n")
close(o)

###
### Summary Plots
###

plot_parameter_convergence!(eki, directory)
plot_error_convergence!(eki, directory)
plot_pairwise_ensembles!(eki, directory)

###
### Hyperparameter optimization
###

# include("hyperparameter_optimization.jl")

###
### CES
###

include("emulate_sample_forward_map.jl")

###
### Sensitivity analysis
###

include("sensitivity_analysis.jl")


text:
#=

IMPORTANT NOTE:

Unfortuantely, Julia does not include any structure within abstract types. This can pose a problem,
because it's not possible to immediately determine what requirements abstract types must meet.

Therefore, any declared abstract type must contain a comment describing those requirements, either
in the docstring or in a comment within the abstract type declaration. Of course, they are not
strictly binding, but make sure that you know why you're breaking the rules if you choose to do so.

=#

"""
    AbstractLattice{D}

Supertype for all lattices of dimension D.

Members consist of a set of primitive lattice vectors and a set of conventional lattice vectors,
each stored in `BasisVectors{D}`.

An `AbstractLattice` must implement the following checks:
  * That the basis vectors are linearly independent and form a right-handed coordinate system.
  * That the set of lattice vectors spanning the largest volume can be expressed in terms of
    integer linear combinations of the primitive basis vectors.
"""
abstract type AbstractLattice{D}
end

"""
    AbstractCrystal{D}

A crystal structure in `D` dimensions, containing information about the lattice, space group, and
atoms contained within the crystal.
"""
abstract type AbstractCrystal{D}
end

"""
    AbstractCrystalData{D}

Supertype for all datasets that may be associated with a crystal.
"""
abstract type AbstractCrystalData{D}
end

"""
    AbstractRealSpaceData{D}

Supertype for crystal data that is defined in real space. This includes atomic coordinates.
"""
abstract type AbstractRealSpaceData{D} <: AbstractCrystalData{D}
end

"""
    AbstractReciprocalSpaceData{D}

Supertype for crystal data that is defined in reciprocal space.
"""
abstract type AbstractReciprocalSpaceData{D} <: AbstractCrystalData{D}
end

"""
    AbstractDensityOfStates

Supertype for all density of states data.
"""
abstract type AbstractDensityOfStates
end

"""
    AbstractKPoints{D}

Supertype for sets of k-points.
"""
abstract type AbstractKPoints{D} <: AbstractReciprocalSpaceData{D}
end

"""
    AbstractPotential

Supertype for atomic potentials.
"""
abstract type AbstractPotential
end

"""
    AbstractPseudopotential

Supertype for atomic pseudopotentials.
"""
abstract type AbstractPseudopotential <: AbstractPotential
end

ext: jl

max_stars_repo_name: brainandforce/Xtal

max_issues_repo_name: brainandforce/Xtal

max_forks_repo_name: brainandforce/Xtal


text:
using GZip

## window_bed_file
#binValues = GZip.open(window_bed_file, "data/K562_tss.bed.gz")

## streaming_cor from a file
f = jldopen("data/small.ChromNet.jld")
data = read(f, "data")
C = streaming_cor([f], chunkSize=3, quiet=true)
@test all(abs(C .- cor(data, vardim=2)) .< 1e-10)

## streaming_cov from an in memory bit array
# @test all(abs(streaming_cov(M, chunkSize=3, quiet=true) .- cov(M')) .< 1e-10)

# ## conditional_cov
# C = cov(M')
# smallC = conditional_cov(C, 1:10, 11:P, 0.0)
# IC = inv(C)
# @test all(abs(smallC .- inv(IC[11:end,11:end])) .< 1e-10)

## cov2cor!
C = cov(data, vardim=2)
@test all(abs(cov2cor!(C) .- cor(data, vardim=2)) .< 1e-10)


text:
using ACE1pack, JuLIP, Test

@info("Test generating solver params")
lsqr_damp = 1e-3
solver = solver_params(type = :lsqr, lsqr_damp = lsqr_damp)
solver = ACE1pack.generate_solver(solver)
@test solver["solver"] == :lsqr
@test solver["lsqr_damp"] == lsqr_damp


text:
using Neo4j
using Mongoc
using UUIDs

#Types
export CloudGraphConfiguration, CloudGraph, CloudVertex, CloudEdge
export BigData, BigDataElement
export BigDataRawType

_mongoDefaultDb = "CloudGraphs"
_mongoDefaultCollection = "bindata"

# Type aliases
BigDataRawType = Union{Vector{UInt8}, Dict{String, Any}, String, Dict{Any, Any}}

mutable struct BigDataElement
    sourceName::String
    sourceId::String
    sourceParams::Dict{String, Any}
    id::String
    description::String
    data::BigDataRawType
    mimeType::String
    neoNodeId::Int
    lastSavedTimestamp::String #UTC DateTime.
    BigDataElement(id::String, desc::String, data::BigDataRawType, neoNodeId::Int; sourceName::String="Mongo", sourceId::String=string(uuid4()), sourceParams::Dict{String, Any}=Dict{String, Any}(), mimeType::String="application/octet-stream", lastSavedTimestamp::String=string(now(Dates.UTC))) = begin
        return new(sourceName, sourceId, sourceParams, id, desc, data, mimeType, neoNodeId, lastSavedTimestamp)
    end
    BigDataElement(dd::Dict{T,Any}, version::String) where {T <: String} = begin
        if(version == "1")
            return new("Mongo", dd["mongoKey"], Dict{String, Any}(), dd["mongoKey"], dd["description"], dd["data"], "application/octet-stream", dd["neoNodeId"], dd["lastSavedTimestamp"])
        elseif(version == "2")
            return new(dd["sourceName"], dd["sourceId"], dd["sourceParams"], dd["id"], dd["description"], dd["data"], dd["mimeType"], dd["neoNodeId"], dd["lastSavedTimestamp"])
        else
            error("BigDataElement version '$version' is not supported.")
        end
    end
end

mutable struct BigData
  isRetrieved::Bool
  isAvailable::Bool
  isExistingOnServer::Bool
  lastSavedTimestamp::String #UTC DateTime.
  version::String
  dataElements::Vector{BigDataElement}
  # This is just for local use, and is not saved directly into the graph.
  BigData() = new(false, false, false, "[N/A]", "2", Vector{BigDataElement}())
  BigData(isRetrieved::Bool, isAvailable::Bool, isExistingOnServer::Bool, lastSavedTimestamp::String, version::String, data::Vector{BigDataElement}) = new(isRetrieved, isAvailable, isExistingOnServer, lastSavedTimestamp, version, data)
  BigData(jsonStr::String) = begin
      dd = JSON.parse(jsonStr)
      bDE = BigDataElement[]
      for (k,v) in dd["dataElements"]
         push!(bDE, BigDataElement(v[1],Vector{UInt8}(),v[2]) )
      end
      version = haskey(dd, "version") ? dd["version"] : "1"
      new(dd["isRetrieved"],dd["isAvailable"],dd["isExistingOnServer"],dd["lastSavedTimestamp"], version, bDE)
    end
end

mutable struct CloudVertex
  packed::Any
  properties::Dict{AbstractString, Any} # UTF8String
  bigData::BigData
  neo4jNodeId::Int
  neo4jNode::Union{Nothing, Neo4j.Node}
  labels::Vector{AbstractString}
  isValidNeoNodeId::Bool
  exVertexId::Int
  isValidExVertex::Bool
  CloudVertex() = new(Union, Dict{AbstractString, Any}(), BigData(), -1, nothing, Vector{AbstractString}(), false, -1, false)
  CloudVertex(packed, properties, bigData::BigData, neo4jNodeId, neo4jNode, isValidNeoNodeId, exVertexId, isValidExVertex; labels::Vector{T}=Vector{String}()) where {T <: AbstractString} = new(packed, properties, bigData, neo4jNodeId, neo4jNode, labels, isValidNeoNodeId, exVertexId, isValidExVertex)
  CloudVertex(packed, properties, bigData::T, neo4jNodeId, neo4jNode, isValidNeoNodeId, exVertexId, isValidExVertex; labels::Vector{T}=Vector{String}()) where {T <: AbstractString} = new(packed, properties, BigData(bigData), neo4jNodeId, neo4jNode, labels, isValidNeoNodeId, exVertexId, isValidExVertex)
end

# A single configuration type for a CloudGraph instance.
struct CloudGraphConfiguration
  neo4jHost::String
  neo4jPort::Int
  neo4jUsername::String
  neo4jPassword::String
  mongoHost::String
  mongoPort::Int
  mongoIsUsingCredentials::Bool
  mongoUsername::String
  mongoPassword::String
end

mutable struct Neo4jInstance
  connection::Neo4j.Connection
  graph::Neo4j.Graph
end

mutable struct MongoDbInstance
  client::Mongoc.Client
  cgBindataCollection #::MongoCollection
end

mutable struct PackedType
  originalType::Type
  packingType::Type
  encodingFunction::Union{Function, Nothing}
  decodingFunction::Union{Function, Nothing}
end

# A CloudGraph instance
mutable struct CloudGraph
  configuration::CloudGraphConfiguration
  neo4j::Neo4jInstance
  mongo::MongoDbInstance
  encodePackedType::Function
  getpackedtype::Function
  decodePackedType::Function
  CloudGraph(configuration, neo4j, mongo, ept, gpt, dpt) = new(configuration, neo4j, mongo, ept, gpt, dpt)
end

mutable struct CloudEdge
  neo4jEdgeId::Int
  neo4jEdge::Union{Nothing,Neo4j.Relationship}
  edgeType::AbstractString #UTF8String
  neo4jSourceVertexId::Int
  SourceVertex::Union{Nothing,CloudGraphs.CloudVertex}  #neo4jSourceVertex::Union{Nothing,Neo4j.Node}
  neo4jDestVertexId::Int
  DestVertex::Union{Nothing,CloudGraphs.CloudVertex}  #neo4jDestVertex::Union{Nothing,Neo4j.Node}
  properties::Dict{AbstractString, Any} # UTF8String
  CloudEdge() = new(-1, nothing, "", -1, nothing, -1, nothing, Dict{AbstractString, Any}())
  # UTF8String
  CloudEdge(vertexSrc::CloudVertex, vertexDest::CloudVertex, edgeType::T; props::Dict{T, Any}=Dict{T, Any}()) where {T <: AbstractString} = new(
    -1, nothing, string(edgeType), # utf8(edgeType)
    vertexSrc.neo4jNodeId,
    vertexSrc, #.neo4jNode,
    vertexDest.neo4jNodeId,
    vertexDest, #.neo4jNode,
    props)
end


text:
module NewtonCotes

export NewtonCotes11pt

function NewtonCotes11pt(f, h)
    ## A function for integration with an 11pt Newton-Cotes rule, 12th order
    # Inputs:
    #   f is the function points to be integrated. length(f)%10 = 1 is required
    #   h is the desired step size
    #
    # Outputs:
    #   I is the scalar found by numerically integrating f with step size h
    #

    if mod(length(f), 10) != 1
        return -1
    end

    n = length(f)
    I = 0

    for i = 1:10:(n-1)
        I += NewtonCotes11ptInterval(f, h, i)
    end

    return I
end

function NewtonCotes11ptInterval(f, h, i)
    # A function for implementing an 11pt Newton-Cotes Rule.
    # The rule is derived from the Lagrange Interpolating Polynomial.
    # Source: Magalhaes & Magalhaes, DOI: https://doi.org/10.3844/jmssp.2010.193.204
    # Inputs:
    #   f is the vector of solution points
    #   h is the step size
    #   i is the lower index of integration. i + 10 will be the upper.
    #
    # Outputs:
    #   X is the result of the integration rule applied to a single
    #   interval of appropriate length
    #

    X = (5/299376)*h*(16067*(f[i] + f[i+10]) + 106300*(f[i+1]+f[i+9]) - 48525*(f[i+2]+f[i+8]) + 272400*(f[i+3]+f[i+7]) - 260550*(f[i+4]+f[i+6]) + 427368*f[i+5])

end

end #module

max_stars_repo_name: jehicken/high-order-and-chaos

max_issues_repo_name: jehicken/high-order-and-chaos

max_forks_repo_name: jehicken/high-order-and-chaos


text:
__precompile__()
module DownloadFunc

#------------------------------------------------------------------#
#For the time being, we need remove_response function from obspy
#This will be replaced by SeisIO modules in the near future.
#Please activate obspy enviroment before launching Julia.
# include("remove_response_obspy.jl")
# using .Remove_response_obspy
#------------------------------------------------------------------#

using SeisIO, Dates

export seisdownload_NOISE, seisdownload_EARTHQUAKE

"""
    seisdownload_NOISE(startid, InputDict::Dict)

Download seismic data, removing instrumental response and saving into JLD2 file.

# Arguments
- `startid`         : start time id in starttimelist
- `InputDict::Dict` : dictionary which contains request information
"""
function seisdownload_NOISE(startid, InputDict::Dict; testdownload::Bool=false)

    #stationlist
    stationlist     = InputDict["stationinfo"]["stationlist"]
    method      	= InputDict["stationinfo"]["stationmethod"]
    src             = InputDict["stationinfo"]["stationsrc"]
    starttime       = InputDict["starttime"]
    endtime         = InputDict["endtime"]
    DL_time_unit    = InputDict["DL_time_unit"]

	#SeisIO getdata option
	if !haskey(InputDict, "get_data_opt")
		InputDict["get_data_opt"] = [true, true, true, true, true]#: [unscale, demean, detrend, taper, ungap]
	end

	if !haskey(InputDict, "savesamplefreq")
		InputDict["savesamplefreq"] = false # default margin: 5 minutes
	end

	if !haskey(InputDict, "download_margin")
		InputDict["download_margin"] = 5 * 60 # default margin: 5 minutes
	end

    #make stlist at all processors

    starttimelist = InputDict["starttimelist"]

    #show progress
    if starttimelist[startid][end-8:end] == "T00:00:00" && !testdownload
        println("start downloading $(starttimelist[startid])")
	elseif testdownload
		print(".")
    end

	dlerror = []

    for i = 1:length(stationlist)
        #---download data---#
        requeststr = stationlist[i]

		#println(starttimelist[startid])

		# including download margin
		starttime = string(DateTime(starttimelist[startid]) - Second(InputDict["download_margin"]))
		dltime = DL_time_unit + 2 * InputDict["download_margin"]

		#println(starttime)
		#println(dltime)

		if InputDict["IsLocationBox"]

	        ex = :(get_data($(method[i]), $(requeststr), s=$(starttime), t=$(dltime), reg=$(InputDict["reg"]),
			 v=$(0), src=$(src[i]), xf=$("$requeststr.$startid.xml"), unscale=$(InputDict["get_data_opt"][1]),
			  demean=$(InputDict["get_data_opt"][2]), detrend=$(InputDict["get_data_opt"][3]),taper=$(InputDict["get_data_opt"][4]),
			  ungap=$(InputDict["get_data_opt"][5]), rr=$(InputDict["IsResponseRemove"])))

	        t_dl = @elapsed Stemp = check_and_get_data(ex, requeststr)
		else

			ex = :(get_data($(method[i]), $(requeststr), s=$(starttime), t=$(dltime),
			 v=$(0), src=$(src[i]), xf=$("$requeststr.$startid.xml"),unscale=$(InputDict["get_data_opt"][1]),
			  demean=$(InputDict["get_data_opt"][2]), detrend=$(InputDict["get_data_opt"][3]),taper=$(InputDict["get_data_opt"][4]),
			  ungap=$(InputDict["get_data_opt"][5]), rr=$(InputDict["IsResponseRemove"])))

			t_dl = @elapsed Stemp = check_and_get_data(ex, requeststr)
		end

		# Check maximum memory allocation
		if sizeof(Stemp)/1024/1024/1024 > InputDict["MAX_MEM_PER_CPU"]
			@warn "maximam allocation of memory per cpu exceeds predescribed MAX_MEM_PER_CPU.
			This may cause transient memory leak, so please track the memory usage." AllocatedMemory_GB=sizeof(Stemp)/1024/1024/1024
		end

		Isdataflag = false
		# manipulate download_margin
		manipulate_tmatrix!(Stemp, starttime, InputDict)

		for j = 1:Stemp.n
			if Stemp.misc[j]["dlerror"] == 0
				Isdataflag = true
				# downsample
				if InputDict["savesamplefreq"] isa Number
					if Stemp.fs[j] > InputDict["savesamplefreq"]
						SeisIO.resample!(Stemp, chans=j, fs=float(InputDict["savesamplefreq"]))
					end
				end
			end
		end

		# if some of SeisChannels in Stemp have a data, save temp file
		if Isdataflag
			ymd = split(starttimelist[startid], r"[A-Z]")
			(y, m, d) = split(ymd[1], "-")
			j = md2j(y, m, d)
			fname_out = join([String(y),
						string(j),
						replace(split(starttimelist[startid], 'T')[2], ':' => '.'),requeststr,
						"FDSNWS",
						src[i],"dat"],
						'.')

			# save as intermediate binary file
			t_write = @elapsed wseis(InputDict["tmppath"]*"/"*fname_out, Stemp)
		end

		if InputDict["IsXMLfileRemoved"] && ispath("$requeststr.$startid.xml")
			rm("$requeststr.$startid.xml")
		end

		push!(dlerror, !Isdataflag)

		#print("[dltime, wtime, fraction of writing]: ")
		#println([t_dl, t_write, t_write/(t_dl+t_write)])
    end

    return dlerror
end

"""
    seisdownload_EARTHQUAKE(startid, InputDict::Dict)

Download seismic data, removing instrumental response and saving into JLD2 file.

# Arguments
- `startid`         : start time id in starttimelist
- `InputDict::Dict` : dictionary which contains request information
"""
function seisdownload_EARTHQUAKE(startid, InputDict::Dict)

	method		    = InputDict["method"]
	event		    = InputDict["event"]
	reg			    = InputDict["reg"]
    pre_filt        = InputDict["pre_filt"]

    #show progress
    if mod(startid, round(0.1*length(event))+1) == 0
        println("start downloading event number: $startid")
    end
    S = SeisData()

    #---download data---#
	for j = 1:length(event[startid]["pickphase"])
		net = event[startid]["pickphase"][j]["net"]
		sta = event[startid]["pickphase"][j]["sta"]
		loc = event[startid]["pickphase"][j]["loc"]
		cha = event[startid]["pickphase"][j]["cha"]
		src = event[startid]["pickphase"][j]["src"]

		starttime = string(event[startid]["pickphase"][j]["starttime"])
		endtime = string(event[startid]["pickphase"][j]["endtime"])

		# make multiple request str

    	requeststr = join([net,sta,loc,cha], ".")

		if InputDict["IsLocationBox"]
			# request with lat-lon box
    		#argv = [method, requeststr, starttime, endtime, InputDict["reg"], 0, src, false, "$requeststr.$startid.xml"]
		    ex = :(get_data($(method), $(requeststr), s=$(starttime), t=$(endtime), reg=$(InputDict["reg"]), v=$(0), src=$(src), xf=$("$requeststr.$startid.xml")))
		    Stemp = check_and_get_data(ex, requeststr)
		else
			# request with lat-lon box
    		#argv = [method, requeststr, starttime, endtime, 0, src, false, "$requeststr.$startid.xml"]
			ex = :(get_data($(method), $(requeststr), s=$(starttime), t=$(endtime), v=$(0), src=$(src), xf=$("$requeststr.$startid.xml")))
		    Stemp = check_and_get_data(ex, requeststr)
		end

	    if Stemp.misc[1]["dlerror"] == 0 && InputDict["IsResponseRemove"]
	        #Remove_response_obspy.remove_response_obspy!(Stemp, "$requeststr.$startid.xml", pre_filt=pre_filt, zeropadlen = float(30*60), output="VEL")
			if InputDict["IsXMLfileRemoved"]
				rm("$requeststr.$startid.xml")
			else
				mkpath("./stationxml")
				mv("$requeststr.$startid.xml", "./stationxml/$requeststr.$startid.xml", force=true)
			end
		else
			if InputDict["IsXMLfileRemoved"]
				rm("$requeststr.$startid.xml")
			else
				mkpath("./stationxml")
				mv("$requeststr.$startid.xml", "./stationxml/$requeststr.$startid.xml", force=true)
			end
	    end

		#fill gap with zero
		SeisIO.ungap!(Stemp, m=true)
		replace!(Stemp.x, NaN=>0)

		append!(S, Stemp)
	end

    return S
end

"""
    check_and_get_data(ex::Expr, requeststr::String)

Download seismic data, removing instrumental response and saving into JLD2 file.

# Arguments
- `ex::Expr`        : expression of get data includin all request information

# Output
- `S::SeisData`     : downloaded SeisData
- `requeststr::String`     : request channel (e.g. "BP.LCCB..BP1")
"""
function check_and_get_data(ex::Expr, requeststr::String)
	try
		#comment out below if you want to print contents of get_data()
		#println(ex)
		S = eval(ex);

		for j = 1:S.n
			if !isempty(S.x[j])
				#download succeeded
				S.misc[j]["dlerror"] = 0
			else
				S.misc[j]["dlerror"] = 1
			end
		end

		return S

	catch y
		#println(y)
		S = SeisData(1)
		S.misc[1]["dlerror"] = 1
		S.id[1] = requeststr
		note!(S, 1, "station is not available for this request.")
		return S
	end
end

"""
    manipulate_tmatrix!(S::SeisData, InputDict::Dict{String,Any})

    manipulate time matrix to remove download margin

"""

function manipulate_tmatrix!(S::SeisData, starttime::String, InputDict::Dict{String,Any})

    #print("before")
    #println(S)
    for i = 1:S.n

		if S.misc[i]["dlerror"] == 1
			continue;
		end

        download_margin = InputDict["download_margin"]
        DL_time_unit    = InputDict["DL_time_unit"]
        #requeststr = "NC.PLO..EHZ" #NC.PCC..EHZ
        requeststr = S.id[i] #NC.PCC..EHZ

		# NOTE: 2020/2/8
		# We decided to use data ungap to avoid data inconsitency
		# due to data gap and time shift.
		if size(S.t[i], 1) > 2
			warning("ungap! is not applied yet to SeisData, which may cause
			data inconsistency. So applying ungap.")
			ungap!(S[i])
		end

        tvec = collect(0:S.t[i][end,1]-1) ./ S.fs[i]
        tlen = trunc(Int, DL_time_unit * S.fs[i])

        si = findfirst(x -> tvec[x] >= download_margin, 1:length(tvec))

		#println([si, download_margin * S.fs[i]])
		#println([string(u2d(S.t[i][1,2] * 1e-6))[1:19],starttime])
		# check if data is within request time window AND start time is equal
		# at the order of second to what is requested
		# rounding subsecond error in downloading
		# NOTE: thie is just used to check the consistency between requested
		# time and truncated time;
		# SeisNoise.phase_shift! is applied when seisxcorr is performed.

		tsync = round(Int, S.t[i][1,2] * 1e-6) * 1e6

		if isnothing(si)
			#downloaded data is only within margin; nodata in the requested time window
			S.misc[i]["dlerror"] = 1
            S.x[i] = zeros(0)

        elseif si < download_margin * S.fs[i] || string(u2d(tsync * 1e-6))[1:19] != starttime
            #println("data missing or starttime not match.")
            S.misc[i]["dlerror"] = 1
            S.x[i] = zeros(0)
        else
            #println("manipulate")
			#===
			Shift time index of si:
			|--------|----------------------------------------|--------|
			  margin              requested data (tlen)         margin
			        [si]                                 [si]+tlen-1

			===#
            ei = trunc(Int, min(si+tlen-1, S.t[i][end,1]))
            x_shifted = zeros(tlen)
            copyto!(x_shifted, S.x[i][si:ei])

			# # manipulate time matrix
			# t_shifted = deepcopy(hcat(S.t[i][:,1] .- (si-1), S.t[i][:,2]))
			# # remove gap outside of time window
			# ri = findall(x -> t_shifted[x, 1] < 1 || t_shifted[x, 1] >  tlen, 1:size(t_shifted, 1))
			# t_shifted = t_shifted[setdiff(1:end, ri), :]
			t_init = [1 trunc(Int, S.t[i][1,2] + float(download_margin)*1e6)]
			t_last = [tlen 0]
			t_shifted = vcat(t_init, t_last)
			S.x[i] = x_shifted
			S.t[i] = t_shifted

        end
    end
    #print("after")
    #println(S)
end

end


text:
function update(
    registry_path = joinpath(homedir(), ".julia", "registries"),
    registry_name = "General"
)
    # gen_stdlib()
    generate(joinpath(registry_path, registry_name))
end


text:
# ---
# title: 1590. Make Sum Divisible by P
# id: problem1590
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Hash Table, Math, Binary Search
# link: <https://leetcode.com/problems/make-sum-divisible-by-p/description/>
# hidden: true
# ---
#
# Given an array of positive integers `nums`, remove the **smallest** subarray
# (possibly **empty** ) such that the **sum** of the remaining elements is
# divisible by `p`. It is **not** allowed to remove the whole array.
#
# Return _the length of the smallest subarray that you need to remove, or_`-1`
# _if it 's impossible_.
#
# A **subarray** is defined as a contiguous block of elements in the array.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [3,1,4,2], p = 6
#     Output: 1
#     Explanation: The sum of the elements in nums is 10, which is not divisible by 6. We can remove the subarray [4], and the sum of the remaining elements is 6, which is divisible by 6.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [6,3,5,2], p = 9
#     Output: 2
#     Explanation: We cannot remove a single element to get a sum divisible by 9. The best way is to remove the subarray [5,2], leaving us with [6,3] with sum 9.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2,3], p = 3
#     Output: 0
#     Explanation: Here the sum is 6. which is already divisible by 3. Thus we do not need to remove anything.
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,2,3], p = 7
#     Output: -1
#     Explanation: There is no way to remove a subarray in order to get a sum divisible by 7.
#
#
# **Example 5:**
#
#
#
#     Input: nums = [1000000000,1000000000,1000000000], p = 3
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 105`
#   * `1 <= nums[i] <= 109`
#   * `1 <= p <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
using ODE
using Base.Test

tol = 1e-2

solvers = [
           ## Non-stiff
           # fixed step
           ODE.ode1,
           ODE.ode2_midpoint,
           ODE.ode2_heun,
           ODE.ode4,
           # Ode.ode4ms,
           # ODE.ode5ms,
           # adaptive
           ODE.ode21, # this fails on Travis with 0.4?! TODO revert once fixed.
           ODE.ode23,
           ODE.ode45_dp,
           ODE.ode45_fe,
           ODE.ode78,

           ## Stiff
           # fixed-step
           # ODE.ode4s_s,
           # ODE.ode4s_kr,
           # adaptive
           ODE.ode23s]

for solver in solvers
    println("using $solver")
    # dy
    # -- = 6 ==> y = 6t
    # dt
    # we need to fix initstep for the fixed-step methods
    t,y=solver((t,y)->6.0, 0., [0:.1:1;], initstep=.1)
    @test maximum(abs(y-6t)) < tol

    # dy
    # -- = 2t ==> y = t.^2
    # dt
    t,y=solver((t,y)->2t, 0., [0:.001:1;], initstep=0.001)
    @test maximum(abs(y-t.^2)) < tol

    # dy
    # -- = y ==> y = y0*e.^t
    # dt
    t,y=solver((t,y)->y, 1., [0:.001:1;], initstep=0.001)
    @test maximum(abs(y-e.^t)) < tol

    t,y=solver((t,y)->y, 1., [1:-.001:0;], initstep=0.001)
    @test maximum(abs(y-e.^(t-1))) < tol

    # dv       dw
    # -- = -w, -- = v ==> v = v0*cos(t) - w0*sin(t), w = w0*cos(t) + v0*sin(t)
    # dt       dt
    #
    # y = [v, w]
    t,y=solver((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;], initstep=0.001)
    ys = hcat(y...).'   # convert Vector{Vector{Float}} to Matrix{Float}
    @test maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol
end

# rober testcase from http://www.unige.ch/~hairer/testset/testset.html
let
    println("ROBER test case")
    function f(t, y)
        ydot = similar(y)
        ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
        ydot[3] = 3.0e7*y[2]*y[2]
        ydot[2] = -ydot[1] - ydot[3]
        ydot
    end
    t = [0., 1e11]
    t,y = ODE.ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8,
                     maxstep=1e11/10, minstep=1e11/1e18)

    refsol = [0.2083340149701255e-07,
              0.8333360770334713e-13,
              0.9999999791665050] # reference solution at tspan[2]
    @test norm(refsol-y[end], Inf) < 2e-10
end
include("interface-tests.jl")

println("All looks OK")