text:
module BIDSMetadata

export institution_name,
       institution_address,
       institutional_department_name,
       manufacturer,
       manufacturer_model_name,
       device_serial_number,
       station_name,
       software_versions,
       magnetic_field_strength,
       receiver_coil_name,
       receive_coil_active_elements,
       gradient_set_type,
       mr_transmit_coil_sequence,
       matrix_coil_mode,
       coil_combination_method,
       pulse_sequence,
       scanning_sequence,
       sequence_variant,
       scan_options,
       echotime,
       inversiontime,
       slicetime,
       slicedirection,
       dwelltime
###
### Institution information
###

"""
    institution_name(x) -> String

The name of the institution in charge of the equipment that produced the
composite instances.
"""
institution_name(x::Any) = ""
dicomtag(::typeof(institution_name)) = (0x0008,0x0080)
jsontag(::typeof(institution_name)) = "InstitutionName"

"""
    institution_address(x) -> String

The address of the institution in charge of the equipment that produced the
composite instances.
"""
institution_address(x::Any) = ""
dicomtag(::typeof(institution_address)) = (0x0008,0x0081)
jsontag(::typeof(institution_address)) = "InstitutionAddress"

"""
    institutional_department_name(x) -> String

The department in the institution in charge of the equipment that produced the
composite instances.
"""
institutional_department_name(x) = ""
dicomtag(::typeof(institutional_department_name)) = (0x0008,0x1040)
jsontag(::typeof(institutional_department_name)) = "InstitutionalDepartmentName"

###
### Scanner Hardware
###
"""
    manufacturer(x) -> String

Manufacturer of the equipment that produced the composite instances.
"""
manufacturer(x::Any) = ""
dicomtag(::typeof(manufacturer)) = (0x0008, 0x0070)
jsontag(::typeof(manufacturer)) = "Manufacturer"

"""
    manufacturer_model_name(x) -> String

Manufacturer's model name of the equipment that produced the composite
instances.
"""
manufacturer_model_name(x::Any) = ""
dicomtag(::typeof(manufacturer_model_name)) = (0x0008, 0x1090)
jsontag(::typeof(manufacturer_model_name)) = "ManufacturersModelName"

"""
    device_serial_number(x) -> Number

The serial number of the equipment that produced the composite instances.
A pseudonym can also be used to prevent the equipment from being identifiable,
so long as each pseudonym is unique within the dataset.
"""
device_serial_number(x::Any) = 0
dicomtag(typeof(device_serial_number)) = (0x0018, 0x1000)
jsontag(::typeof(device_serial_number)) = "DeviceSerialNumber"

"""
    station_name(x) -> String

Institution defined name of the machine that produced the composite instances.
"""
station_name(x::Any) = ""
dicomtag(::typeof(station_name)) = (0x0008, 0x1010)
jsontag(::typeof(station_name)) = "StationName"

"""
    software_versions(x) -> Vector

Manufacturer’s designation of software version of the equipment that produced
the composite instances.
"""
software_versions(x::Any) = ""
dicomtag(::typeof(software_versions)) = (0x0018, 0x1020)
jsontag(::typeof(software_versions)) = "SoftwareVersions"

"""
    magnetic_field_strength(x) -> Float64

Nominal field strength of MR magnet in Tesla.
"""
magnetic_field_strength(x::Any) = Float64
dicomtag(::typeof(magnetic_field_strength)) = (0x0018,0x0087)
jsontage(::typeof(magnetic_field_strength)) = "MagneticFieldStrength"

"""
    receiver_coil_name(x) -> String

Information describing the receiver coil. Corresponds to DICOM Tag
Receive Coil Name, although not all vendors populate that DICOM Tag, in which
case this field can be derived from an appropriate private DICOM field.
"""
receiver_coil_name(x) = ""
dicomtag(::typeof(receiver_coil_name)) = (0x0018, 0x1250)
jsontag(::typeof(receiver_coil_name)) = "ReceiverCoilName"

"""
    receive_coil_active_elements(x)

Information describing the active/selected elements of the receiver coil. This
doesn’t correspond to a tag in the DICOM ontology. The vendor-defined
terminology for active coil elements can go in this field. As an example, for
Siemens, coil channels are typically not activated/selected individually, but
rather in pre-defined selectable "groups" of individual channels, and the list
of the "groups" of elements that are active/selected in any given scan
populates the Coil String entry in Siemens’ private DICOM fields (e.g., HEA;HEP
for the Siemens standard 32 ch coil when both the anterior and posterior groups
are activated). This is a flexible field that can be used as most appropriate
for a given vendor and coil to define the "active" coil elements. Since
individual scans can sometimes not have the intended coil elements selected, it
is preferable for this field to be populated directly from the DICOM for each
individual scan, so that it can be used as a mechanism for checking that a
given scan was collected with the intended coil elements selected.

* Support suggestion: recommended
"""
function receive_coil_active_elements end
jsontag(::typeof(receive_coil_active_elements)) = "ReceiveCoilActiveElements"

"""
    gradient_set_type(x)

It should be possible to infer the gradient coil from the scanner model. If
not, e.g. because of a custom upgrade or use of a gradient insert set, then the
specifications of the actual gradient coil should be reported independently.

* Support suggestion: recommended
"""
function gradient_set_type end
jsontag(::typeof(gradient_set_type)) = "GradientSetType"

"""
    mr_transmit_coil_sequence(x)

This is a relevant field if a non-standard transmit coil is used.
"""
function mr_transmit_coil_sequence end
dicomtag(::typeof(mr_transmit_coil_sequence)) = (0x0018, 0x9049)
jsontag(::typeof(mr_transmit_coil_sequence)) = "MRTransmitCoilSequence"

"""
    matrix_coil_mode

A method for reducing the number of independent channels by combining in analog
the signals from multiple coil elements. There are typically different default
modes when using un-accelerated or accelerated (e.g. GRAPPA, SENSE) imaging.
"""
function matrix_coil_mode end
jsontag(::typeof(matrix_coil_mode)) = "MatrixCoilMode"

"""
    coil_combination_method

Almost all fMRI studies using phased-array coils use root-sum-of-squares (rSOS)
combination, but other methods exist. The image reconstruction is changed by
the coil combination method (as for the matrix coil mode above), so anything
non-standard should be reported
"""
function coil_combination_method end
jsontag(::typeof(coil_combination_method)) = "CoilCombinationMethod"

# TODO does this mean JSON is in seconds and DICOM is in ms
"""
    echotime(x) -> Seconds

The echo time (TE) for the acquisition.

This parameter is REQUIRED if corresponding fieldmap data is present or the
data comes from a multi echo sequence.

please note that the DICOM term is in milliseconds not seconds
"""
echotime(x::Any) = Float64
dicomtag(::typeof(echotime)) = (0x0018, 0x0081)
jsontag(::typeof(echotime)) = "EchoTime"

"""
    inversiontime

The inversion time (TI) for the acquisition, specified in seconds. Inversion
time is the time after the middle of inverting RF pulse to middle of excitation
pulse to detect the amount of longitudinal magnetization.

please note that the DICOM term is in milliseconds not seconds
"""
function inversiontime end
dicomtag(::typeof(inversiontime)) = (0x0018, 0x0082)
jsontag(::typeof(inversiontime)) = "InversionTime"

"""
    slicetime

The time at which each slice was acquired within each volume (frame) of the
acquisition. Slice timing is not slice order -- rather, it is a list of times
(in JSON format) containing the time (in seconds) of each slice acquisition in
relation to the beginning of volume acquisition. The list goes through the
slices along the slice axis in the slice encoding dimension (see below). Note
that to ensure the proper interpretation of the `SliceTiming` field, it is
important to check if the OPTIONAL `SliceEncodingDirection` exists. In
particular, if `SliceEncodingDirection` is negative, the entries in
`SliceTiming` are defined in reverse order with respect to the slice axis
(i.e., the final entry in the `SliceTiming` list is the time of acquisition of
slice 0). This parameter is REQUIRED for sparse sequences that do not have the
`DelayTime` field set. In addition without this parameter slice time correction
will not be possible.
"""
function slicetiming end

"""
    slicedirection

Possible values: `i`, `j`, `k`, `i-`, `j-`, `k-` (the axis of the NIfTI data
along which slices were acquired, and the direction in which `SliceTiming` is
defined with respect to). `i`, `j`, `k` identifiers correspond to the first,
second and third axis of the data in the NIfTI file. A `-` sign indicates that
the contents of `SliceTiming` are defined in reverse order - that is, the first
entry corresponds to the slice with the largest index, and the final entry
corresponds to slice index zero. When present, the axis defined by
`SliceEncodingDirection` needs to be consistent with the ‘slice_dim’ field in
the NIfTI header. When absent, the entries in `SliceTiming` must be in the
order of increasing slice index as defined by the NIfTI header.                                                                                |
"""
function slicedirection end

"""
    dwelltime -> Second

Actual dwell time (in seconds) of the receiver per point in the readout
direction, including any oversampling. For Siemens, this corresponds to DICOM
field (0019,1018) (in ns). This value is necessary for the optional readout
distortion correction of anatomicals in the HCP Pipelines. It also usefully
provides a handle on the readout bandwidth, which isn’t captured in the other
metadata tags. Not to be confused with `EffectiveEchoSpacing`, and the frequent
mislabeling of echo spacing (which is spacing in the phase encoding direction)
as "dwell time" (which is spacing in the readout direction).                                                                                |
"""
dwelltime(x::Any) = Second(1)
dicomtag(::typeof(dwilltime)) = (0x0019, 0x1018)
jsontag(::typeof(dwilltime)) = "DwellTime"

"""
    flipangle -> Int

Flip angle for the acquisition in degrees.
"""
flipangle(x::Any) = 0
dicomtag(::typeof(flipangle)) = (0x0018,0x1314)
jsontag(::typeof(flipangle)) = "FlipAngle"

"""
    multiband_acceleration_factor

The multiband factor, for multiband acquisitions.
"""
function multiband_acceleration_factor end
jsontag(::typeof(multiband_acceleration_factor)) = "MultibandAccelerationFactor"

"""
    negative_contrast -> Bool
"""
function negative_contrast end
jsontag(::typeof(negative_contrast)) = "NegativeContrast"

"""
    anatomical_landmarks(x) -> Pair

Key:value pairs of any number of additional anatomical landmarks and their
coordinates in voxel units (where first voxel has index 0,0,0) relative to the
associated anatomical MRI, (e.g. {"AC": [127,119,149], "PC": [128,93,141],
"IH": [131,114,206]}, or {"NAS": [127,213,139], "LPA": [52,113,96],
"RPA": [202,113,91]}).
"""
function anatomical_landmarks end
jsontag(::typeof(anatomical_landmarks)) = "AnatomicalLandmarkCoordinates"

"""
    nshots

The number of RF excitations need to reconstruct a slice or volume. Please mind
that this is not the same as Echo Train Length which denotes the number of
lines of k-space collected after an excitation.                                                                                |
"""

"""
    parallel_reduction_factor_in_plane

The parallel imaging (e.g, GRAPPA) factor. Use the denominator of the fraction
of k-space encoded for each slice. For example, 2 means half of k-space is
encoded. Corresponds to DICOM Tag 0018, 9069 `Parallel Reduction Factor
In-plane`.                                                                                |

    parallel_reduction_factor_inplane

Measurement time reduction factor expressed as ratio of original and reduced
measurement time for the in-plane direction. Required if Frame Type (0008,9007)
Value 1 of this frame is ORIGINAL and Parallel Acquisition (0018,9077) equals
YES. Otherwise may be present if Frame Type (0008,9007) Value 1 is DERIVED and
Parallel Acquisition (0018,9077) equals YES.
"""
function parallel_reduction_factor_inplane end
dicomtag(::typeof(parallel_reduction_factor_inplane)) = (0x0018,0x9069)

"""
    parallel_acquisition_technique

The type of parallel imaging used (e.g. GRAPPA, SENSE). Corresponds to DICOM
Tag 0018, 9078 `Parallel Acquisition Technique`.                                                                                |
"""

"""
    partial_fourier

The fraction of partial Fourier information collected. Corresponds to DICOM
Tag 0018, 9081 `Partial Fourier`.                                                                                |
"""
function partial_fourier end

"""
    partial_fourier_direction

The direction where only partial Fourier information was collected. Corresponds to DICOM Tag 0018, 9036 `Partial Fourier Direction`.                                                                                |
"""

"""
    phase_encoding_direction

Possible values: `i`, `j`, `k`, `i-`, `j-`, `k-`. The letters `i`, `j`, `k`
correspond to the first, second and third axis of the data in the NIFTI file.
The polarity of the phase encoding is assumed to go from zero index to maximum
index unless `-` sign is present (then the order is reversed - starting from
the highest index instead of zero). `PhaseEncodingDirection` is defined as the
direction along which phase is was modulated which may result in visible
distortions. Note that this is not the same as the DICOM term
`InPlanePhaseEncodingDirection` which can have `ROW` or `COL` values. This
parameter is REQUIRED if corresponding fieldmap data is present or when using
multiple runs with different phase encoding directions (which can be later used
for field inhomogeneity correction).
"""
function phase_encoding_direction end

"""
    effective_echo_spacing

The "effective" sampling interval, specified in seconds, between lines in the
phase-encoding direction, defined based on the size of the reconstructed image
in the phase direction. It is frequently, but incorrectly, referred to as
"dwell time" (see `DwellTime` parameter below for actual dwell time). It is
required for unwarping distortions using field maps. Note that beyond just
in-plane acceleration, a variety of other manipulations to the phase encoding
need to be accounted for properly, including partial fourier, phase
oversampling, phase resolution, phase field-of-view and interpolation.2
This parameter is REQUIRED if corresponding fieldmap data is present.                                                                                |

2Conveniently, for Siemens’ data, this value is easily obtained as
1/\[`BWPPPE` \* `ReconMatrixPE`\], where BWPPPE is the
"`BandwidthPerPixelPhaseEncode` in DICOM tag (0019,1028) and ReconMatrixPE is
the size of the actual reconstructed data in the phase direction (which is NOT
reflected in a single DICOM tag for all possible aforementioned scan
manipulations). See [here](https://lcni.uoregon.edu/kb-articles/kb-0003) and
[here](https://github.com/neurolabusc/dcm_qa/tree/master/In/TotalReadoutTime)

"""
function effective_echo_spacing end

"""
    total_readout_time

This is actually the "effective" total readout time , defined as the readout
duration, specified in seconds, that would have generated data with the given
level of distortion. It is NOT the actual, physical duration of the readout
train. If `EffectiveEchoSpacing` has been properly computed, it is just
`EffectiveEchoSpacing * (ReconMatrixPE - 1)`.3 .

* This parameter is
REQUIRED if corresponding "field/distortion" maps acquired with opposing phase
encoding directions are present (see 8.9.4).                                                                                |

3We use the "FSL definition", i.e, the time between the center of the
first "effective" echo and the center of the last "effective" echo.
"""
function total_readout_time end

"""
    ntemporal_positions

Total number of temporal positions prescribed.
"""
function ntemporal_positions end
dicomtag(::typeof(ntemporal_positions)) = (0x0020, 0x0105)

"""
    temporal_resolution

Time delta between Images in a dynamic of functional set of Images.
"""
function temporal_resolution end
dicomtag(::typeof(temporal_resolution)) = (0x0020, 0x0110)

# TODO Pulse sequence type could probably have specific type parameters but for
# now this is just going to return String
"""
    pulse_sequence(x) -> String

A general description of the pulse sequence used for the scan (i.e. MPRAGE,
Gradient Echo EPI, Spin Echo EPI, Multiband gradient echo EPI).                                                                                |
"""
pulse_sequence(x::Any) = ""

"""
    scanning_sequence

Description of the type of data acquired. Corresponds to DICOM Tag 0018, 0020
`Scanning Sequence`.                                                                                |
"""
scanning_sequence(x::Any) = ""

"""
    sequence_variant

Variant of the ScanningSequence. Corresponds to DICOM Tag 0018, 0021 `Sequence
Variant`.                                                                                |
"""
function sequence_variant end

"""
    scan_options

Parameters of ScanningSequence. Corresponds to DICOM Tag 0018, 0022
`Scan Options`.                                                                                |
"""
function scan_options end

"""
    sequence_name

Manufacturer’s designation of the sequence name. Corresponds to DICOM Tag 0018,
0024 `Sequence Name`.                                                                                |
"""
function sequence_name end

"""
    pulse_sequence_details

Information beyond pulse sequence type that identifies the specific pulse
sequence used (i.e. "Standard Siemens Sequence distributed with the VB17
software," "Siemens WIP ### version #.##," or "Sequence written by X using a
version compiled on MM/DD/YYYY"). |
"""
function pulse_sequence_details end

"""
    nonlinear_gradient_correction(x) -> Bool

Stating if the image saved has been corrected for gradient nonlinearities by
the scanner sequence.                                                                                |
"""
function nonlinear_gradient_correction(x) end

# TODO does this mean JSON is in seconds and DICOM is in ms
"""
    echotime(x) -> Millisecond

The echo time (TE) for the acquisition.

This parameter is REQUIRED if corresponding fieldmap data is present or the
data comes from a multi echo sequence.

please note that the DICOM term is in milliseconds not seconds
"""
echotime(x::Any) = Millisecond(1)
dicomtag(::typeof(echotime)) = (0x0018, 0x0081)
jsontag(::typeof(echotime)) = "EchoTime"

"""
    inversiontime(x) -> Second

The inversion time (TI) for the acquisition, specified in seconds. Inversion
time is the time after the middle of inverting RF pulse to middle of excitation
pulse to detect the amount of longitudinal magnetization. Corresponds to DICOM
Tag 0018, 0082 `Inversion Time` (please note that the DICOM term is in
milliseconds not seconds).                                                                                |
"""
inversiontime(x::Any) = Second(1)
dicomtag(::typeof(inversiontime)) = (0x0018, 0x0082)
jsontag(::typeof(inversiontime)) = "InversionTime"

"""
    slicetime

The time at which each slice was acquired within each volume (frame) of the
acquisition. Slice timing is not slice order -- rather, it is a list of times
(in JSON format) containing the time (in seconds) of each slice acquisition in
relation to the beginning of volume acquisition. The list goes through the
slices along the slice axis in the slice encoding dimension (see below). Note
that to ensure the proper interpretation of the `SliceTiming` field, it is
important to check if the OPTIONAL `SliceEncodingDirection` exists. In
particular, if `SliceEncodingDirection` is negative, the entries in
`SliceTiming` are defined in reverse order with respect to the slice axis
(i.e., the final entry in the `SliceTiming` list is the time of acquisition of
slice 0). This parameter is REQUIRED for sparse sequences that do not have the
`DelayTime` field set. In addition without this parameter slice time correction
will not be possible.
"""
function slicetiming end

"""
    slicedirection

Possible values: `i`, `j`, `k`, `i-`, `j-`, `k-` (the axis of the NIfTI data
along which slices were acquired, and the direction in which `SliceTiming` is
defined with respect to). `i`, `j`, `k` identifiers correspond to the first,
second and third axis of the data in the NIfTI file. A `-` sign indicates that
the contents of `SliceTiming` are defined in reverse order - that is, the first
entry corresponds to the slice with the largest index, and the final entry
corresponds to slice index zero. When present, the axis defined by
`SliceEncodingDirection` needs to be consistent with the ‘slice_dim’ field in
the NIfTI header. When absent, the entries in `SliceTiming` must be in the
order of increasing slice index as defined by the NIfTI header.                                                                                |
"""
function slicedirection end

"""
    dwelltime -> Nanosecond

Actual dwell time (in seconds) of the receiver per point in the readout
direction, including any oversampling. For Siemens, this corresponds to DICOM
field (0019,1018) (in ns). This value is necessary for the optional readout
distortion correction of anatomicals in the HCP Pipelines. It also usefully
provides a handle on the readout bandwidth, which isn’t captured in the other
metadata tags. Not to be confused with `EffectiveEchoSpacing`, and the frequent
mislabeling of echo spacing (which is spacing in the phase encoding direction)
as "dwell time" (which is spacing in the readout direction).                                                                                |
"""
dwelltime(x::Any) = Nanosecond(1)
dicomtag(::typeof(dwilltime)) = (0x0019, 0x1018)
jsontag(::typeof(dwilltime)) = "DwellTime"

"""
    flipangle

Flip angle for the acquisition
"""
function flipangle end
dicomtag(::typeof(flipangle)) = (0x0018,0x1314)
jsontag(::typeof(flipangle)) = "FlipAngle"

"""
    multiband_acceleration_factor

The multiband factor, for multiband acquisitions.
"""
function multiband_acceleration_factor end
jsontag(::typeof(multiband_acceleration_factor)) = "MultibandAccelerationFactor"

"""
    negative_contrast -> Bool
"""
function negative_contrast end
jsontag(::typeof(negative_contrast)) = "NegativeContrast"

"""
    contrast_bolus_ingredient(x) = ""

Active ingredient of agent. Values MUST be one of: IODINE, GADOLINIUM,
CARBON DIOXIDE, BARIUM, XENON.
"""
contrast_bolus_ingredient(x::Any) = ""
dicomtag(::typeof(contrast_bolus_ingredient)) = (0x0018,0x1048)
jsontag(::typeof(contrast_bolus_ingredient)) = "ContrastBolusIngredient"

"""
    reptime(x) -> Float64

The time in seconds between the beginning of an acquisition of one volume and
the beginning of acquisition of the volume following it (TR). Please note that
this definition includes time between scans (when no data has been acquired) in
case of sparse acquisition schemes. This value needs to be consistent with the
pixdim[4] field (after accounting for units stored in xyzt_units field) in the
NIfTI header. This field is mutually exclusive with VolumeTiming and is derived
from DICOM Tag 0018, 0080 and converted to seconds.
"""
reptime(x::Any) -> one(Float64)
dicomtag(::typeof(reptime)) = (0x0018,0x0080)
jsontag(::typeof(reptime)) = "RepetitionTime"

"""
    voltime(x)

The time at which each volume was acquired during the acquisition. It is
described using a list of times (in JSON format) referring to the onset of each
volume in the BOLD series. The list must have the same length as the BOLD
series, and the values must be non-negative and monotonically increasing. This
field is mutually exclusive with RepetitionTime and DelayTime. If defined, this
requires acquisition time (TA) be defined via either SliceTiming or
AcquisitionDuration be defined.
"""
function voltime end
jsontag(::typeof(voltime)) = "VolumeTiming"

"""
    nvol_scanner_discarded(x) -> Int

Number of volumes ("dummy scans") discarded by the scanner (as opposed to those
discarded by the user post hoc) before saving the imaging file. For example, a
sequence that automatically discards the first 4 volumes before saving would
have this field as 4. A sequence that doesn't discard dummy scans would have
this set to 0. Please note that the onsets recorded in the _event.tsv file
should always refer to the beginning of the acquisition of the first volume in
the corresponding imaging file - independent of the value of
NumberOfVolumesDiscardedByScanner field.
"""
nvol_scanner_discarded(x::Any) = 0
jsontag(::typeof(nvol_scanner_discarded)) = "NumberOfVolumesDiscardedByScanner"

"""
    nvol_user_discarded(x) -> Int

Number of volumes ("dummy scans") discarded by the user before including the
file in the dataset. If possible, including all of the volumes is strongly
recommended. Please note that the onsets recorded in the _event.tsv file
should always refer to the beginning of the acquisition of the first volume in
the corresponding imaging file - independent of the value of
NumberOfVolumesDiscardedByUser field.
"""
nvol_user_discarded(x::Any) = 0
jsontag(::typeof(nvol_user_discarded)) = "NumberOfVolumesDiscardedByUser"

"""
    taskname(x) -> String

Name of the task. No two tasks should have the same name. Task label (task-)
included in the file name is derived from this field by removing all non
alphanumeric ([a-zA-Z0-9]) characters. For example task name faces n-back will
corresponds to task label facesnback. A RECOMMENDED convention is to name
resting state task using labels beginning with rest.
"""
taskname(x::Any) = ""

"""
    taskinstructions(x) -> String
"""
taskinstructions(x::Any) = ""

"""
    taskdescription(x) -> String
"""
taskdescription(x::Any) = ""

"""
    cognitive_atlas_id(x)
"""
function cognitive_atlas_id end
const cogatlas = cognitive_atlas_id

"""
    cognitive_atlas_paradigm_ontology
"""
function cognitive_atlas_paradigm_ontology end
const cogpo = cognitive_atlas_paradigm_ontology

end


text:
module DynamicGrids
# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    # Use [`XX`](@ref) in the docs but not the readme
    text = replace(read(path, String), r"`(\w+\w)`" => s"[`\1`](@ref)")
    # Use doctests
    replace(text, "```julia" => "```jldoctest")
end DynamicGrids

using Colors,
      ConstructionBase,
      Crayons,
      DocStringExtensions,
      FieldDefaults,
      FieldMetadata,
      FieldDocTables,
      FileIO,
      Mixers,
      OffsetArrays,
      REPL,
      Setfield,
      UnicodeGraphics

using Base: tail
using Lazy: @forward

import Base: show, getindex, setindex!, lastindex, size, length, push!, append!,
              broadcast, broadcast!, similar, eltype, iterate

import FieldMetadata: @description, description, @limits, limits,
                      @flattenable, flattenable, default

export sim!, resume!, replay

export savegif, show_frame, ruletypes

export distances, broadcastable_indices, sizefromradius

export AbstractSimData, SimData, MultiSimData

export Rule, PartialRule, NeighborhoodRule, PartialNeighborhoodRule, CellRule, Chain

export Interaction, NeighborhoodInteraction, CellInteraction,
       PartialInteraction, PartialNeighborhoodInteraction

export AbstractRuleset, Ruleset, MultiRuleset

export Life

export Neighborhood, AbstractRadialNeighborhood, RadialNeighborhood,
       AbstractCustomNeighborhood, CustomNeighborhood, LayeredCustomNeighborhood,
       VonNeumannNeighborhood

export RemoveOverflow, WrapOverflow

export Output, GraphicOutput, ImageOutput, ArrayOutput, REPLOutput

export FrameProcessor, ColorProcessor, Greyscale, Grayscale

export CharStyle, Block, Braile

const FIELDDOCTABLE = FieldDocTable((:Description, :Default, :Limits),
                                    (description, default, limits);
                                    truncation=(100,40,100))

# Documentation templates
@template TYPES =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

include("rules.jl")
include("interactions.jl")
include("chain.jl")
include("rulesets.jl")
include("simulationdata.jl")
include("outputs/output.jl")
include("outputs/graphic.jl")
include("outputs/image.jl")
include("outputs/array.jl")
include("outputs/repl.jl")
include("interface.jl")
include("framework.jl")
include("sequencerules.jl")
include("maprules.jl")
include("mapinteractions.jl")
include("neighborhoods.jl")
include("utils.jl")
include("life.jl")

end

ext: jl


text:
module TestGrid

using Test
import LinearAlgebra: norm, norm2
import GraphIdx
import GraphIdx.Grid: Pixel, line_D, incmat
import GraphIdx.Grid: GridGraph, iter_edges, iter_edges_pixel
import GraphIdx: num_edges, num_nodes

@testset "Pixel                          " begin
    p0 = Pixel(0, 0)
    @test norm(p0) == 0.0
    @test norm2(p0) == 0.0

    p1 = Pixel(1, 1)
    @test norm2(p1) == 2.0
    @test norm(p1) ≈ √2

    p3 = Pixel(1, 3)
    @test norm2(p3) == 10
    @test norm(p3)  ≈ √10
end

@testset "line_D                         " begin
    D = line_D(3)
    @test D ≈ [1  -1   0
               0  +1  -1]
    D = line_D(5)
    @test D ≈ [1  -1   0   0   0
               0  +1  -1   0   0
               0   0  +1  -1   0
               0   0   0  +1  -1]

    for n in [11, 13, 42]
        @test line_D(n) == incmat(n, 1)
        @test line_D(n) == incmat(1, n)
    end
end

@testset "GridGraph: 3x2                 " begin
    g = GridGraph(3, 2)
    @test (g |> collect).edges == [(1, 2), (4, 5), (2, 3), (5, 6), (1, 4), (2, 5), (3, 6)]
end

@testset "GridGraph: 2x3                 " begin
    g = GridGraph(2, 3)
    m = 2*2 + 3
    n = 6
    @test num_nodes(g) == n
    @test num_edges(g) == m
    @test (g |> collect).edges == [
        (1, 2), (3, 4), (5, 6), (1, 3), (3, 5), (2, 4), (4, 6)
    ]
    @testset "iter edge pixels" begin
        local k = 0
        iter_edges_pixel(g) do i1, j1, i2, j2, lam::Float64
            k += 1
            @test lam ≈ 1
            if k == 1
                @test (i1, j1) == (1, 1)
                @test (i2, j2) == (2, 1)
            elseif k == 4
                @test (i1, j1) == (1, 1)
                @test (i2, j2) == (1, 2)
            elseif k == 5
                @test (i1, j1) == (1, 2)
                @test (i2, j2) == (1, 3)
            end
        end
        @test k == m
    end

    @testset "iter edges" begin
        local k = 0
        iter_edges(g) do u, v, lam::Float64
            k += 1
            if k == 1
                @test u == 1
                @test v == 2
            elseif k == 4
                @test u == 1
                @test v == 3
            elseif k == 5
                @test u == 3
                @test v == 5
            else
                @test u < v
            end
            @test lam ≈ 1
        end
        @test k == m
    end
end

@testset "GridGraph: 3x7                 " begin
    """
     1  4  7  10  13  16  19
     2  5  8  11  14  17  20
     3  6  9  12  15  18  21
    """
    g = GridGraph(3, 7)
    edge_graph = g |> collect
    @test edge_graph.edges == [
        (1, 2),
        (4, 5),
        (7, 8),
        (10, 11),
        (13, 14),
        (16, 17),
        (19, 20),
        (2, 3),
        (5, 6),
        (8, 9),
        (11, 12),
        (14, 15),
        (17, 18),
        (20, 21),
        (1, 4),
        (4, 7),
        (7, 10),
        (10, 13),
        (13, 16),
        (16, 19),
        (2, 5),
        (5, 8),
        (8, 11),
        (11, 14),
        (14, 17),
        (17, 20),
        (3, 6),
        (6, 9),
        (9, 12),
        (12, 15),
        (15, 18),
        (18, 21)
    ]
    idx = edge_graph |> e -> GraphIdx.IncidenceIndex(GraphIdx.num_nodes(g), e.edges)
    @test idx[1] == [(2, 1), (4, 15)]
    @test idx[5] == [(4, 2), (6, 9), (2, 21), (8, 22)]
    @testset "compare to enumerate(GridGraph)" begin
        edges = Set{NTuple{2,Int}}()
        GraphIdx.enumerate_edges(g) do e::Int, u::Int, v::Int
            push!(edges, (u, v))
        end
        @test edges == Set(edge_graph.edges)
    end
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# MemoryFunction.jl
# References:
# [1] R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Mobility of slow electrons in a polar crystal, PhysicalReview127, 1004 (1962) https://doi.org/10.1103/PhysRev.127.1004.
# [2] Devreese, J. De Sitter, and M. Goovaerts, “Optical absorption of polarons in thefeynman-hellwarth-iddings-platzman approximation,”Phys. Rev. B, vol. 5, pp. 2367–2381, Mar 1972.
# [3] F. Peeters and J. Devreese, “Theory of polaron mobility,” inSolid State Physics,pp. 81–133, Elsevier, 1984.

"""
    polaron_memory_function(Ω, β, α, v, w; ω = 1.0, rtol = 1e-3)

Calculate the memory function χ(Ω) of the polaron at finite temperatures (equation (35) in FHIP 1962) for a given frequency Ω.

Includes conditionals for zero-temperature and DC limits.

# Arguments
- `Ω::Float64`: electric field frequency (2π THz).
- `β::Union{Float64, Vector{Float64}` reduced thermodynamic temperatures (unitless). Either a single value for one phonon frequency or a vector of values for multiple phonon frequencies.
- `α::Union{Float64, Vector{Float64}}`: Frohlich alpha coupling parameter (unitless). Either a single value for one phonon frequency or a vector of values for multiple phonon frequencies.
- `v::Union{Float64, Vector{Float64}}`: variational 'v' parameters that minimise the polaron free energy (unitless).
- `w::Union{Float64, Vector{Float64}}`: variational 'w' parameters that minimise the polaron free energy (unitless).
- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz). Predefined as `ω = 1.0` for a single mode in polaron units.
- `rtol = 1e-3`: specifies the relative error tolerance for the QuadGK integral.

See R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Mobility of slow electrons in a polar crystal, PhysicalReview127, 1004 (1962). https://doi.org/10.1103/PhysRev.127.1004.

See also [`polaron_memory_function_thermal`](@ref), [`polaron_memory_function_athermal`](@ref), [`polaron_memory_function_dc`](@ref).
"""
function polaron_memory_function(Ω, β, α, v, w; ω = 1.0, rtol = 1e-3)

    # Zero temperature and frequency is just zero.
	if Ω == 0 && any(x -> x == Inf64, β)
		return 0.0 + im * 0.0

    # DC zero frequency limit at finite temperatures.
	elseif Ω == 0 && any(x -> x != Inf64, β)
		return polaron_memory_function_dc(β, α, v, w; ω = ω, rtol = rtol)

    # Zero temperature limit at AC finite frequencies.
	elseif Ω != 0 && any(x -> x == Inf64, β)
		return polaron_memory_function_athermal(Ω, α, v, w; ω = ω, rtol = rtol)

    # Finite temperatures and frequencies away from zero limits.
	elseif Ω != 0 && any(x -> x != Inf64, β)
		return polaron_memory_function_thermal(Ω, β, α, v, w; ω = ω, rtol = rtol)

    # Any other frequencies or temperatures (e.g. negative or complex) prints error message.
	else
		println("Photon frequency Ω and thermodynamic temperature β must be ≥ 0.0.")
	end
end

"""
----------------------------------------------------------------------
Multiple Branch Polaron Memory Function and Complex Conductivity
----------------------------------------------------------------------

This section of the code is dedicated to calculating the polaron memory function and complex conductivity, generalised from FHIP's expression to the case where multiple phonon modes are present in the material.
"""

"""
    function polaron_memory_function_thermal(Ω, β, α, v, w; ω = 1.0, rtol = 1e-3)

Calculate polaron complex memory function at finite temperature inclusive of multiple phonon branches j, each with angular frequency ω[j] (2π THz).

# Arguments
- `Ω::Float64`: is the frequency (THz) of applied electric field.
- `β::Vector{Float64}`: is an array of reduced thermodynamic betas, one for each phonon frequency ω[j].
- `α::Vector{Float64}`: is an array of decomposed Frohlich alphas, one for each phonon frequency ω[j].
- `v::Union{Float64, Vector{Float64}}`: is a vector or single value of the 'v' variational parameter(s).
- `w::Union{Float64, Vector{Float64}}`: is a vector or single value of the 'w' variational parameter(s).
- `m_eff::Float64`: is the bare-band mass of the particle (typically electron / hole, in units of electron mass m_e).
- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz). Predefined as `ω = 1.0` for a single mode in polaron units.
- `rtol = 1e-3`: specifies the relative error tolerance for the QuadGK integral.

See also [`polaron_memory_function`](@ref).
"""
function polaron_memory_function_thermal(Ω, β, α, v, w; ω = 1.0, rtol = 1e-3)

    # FHIP1962, page 1009, eqn (36).
    S(t, β) = (exp(1im * t) + exp(-1im * t - β)) / (1 - exp(-β)) / D_j(-1im * t, β, v, w)^(3 / 2)

    # FHIP1962, page 1009, eqn (35a).
    integrand(t, β, Ω) = (1 - exp(1im * Ω * 2π * t)) * imag(S(t, β))

    integral = map((x, y) -> quadgk(t -> integrand(t, x, Ω / y), 0.0, Inf64, rtol = rtol)[1], β, ω)

     memory = sum(2 .* α .* ω .^2 .* integral ./ (3 * √π * Ω * 2π))

    return memory
end

"""
    function polaron_memory_function_athermal(Ω, α, v, w; ω = 1.0, rtol = 1e-3)

Calculate polaron complex memory function at zero-temperature inclusive of multiple phonon branches j, each with angular frequency ω[j] (2π THz).

# Arguments
- `Ω::Float64`: is the frequency (THz) of applied electric field.
- `α::Vector{Float64}`: is an array of decomposed Frohlich alphas, one for each phonon frequency ω[j].
- `v::Union{Float64, Vector{Float64}}`: is a vector or single value of the 'v' variational parameter(s).
- `w::Union{Float64, Vector{Float64}}`: is a vector or single value of the 'w' variational parameter(s).
- `m_eff::Float64`: is the bare-band mass of the particle (typically electron / hole, in units of electron mass m_e).
- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz). Predefined as `ω = 1.0` for a single mode in polaron units.
- `rtol = 1e-3`: specifies the relative error tolerance for the QuadGK integral.

See also [`polaron_memory_function`](@ref).
"""
function polaron_memory_function_athermal(Ω, α, v, w; ω = 1.0, rtol = 1e-3)
    # FHIP1962, page 1009, eqn (36).
    S(t) = exp(1im * t) / D_j(-1im * t, v, w)^(3 / 2)

    # FHIP1962, page 1009, eqn (35a).
    integrand(t, Ω) = (1 - exp(1im * 2π * Ω * t)) * imag(S(t))

    memory = 0.0

    for j in 1:length(ω)

        integral = quadgk(t -> integrand(t, Ω / ω[j]), 0.0, 1 / rtol, rtol = rtol)[1]
        memory += 2 * α[j] * ω[j]^2 * integral / (3 * √π * Ω * 2π)

    end

    return memory
end

"""
    function polaron_memory_function_dc(β, α, v, w; ω = 1.0, rtol = 1e-3)

Calculate zero-frequency polaron complex memory function at finite temperature inclusive of multiple phonon branches j, each with angular frequency ω[j] (2π THz).

# Arguments
- `β::Vector{Float64}`: is an array of reduced thermodynamic betas, one for each phonon frequency ω[j].
- `α::Vector{Float64}`: is an array of decomposed Frohlich alphas, one for each phonon frequency ω[j].
- `v::Union{Float64, Vector{Float64}}`: is a vector or single value of the 'v' variational parameter(s).
- `w::Union{Float64, Vector{Float64}}`: is a vector or single value of the 'w' variational parameter(s).
- `m_eff::Float64`: is the bare-band mass of the particle (typically electron / hole, in units of electron mass m_e).
- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz). Predefined as `ω = 1.0` for a single mode in polaron units.
- `rtol = 1e-3`: specifies the relative error tolerance for the QuadGK integral.

See also [`polaron_memory_function`](@ref).
"""
function polaron_memory_function_dc(β, α, v, w; ω = 1.0, rtol = 1e-3)
    # FHIP1962, page 1009, eqn (36).
    S(t, β) = (exp(1im * t) + exp(-1im * t - β)) / (1 - exp(-β)) / D_j(-1im * t, β, v, w)^(3 / 2)

    # FHIP1962, page 1009, eqn (35a).
    integrand(t, β) = -im * t * imag(S(t, β))

    integral = map(x -> quadgk(t -> integrand(t, x), 0.0, 1 / rtol, rtol = rtol)[1], β)

    memory = sum(2 .* α .* ω .* integral ./ (3 * √π))

    return memory
end


text:
using WallTimeProgress
tracker = WallTimeTracker(10)   # output every 10 periods
increment!(tracker)
sleep(1)
increment!(tracker, 9)
count(tracker)


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

using Markdown
using InteractiveUtils

# ╔═╡ e637f8ae-7e0b-11eb-31d9-c5343065be69
md"""
The Heston volatility model follows

$dS_{t} = S_{t} = (r - q)dt + \sqrt{V_{t}} dW_{t}$
$dV_{t} = k(\nu - V_{t}) dt + \lambda \sqrt{\nu_{t}}dZ_{t}$
$Cov(W_{t}, Z_{t}) = \rho dt, \rho \in (-1, 1)$
+ Initial volatility : $\sigma_{0} > 0$
+ Equilibrium volatility : $\nu$
+ Speed of mean reversion $k$
+ vol of vol: $\lambda$ . If $2k\nu > \lambda^2$ variance cannot reach zero (Feller condition)
+ Variance $V_{t}$  is is non-centrally chi-squared distributed, \sqrt{V_{t}} is Rayleigh Distributed
+ Log-Price: $\phi(u, t) := \mathbb{E}[exp(iu \log(S_{t}) | S_{0}, \sigma_{0}^{2}]$

\begin{align}
\phi(u, t)  &=  exp(iu(logS0+ (r−q)t))×exp(ηκλ−2((κ−ρλiu+d)t−2 \\
            & \times log((1−g1edt)/(1−g1))))×exp(σ20λ−2(κ−ρλiu+d)(1−edt)/(1−g1edt))\\
\end{align}

Fast-Reversion Limit of the Heston Model

https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2418631

"""

# ╔═╡ Cell order:
# ╠═e637f8ae-7e0b-11eb-31d9-c5343065be69


text:
@with_kw mutable struct DecreasingParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    γ::T = 1.
end

struct Decreasing <: AbstractStep
    params::DecreasingParameters{Float64}

    function (::Type{Decreasing})(; kwargs...)
        return new(DecreasingParameters(; kwargs...))
    end
end

function stepsize(decreasing::Decreasing,klocal::Integer,kglobal::Integer,fval::Real,x::AbstractVector,g::AbstractVector)
    @unpack γ = decreasing.params
    return γ/√kglobal
end


text:
# include("cooling.jl")
include("adiabatic.jl")

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
# ---
# title: 828. Count Unique Characters of All Substrings of a Given String
# id: problem828
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Two Pointers
# link: <https://leetcode.com/problems/count-unique-characters-of-all-substrings-of-a-given-string/description/>
# hidden: true
# ---
#
# Let's define a function `countUniqueChars(s)` that returns the number of
# unique characters on `s`, for example if `s = "LEETCODE"` then `"L"`,
# `"T"`,`"C"`,`"O"`,`"D"` are the unique characters since they appear only once
# in `s`, therefore `countUniqueChars(s) = 5`.
#
# On this problem given a string `s` we need to return the sum of
# `countUniqueChars(t)` where `t` is a substring of `s`. Notice that some
# substrings can be repeated so on this case you have to count the repeated ones
# too.
#
# Since the answer can be very large, return the answer modulo `10 ^ 9 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "ABC"
#     Output: 10
#     Explanation: All possible substrings are: "A","B","C","AB","BC" and "ABC".
#     Evey substring is composed with only unique letters.
#     Sum of lengths of all substring is 1 + 1 + 1 + 2 + 2 + 3 = 10
#
#
# **Example 2:**
#
#
#
#     Input: s = "ABA"
#     Output: 8
#     Explanation: The same as example 1, except countUniqueChars("ABA") = 1.
#
#
# **Example 3:**
#
#
#
#     Input: s = "LEETCODE"
#     Output: 92
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 10^4`
#   * `s` contain upper-case English letters only.
#
#
## @lc code=start
using LeetCode

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


text:
#######################
# seeded reverse pass #
#######################

# derivative (input is an scalar, output is an array/scalar) #
#------------------------------------------------------------#

function seeded_reverse_pass!(result::AbstractArray, output::AbstractArray, input::TrackedReal, tape)
    result_vector = reshape(result, length(output))
    for i in eachindex(output)
        result_vector[i] = seeded_reverse_pass!(output[i], input, tape)
    end
    return result
end

function seeded_reverse_pass!(output::TrackedReal, input::TrackedReal, tape)
    pull_value!(output)
    unseed!(input)
    seed!(output)
    reverse_pass!(tape)
    return deriv(input)
end

# gradient (input is an array, output is a scalar) #
#--------------------------------------------------#

function seeded_reverse_pass!(result, output::TrackedReal, input, tape)
    pull_value!(output)
    unseed!(input)
    seed!(output)
    reverse_pass!(tape)
    extract_result!(result, output, input)
    return result
end

function seeded_reverse_pass!(result, output::Number, input, tape)
    extract_result!(result, output)
    return result
end

# jacobian (input and output are both arrays) #
#---------------------------------------------#

function seeded_reverse_pass!(result::AbstractArray, output::AbstractArray, input::TrackedArray, tape)
    result_matrix = reshape(result, length(output), length(input))
    input_deriv = deriv(input)
    pull_value!(output)
    for i in eachindex(output)
        unseed!(input)
        seed!(output, i)
        reverse_pass!(tape)
        for j in eachindex(input)
            result_matrix[i, j] = input_deriv[j]
        end
        unseed!(output, i)
    end
    return result
end

function seeded_reverse_pass!(result::DiffResult, output::AbstractArray, input::TrackedArray, tape)
    seeded_reverse_pass!(DiffResults.jacobian(result), output, input, tape)
    extract_result_value!(result, output)
    return result
end

function seeded_reverse_pass!(result::Tuple, output::AbstractArray, input::Tuple, tape)
    for i in eachindex(result)
        seeded_reverse_pass!(result[i], output, input[i], tape)
    end
    return result
end

#####################
# result extraction #
#####################

function extract_result!(result::Tuple, output, input::Tuple)
    for i in eachindex(result)
        extract_result!(result[i], output, input[i])
    end
    return result
end

function extract_result!(result::Tuple, output)
    for i in eachindex(result)
        extract_result!(result[i], output)
    end
    return result
end

function extract_result!(result::AbstractArray, output::TrackedReal, input::TrackedArray)
    copyto!(result, deriv(input))
    return result
end

function extract_result!(result::DiffResult, output::TrackedReal, input::TrackedArray)
    result = DiffResults.value!(result, value(output))
    copyto!(DiffResults.gradient(result), deriv(input))
    return result
end

function extract_result!(result::AbstractArray, output::Number)
    fill_zeros!(result)
    return result
end

function extract_result!(result::DiffResult, output::Number)
    result = DiffResults.value!(result, output)
    fill_zeros!(DiffResults.gradient(result))
    return result
end

function extract_result_value!(result::Tuple, output)
    for i in eachindex(result)
        extract_result_value!(result[i], output)
    end
    return result
end

function extract_result_value!(result::DiffResult, output::AbstractArray)
    result = DiffResults.value!(value, result, output)
    return result
end

function extract_result_value!(result::DiffResult, output::TrackedArray)
    result = DiffResults.value!(result, value(output))
    return result
end

function extract_result_value!(result::AbstractArray, output::AbstractArray)
    map!(value, result, output)
    return result
end

function extract_result_value!(result::AbstractArray, output::TrackedArray)
    copyto!(result, value(output))
    return result
end

fill_zeros!(result::AbstractArray) = fill!(result, zero(eltype(result)))

#######################
# result construction #
#######################

construct_result(output::AbstractArray, input::Tuple) = map(x -> construct_result(output, x), input)

construct_result(output::AbstractArray, input::TrackedArray) = similar(deriv(input), length(output), length(input))

construct_result(input::TrackedArray) = similar(deriv(input))

construct_result(input::Tuple) = map(construct_result, input)


text:
using RegionTrees
using StaticArrays
using PyPlot
using Statistics
using LinearAlgebra
using SparseArrays
using VectorizedRoutines
using ProgressBars
using MAT

struct ray3d{T}
    origin::T
    dir::T
    inv_direction::T
    sign::Union{Array{Int64,1},SArray{Tuple{3},Int64,1,3}}
end

function ray3d(p1::T, p2::T) where {T}
    dir::T = p2 - p1
    invdir::T = 1.0 ./ dir
    sign = (invdir .< 0) .* 1
    return ray3d(p1, dir, invdir, sign)
end

#An Efficient and Robust Ray–Box Intersection Algorithm by
#Amy Williams Steve Barrus R.Keith Morley Peter Shirley.
@inline @inbounds function intersects(
    r::ray3d,
    b::HyperRectangle{3,Float64},
)::Float64
    Q = [b.origin, b.origin + b.widths]
    tmin = (Q[r.sign[1]+1][1] - r.origin[1]) * r.inv_direction[1]
    tmax = (Q[1-r.sign[1]+1][1] - r.origin[1]) * r.inv_direction[1]
    tymin = (Q[r.sign[2]+1][2] - r.origin[2]) * r.inv_direction[2]
    tymax = (Q[1-r.sign[2]+1][2] - r.origin[2]) * r.inv_direction[2]
    if ((tmin > tymax) || (tymin > tmax))
        return -1.0
    end
    if (tymin > tmin)
        tmin = tymin
    end
    if (tymax < tmax)
        tmax = tymax
    end
    tzmin = (Q[r.sign[3]+1][3] - r.origin[3]) * r.inv_direction[3]
    tzmax = (Q[1-r.sign[3]+1][3] - r.origin[3]) * r.inv_direction[3]

    if ((tmin > tzmax) || (tzmin > tmax))
        return -1.0
    end
    if (tzmin > tmin)
        tmin = tzmin
    end
    if (tzmax < tmax)
        tmax = tzmax
    end
    if (~isfinite(tmax) || ~isfinite(tmin))
        error("An ambiguous ray trajectory detected: ray touches the voxel at its boundary. ")
    end
    return norm((r.origin + r.dir * tmin) - (r.origin + r.dir * tmax))

end

@inline function intersects(a::HyperRectangle{N,Float64}, b::HyperRectangle{N,Float64}) where {N}
    maxA = a.origin + a.widths
    minA = a.origin
    maxB = b.origin + b.widths
    minB = b.origin
    return (all(maxA .>= minB) && all(maxB .> minA))

end

function inittree3d(n)
    empty::Array{Int64,1} = []
    oct = Cell(SVector(-1.0, -1, -1), SVector(2.0, 2, 2), empty)
    for _ = 1:n
        for leaf in allleaves(oct)
            split!(leaf)
        end
    end
    for node in allcells(oct)
        node.data = copy(node.data)
    end
    return oct
end

function printdata(tree::Cell)
    for i = 1:length(tree.children)
        if (tree[i].children !== nothing)
            for j = 1:length(tree[i].children)
                println(i, " ", j, " ", tree[i][j].data)
            end
        else
            println(i, " ", tree[i].data)
        end
    end
end

function inserttotree!(
    box::HyperRectangle{N,Float64},
    tree::Cell,
    index::Int64,
) where {N}
    if (isleaf(tree))
        if (intersects(box, tree.boundary))
            push!(tree.data, index)
        end
    elseif (intersects(box, tree.boundary))
        for i = 1:length(tree.children)
            inserttotree!(box, tree[i], index)
        end
    end
    return nothing
end

function possiblevoxels(r::ray3d, tree::Cell)::Vector{Int64}
    if (isleaf(tree) && intersects(r, tree.boundary) > 0.0)
        return tree.data

    elseif (~(tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblevoxels(r, tree[i]))
        end
        return v

    elseif ((tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblevoxels(r, tree[i]))
        end
        return collect(Set(v))

    else
        return []
    end

end

function voxelray(
    r::ray3d,
    vv::Vector{HyperRectangle{3,Float64}},
    checklist::Vector{Int64},
)
    N = length(checklist)
    indices = Vector{Int64}()
    tlens = Vector{Float64}()
    for i = 1:N
        q = intersects(r, vv[checklist[i]])
        if (q>0.0)
            push!(tlens, q)
            push!(indices, checklist[i])
        end
    end
    return (indices, tlens)
end

@inline function rodrigues(k::T,v::T,theta::Float64)::T where {T}
    kn = k/norm(k)
    return v*cos(theta) + cross(kn,v)*sin(theta) + kn*dot(kn,v)*(1.0-cos(theta))
end

function drawrays(N, theta, fii)
    @assert length(theta) == length(fii)
    Na = length(theta)
    if (~isa(theta, Union{SArray,Array}))
        theta = [theta]
        fii = [fii]
    end
    r = 0.5 * sqrt(12)
    v1 = [1.0, 0, 0]
    v2 = [0, 1.0, 0]
    v3 = [0, 0, 1.0]
    span = range(-r, stop = r, length = N)
    using3D()
    fig = figure()
    ax = fig.add_subplot(111, projection = "3d")
    for a = 1:Na
        x1 = r * sin(theta[a]) * cos(fii[a])
        y1 = r * sin(theta[a]) * sin(fii[a])
        z1 = r * cos(theta[a])
        d = [x1, y1, z1]
        fullangle = acos(dot(d, v1) / (norm(d) * norm(v1)))
        k = cross(v1, [x1, y1, z1])
        mainaxle = rodrigues(k, v1, fullangle)
        axle2 = rodrigues(k, v2, fullangle)
        axle3 = rodrigues(k, v3, fullangle)
        center = @SVector [mainaxle[1] * r, mainaxle[2] * r, mainaxle[3] * r]
        for i = 1:N
            for j = 1:N
                p1 = center + axle2 * span[j] + axle3 * span[i]
                p2 = p1 - 2 * mainaxle * r
                ax.plot3D([p1[1], p2[1]], [p1[2], p2[2]], [p1[3], p2[3]])

            end
        end
    end

end

function constructmatrix(
    tree::Cell,
    vv::Vector{HyperRectangle{3,Float64}},
    N::Int64,
    theta,
    fii,
)
    @assert length(theta) == length(fii)
    v1 = [1.0, 0, 0]
    v2 = [0, 1.0, 0]
    v3 = [0, 0, 1.0]
    r = 0.5 * sqrt(12)

    Na = length(theta)
    if (~isa(theta, Union{SArray,Array}))
        theta = [theta]
        fii = [fii]
    end

    span = range(-r, stop = r, length = N)
    Ntotal = Na * N^2
    Nth = Threads.nthreads()
    rows = Vector{Vector{Int64}}(undef, Nth)
    cols = Vector{Vector{Int64}}(undef, Nth)
    vals = Vector{Vector{Float64}}(undef, Nth)
    for p = 1:Nth
        rows[p] = []
        cols[p] = []
        vals[p] = []
    end

    t = time()
    if (Na > 1)
        pb = ProgressBar(1:Na)
    else
        pb = 1:1
    end
    for a in pb
        x1 = r * sin(theta[a]) * cos(fii[a])
        y1 = r * sin(theta[a]) * sin(fii[a])
        z1 = r * cos(theta[a])
        d = [x1, y1, z1]
        fullangle = acos(dot(d, v1) / (norm(d) * norm(v1)))
        k = cross(v1, [x1, y1, z1])
        mainaxle = rodrigues(k, v1, fullangle)
        axle2 = rodrigues(k, v2, fullangle)
        axle3 = rodrigues(k, v3, fullangle)
        center = @SVector [mainaxle[1] * r, mainaxle[2] * r, mainaxle[3] * r]
        Threads.@threads for i = 1:N
            for j = 1:N
                rayindex = (a - 1) * N^2 + (i - 1) * N + j
                p1 = center + axle2 * span[j] + axle3 * span[i]
                p2 = p1 - 2 * mainaxle * r
                or = ray3d(p1, p2)
                checklist = possiblevoxels(or, tree)
                indices, tlens = voxelray(or, vv, checklist)
                # append!(rows,rayindex*ones(Int64,length(tlens)))
                # append!(cols,indices)
                # append!(vals,tlens)

                append!(rows[Threads.threadid()],rayindex * ones(Int64, length(tlens)))
                append!(cols[Threads.threadid()], indices)
                append!(vals[Threads.threadid()], tlens)

            end
        end
    end
    rows = vcat(rows...)
    cols = vcat(cols...)
    vals = vcat(vals...)
    M = sparse(rows, cols, vals, Ntotal, length(vv))
    println("Matrix constructed in ", time() - t, " seconds.")
    return M
end

function setupvoxels(sizex::Int64,sizey::Int64,sizez::Int64,octsize::Int64)
    oct = inittree3d(octsize)
    l = Array([[1.0, 0.0, 0.0] [0,0,1] [1, 1, 0] [0, 0, 0] [0, 1, 1] [1, 0, 1] [1, 1, 1] [0, 1, 0]]')
    pm = mean(l, dims = 1)
    size = [sizex, sizey, sizez]
    cellsize = 2.0 ./ size
    width = @SVector[cellsize[1], cellsize[2], cellsize[3]]
    for i = 1:3
        l[:, i] =cellsize[i] .* (l[:, i] .- pm[i]) .- cellsize[i] / 2.0 * size[i] .+ 0.5 * cellsize[i]
    end
    voxelvector = Vector{HyperRectangle{3,Float64}}(undef, sizex * sizey * sizez)
    t = time()
    for k = 1:sizez
        for j = 1:sizey
            for i = 1:sizex
                ll = l .+ reshape(cellsize .* [i - 1, j - 1, k - 1], 1, 3)
                origin = @SVector[
                    minimum(ll[:, 1]),
                    minimum(ll[:, 2]),
                    minimum(ll[:, 3]),
                ]
                voxel = HyperRectangle(origin, width)
                index = j + (i - 1) * sizey + (k - 1) * sizey * sizex
                inserttotree!(voxel, oct, index)
                voxelvector[index] = voxel
            end
        end
    end
    println("Voxels initialized in ", time() - t, " seconds.")
    return (oct, voxelvector)
end

function sub2ind(sizes,i::Int64,j::Int64,k::Int64)
    @assert i > 0 && i<= sizes[1]
    @assert j > 0 && j<= sizes[2]
    @assert k > 0 && k<= sizes[3]
    return (k-1)sizes[1]*sizes[2] + sizes[1]*(j-1) + i
end

function test()
    Nx = 100; Ny = 100; Nz = 100; Os = 6; Ns = 100
    (oct,voxelvector)=setupvoxels(Nx,Ny,Nz,Os)
    theta = [0,pi/2]; fii = [0, pi/2]
    M=constructmatrix(oct,voxelvector,Ns,theta,fii)

    l = Array([[1.0, 0.0, 0.0] [0,0,1] [1, 1, 0] [0, 0, 0] [0, 1, 1] [1, 0, 1] [1, 1, 1] [0, 1, 0]]')
    pm = mean(l, dims = 1)
    size = [Nx,Ny,Nz]
    cellsize = 2.0 ./ size
    qx=Vector(range(-1 + cellsize[1]/2.0,stop=1- cellsize[1]/2.0,length=Nx))
    qy=Vector(range(-1+ cellsize[2]/2.0,stop=1- cellsize[2]/2.0,length=Ny))
    qz=Vector(range(-1+ cellsize[3]/2.0,stop=1- cellsize[3]/2.0,length=Nz))

    mgrid = Matlab.meshgrid(qx,qy,qz)
    x = mgrid[1]; y = mgrid[2]; z = mgrid[3]

    object = @. 5.0.*(x^2/1^2 + (y)^2/0.5^2 + z^2/1^2 < 0.25)
    fig = plt.figure()
    imshow(object[:,:,Nz÷2+1])
    objf = object[:]
    sino = M*objf
    sino1 = sino[1:Ns^2]
    sinor1 = reshape(sino1,Ns,Ns)
    fig = plt.figure()
    imshow(sinor1)

    fig = plt.figure()
    sino2 = sino[Ns^2+1:end]
    sinor2 = reshape(sino2,Ns,Ns)
    imshow(sinor2)
    return  nothing

end

test()

max_stars_repo_name: suurj/radon

max_issues_repo_name: suurj/radon

max_forks_repo_name: suurj/radon


text:
abstract type AbstractWineQuality <: Name end

struct RedWineQuality <: AbstractWineQuality end
struct WhiteWineQuality <: AbstractWineQuality end

# mandatory methods
format(::Type{<:AbstractWineQuality}) = TabularData
task(::Type{<:AbstractWineQuality}) = MultiClass
source(::Type{<:AbstractWineQuality}) = "https://archive.ics.uci.edu/ml/datasets/wine+quality"

function downloadlink(::Type{RedWineQuality})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv"
end

function downloadlink(::Type{WhiteWineQuality})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-white.csv"
end

function preprocess(N::Type{<:AbstractWineQuality})
    return path -> preprocess_csv(N, path, :train; istarget = 12, header = true)
end

nclasses(::Type{RedWineQuality}) = 6
nclasses(::Type{WhiteWineQuality}) = 7
classes(::Type{RedWineQuality}) = collect(3:8)
classes(::Type{WhiteWineQuality}) = collect(3:9)
classes_orig(::Type{RedWineQuality}) = collect(3:8)
classes_orig(::Type{WhiteWineQuality}) = collect(3:9)
nattributes(::Type{<:AbstractWineQuality}) = (11, )
ninstances(::Type{RedWineQuality}) = (1599, 0, 0)
ninstances(::Type{WhiteWineQuality}) = (4898, 0, 0)

# optional methods
function checksum(::Type{RedWineQuality})
    return "4a402cf041b025d4566d954c3b9ba8635a3a8a01e039005d97d6a710278cf05e"
end
function checksum(::Type{WhiteWineQuality})
    return "76c3f809815c17c07212622f776311faeb31e87610d52c26d87d6e361b169836"
end

name(::Type{<:AbstractWineQuality}) = "Wine Quality Data Set"
author(::Type{<:AbstractWineQuality}) = ["Paulo Cortez", "A. Cerdeira", "F. Almeida", "T. Matos", "J. Reis"]

function citation(::Type{<:AbstractWineQuality})
    return """
    @article{cortez2009modeling,
        title={Modeling wine preferences by data mining from physicochemical properties},
        author={Cortez, Paulo and Cerdeira, Ant{\'o}nio and Almeida, Fernando and Matos, Telmo and Reis, Jos{\'e}},
        journal={Decision support systems},
        volume={47},
        number={4},
        pages={547--553},
        year={2009},
        publisher={Elsevier}
    }
    """
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
function apply(obj::Transform, arg0::Hyperplane)
    return jcall(obj, "apply", Hyperplane, (Hyperplane,), arg0)
end

function apply(obj::Transform, arg0::Point)
    return jcall(obj, "apply", Point, (Point,), arg0)
end

function apply(obj::Transform, arg0::SubHyperplane, arg1::Hyperplane, arg2::Hyperplane)
    return jcall(obj, "apply", SubHyperplane, (SubHyperplane, Hyperplane, Hyperplane), arg0, arg1, arg2)
end


text:
import Oracle
using Dates

@assert isfile(joinpath(@__DIR__, "..", "test","credentials.jl")) """
Before running tests, create a file `test/credentials.jl` with the content:

username = "your-username"
password = "your-password"
connect_string = "your-connect-string"
auth_mode = Oracle.ORA_MODE_AUTH_DEFAULT # or Oracle.ORA_MODE_AUTH_SYSDBA if user is SYSDBA
"""
include(joinpath(@__DIR__, "..", "test","credentials.jl"))

const NUM_ITERATIONS = 10

function bench()
    conn = Oracle.Connection(username, password, connect_string, auth_mode=auth_mode)

    try

        Oracle.execute(conn, "CREATE TABLE TB_BIND_BY_NAME ( ID NUMBER(15,0) NULL, FLT NUMBER(15,4) NULL, STR VARCHAR(255) NULL, DT DATE NULL)")

        stmt = Oracle.Stmt(conn, "INSERT INTO TB_BIND_BY_NAME ( ID, FLT, STR, DT ) VALUES ( :id, :flt, :str, :dt )")
        @assert stmt.bind_count == 4

        for i in 1:10
            stmt[:id] = 1 + i
            stmt[:flt] = 10.23 + i
            stmt[:str] = "🍕 $i"
            stmt[:dt] = Date(2018,12,31) + Dates.Day(i)
            Oracle.execute(stmt)
        end
        Oracle.commit(conn)

        let
            row_number = 1
            Oracle.query(conn, "SELECT * FROM TB_BIND_BY_NAME") do cursor
                for row in cursor
                    @assert row["ID"] == 1 + row_number
                    @assert row["FLT"] == 10.23 + row_number
                    @assert row["STR"] == "🍕 $row_number"
                    @assert row["DT"] == Date(2018,12,31) + Dates.Day(row_number)

                    row_number += 1
                end
            end
        end

        Oracle.execute(conn, "DELETE FROM TB_BIND_BY_NAME")

        stmt[:id, Int] = missing
        stmt[:flt, Float64] = missing
        stmt[:str, String] = missing
        stmt[:dt, Date] = missing

        Oracle.execute(stmt)
        Oracle.commit(conn)

        let
            row_number = 0
            Oracle.query(conn, "SELECT * FROM TB_BIND_BY_NAME") do cursor
                for row in cursor
                    @assert ismissing(row["ID"])
                    @assert ismissing(row["FLT"])
                    @assert ismissing(row["STR"])
                    @assert ismissing(row["DT"])
                    row_number += 1
                end
            end

            @assert row_number == 1
        end

        Oracle.close(stmt)
        Oracle.execute(conn, "DROP TABLE TB_BIND_BY_NAME")

    finally
        Oracle.close(conn)
    end
end

function main()
    for iter in 1:NUM_ITERATIONS
        bench()
        println("Iteration $iter.")
        GC.gc()
    end
end

main()


text:
using Printf
using Dates
using Base.Threads

const ENDProc =UInt(0)
const MERGEDProc =UInt(1)
const TERMINATEDProc =UInt(2)
#=
 fffTffffffff TTTT fffff TTTTTTT ....
=#
#=
    按照空间 ，即光程从1 到光缆末端进行空间上的合并、过滤，超限处理检查。
    1、过滤： 在合并前，如果 连续告警的区域的空间长度小于 “过滤长度”，则认为该区域为 无效报警区域。
    2、合并： 在两个告警区域之间的 非报警区域 的长度，如果小于 合并间隔，则认为两个告警区域是联通的。
    3、告警的上下限： 在处理完告警的过滤、合并之后，如果连续的告警区域的长度不在 告警上下限之间，则认为属于无效告警。

    输入数据来自于 timeSchmidt的输出，即每段上的报警状态： curAlarmsinSEG,defined in init.jl
    告警参数配置来自于： glbSegAlarmParasDictObj
    输出数据保存在：
        NameTuple ( fiberSegID  ,t0Uint,persistTime)
        # CurAlarmsInSpace = Array{ NamedTupe, 1}
    简化实现：

=#
alarmTuple = NamedTuple{(:fiberSegID ,:t0Uint,:persistTime),Tuple{Int,UInt,UInt} }
historyAlarms = Dict{ UInt,alarmTuple}

function getchn2BegID()
    return UInt( CalChn1_2Obj.CfgInfoObj.chn2SegBegIdx );
    # Chn2Seg2 = UInt(cfgObj.SegNumber) )
    # return 4097)
end
function getSegNUmbers()
    return UInt( CalChn1_2Obj.CfgInfoObj.SegNumber  )    #UInt(8192)
end
mutable struct AlarmPerSeg
    curdtUInt::UInt
    curFrame::UInt
    chan2BegID::UInt
    SegNumbers::UInt
    alarmStatus_1::Array{Bool}
    alarmStatus_2::Array{Bool}
    function AlarmPerSeg()
        curdtUInt=UInt(0)
        chan2BegID = getchn2BegID()
        SegNumbers = getSegNUmbers()
        alarmStatus_1 = Array{Bool}(undef,SegNumbers)
        alarmStatus_1 .= false
        alarmStatus_2 = Array{Bool}(undef,SegNumbers)
        alarmStatus_2 .= false
        new(curdtUInt,0,chan2BegID, SegNumbers,alarmStatus_1,alarmStatus_2 )
    end
end
SegAlarms = nothing
function InitSegAlarms()
    global SegAlarms
    SegAlarms = AlarmPerSeg()
end
#=
    将报警状态拷贝到本段的缓冲中
=#
function l0_handle(  )
    SegAlarms.curdtUInt = curAlarmsinSEG[1].L0time
    SegAlarms.curFrame= curAlarmsinSEG[1].L0Frame

    for i = 1: length( curAlarmsinSEG )
        alr = curAlarmsinSEG[ i ]
        SegAlarms.alarmStatus_1[i] = alr.final_Status
    end
end

function getNextAlarmReg( nextSegAlrs,seg1::UInt,segEnd::UInt)
    s1=UInt(0)
    s2=UInt(0)
    s3=UInt(0)
    if  seg1 >= segEnd || seg1  == 0
        return (s1,s2,s3 )
    end
    for i::UInt = seg1 :segEnd
        isAlr = nextSegAlrs[i]
        # @info isAlr,"-------"
        if isAlr == true
            if s1 == UInt(0)
                s1 = i
                s2 = i
            else
                s2 = i
            end
        end
        if s1 != 0 && isAlr == false # 找到了下一 连续告警区域的末端
            s3 = i
            break
        end
    end
    # @info s1,s2 ,seg1 ,segEnd,"-----------------------------------"
    return (s1,s2,s3) # ( 入侵区域起点，入侵区域终点，下一探测区域的起点:0 表示结束了，其它表示有效的下标)
end

#=
    进行过滤处理
    lastAlrSegReg 包括上一区域的起点、终点 SegID.

    如果本段的下一告警区域与 lastAlrSegReg相隔
        1 超过 合并间隔， 则认为上一区域的告警为一 独立告警区域，返回
            （ Merge = false,)+lastAlrSegReg
        2 未超过 合并间隔， 则认为上一区域的告警与本区域告警合并，返回
        （ Merge = true,)+lastAlrSegReg

=#
function l10_handle( s1::UInt,s2::UInt,s3::UInt,lastSegBeg,lastSegEnd ,mergeL0)
    # 以下进行可否 合并 两个区域的计算。
    # 1、 如果 s0 处的 过滤长度，大于s2-s1,则进行过滤，不能合并。
    mergeL0,FilterL0 = getSegMicroSpaceParaValue( glbSegAlarmParasDictObj , Int(s1) )
        lastmergeL,lastFilterL = getSegMicroSpaceParaValue( glbSegAlarmParasDictObj , Int(lastSegEnd ) )
        intervalSegNum = s1 - lastSegEnd
        # 2、 如果 s0 处的 合并长度，或者 lastSegEnd 的合并长度，大于 两个段间的间隔，则可进行合并，返回合并后的区域。
        #    否则不能合并
        if  intervalSegNum  < mergeL0 || intervalSegNum < lastmergeL
            return ( MERGEDProc ,lastSegBeg , s2, s3 )
        else
            return ( TERMINATEDProc, lastSegBeg, lastSegEnd,s3 )
        end
end
function l1_handle( nextSegAlrs,seg1::UInt,segEnd::UInt,lastAlrSegReg::Tuple )
    lastSegEnd = lastAlrSegReg[2] # 上一区域的终点SegID。
    lastSegBeg =  lastAlrSegReg[1]
    while true
        s1 , s2 , s3 = getNextAlarmReg( nextSegAlrs,seg1,segEnd)
        if s1 == 0  ### 处理结束，到达颗粒的末端了
            return (ENDProc ,lastAlrSegReg[1],lastAlrSegReg[2],s3 )
        end
        mergeL0,FilterL0 = getSegMicroSpaceParaValue( glbSegAlarmParasDictObj , Int(s1) )
        if s2-s1 > FilterL0 # 是一个连续的新报警区域
            if lastSegEnd == 0 || lastSegBeg == 0 # 一个全新的区域
            # 执行至此，已经找到了下一告警区域、或者已经遍历nextSegAlrs，但未找到发生报警区域。
                return ( MERGEDProc ,s1, s2, s3 )
            else
                return l10_handle(s1,s2,s3,lastSegBeg,lastSegEnd,mergeL0)
            end
        else
            seg1 = s3
            continue
        end
    end
    #  end
end

function AlarmRegAna( s_1::UInt ,s_2::UInt )
    maxWidth1 ,minWidth1 = getSegMacroSpaceParaValue(glbSegAlarmParasDictObj, Int( s_1))
    maxWidth2, minWidth2 = getSegMacroSpaceParaValue(glbSegAlarmParasDictObj, Int( s_2))
    minwidth = floor( UInt,(minWidth1+minWidth2) /2 )
    maxwidth = floor( UInt,(maxWidth1+maxWidth2) /2 )
    width = UInt(s_2 -s_1 +1)
    # @warn "$minwidth , $maxwidth $width"
    if width >= minwidth &&  width <= maxwidth
        if true #WindID.isinWindRng(s_1,s_2  ) == false
            AlarmMerge.mergeProc(s_1,s_2,SegAlarms.curdtUInt) ## 进行大尺度空间、时间合并，产生上报的报警
            EventL1.ReportEventL1(s_1,s_2, SegAlarms.curdtUInt,SegAlarms.curFrame ,"Windoff")
        else
            EventL1.ReportEventL1(s_1,s_2, SegAlarms.curdtUInt,SegAlarms.curFrame ,"Windon")
            @warn "filterted by wind : $s_1,$s_2"
        end
    else
        if width < minwidth
            #EventL1.ReportEventL1(s_1,s_2, SegAlarms.curdtUInt,SegAlarms.curFrame ,"RNGNarrow")
        elseif width > maxwidth
            EventL1.ReportEventL1(s_1,s_2, SegAlarms.curdtUInt,SegAlarms.curFrame ,"RNGWide")
        end

    end
end
#=

=#
function l2_handle(segID0::UInt,SegID1::UInt )
    nextSegAlrs = SegAlarms.alarmStatus_1
    seg1 = UInt( segID0 )
    segEnd = UInt( SegID1 )
    # @info "$segEnd,-------------"
    # nextSegAlrs = alarms
    s_1 ,s_2 ,s_3  = getNextAlarmReg( nextSegAlrs,seg1,segEnd)
    s1 = s_1
    s2 = s_2
    seg1 =s_3
    if s1 == 0 && s2 == 0
        return  # 未见报警
    end
    doesMerged = TERMINATEDProc
    while true #开始处理过滤、合并
        # if s3 != 0 ## 到了末端，不需要 继续
        # @info "$doesMerged ,$s_1 ,$s_2 ,$s_3 "
        doesMerged ,s_1 ,s_2 ,s_3 = l1_handle( nextSegAlrs,seg1,segEnd,(s1,s2) )
        # end
        if doesMerged == ENDProc
            if  s_1 != 0
                AlarmRegAna( UInt(s_1) ,UInt(s_2) )
                break
            elseif s_1 == 0
                break
            end
        elseif doesMerged == TERMINATEDProc  && s_1 != 0 ## 一个独立的报警区域已经识别到位，保存在s_1 S_2中
            AlarmRegAna( UInt(s_1) ,UInt(s_2) )
            # s1 = 0
            # s2 = 0
            if s_3 == 0
                break
            end
            seg1 =s_3
            s_1 ,s_2 ,s_3 = getNextAlarmReg( nextSegAlrs,seg1,segEnd)
            s1 = s_1
            s2 = s_2
            seg1 =s_3
            # @info doesMerged ,s_1 ,s_2 ,s_3
        elseif doesMerged == MERGEDProc ## 已经合并了一个新区域，需要继续合并下一个区域
            s1 = s_1
            s2 = s_2
            seg1 = s_3
            # @info doesMerged ,s_1 ,s_2 ,s_3
        end
    end
    # @info "end space "
    return
end

function updateSpaceStatus( )
    l0_handle( )

    segID0 = UInt(1); segID1 = SegAlarms.chan2BegID -1 ;

    # @info  "$segID0,$segID1"

    segID2 = segID1 +1 ; segID3 = SegAlarms.SegNumbers ;

    a= Base.Threads.@spawn   l2_handle(segID0,segID1 )
    b= Base.Threads.@spawn  l2_handle(segID2,segID3 )
    wait(a)
    wait(b)

end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
## Import packages

using RadiativeTransfer, RadiativeTransfer.vSmartMOM
using Statistics, Plots

include("6SV1_R_trues.jl")
R_modeled_all = zeros(6, 3, 3, 16);
R_deltas_all = zeros(6, 3, 3, 16);
parameters = parameters_from_yaml("/home/rjeyaram/RadiativeTransfer/test/6SV1_1.yaml");

## Function to perform test, given conditions

function test_against_6SV1(case_i, azs, szas, λ, τ, ρ)
    for sza_i in 1:length(szas)
        for az_i in 1:length(azs)
            R_true = R_trues[case_i]
            R_modeled = view(R_modeled_all, case_i, :, :, :)
            R_deltas = view(R_deltas_all, case_i, :, :, :)
            parameters.spec_bands = [1e7/λ (1e7/λ + 1)]
            parameters.vaz = repeat([azs[az_i]], 16)
            parameters.sza = szas[sza_i]
            parameters.brdf = [RadiativeTransfer.vSmartMOM.LambertianSurfaceScalar(ρ * π)]
            model = model_from_parameters(parameters);
            model.τ_rayl[1] .= τ
            R_modeled[sza_i, az_i, :] = vSmartMOM.rt_run(model, i_band=1)[:,1,1] / model.quad_points.μ₀
            R_deltas[sza_i, az_i, :] = abs.(R_true[sza_i][az_i] - R_modeled[sza_i, az_i, :]) ./ R_true[sza_i][az_i]
        end
    end

    println("Maximum discrepancy: ", maximum(R_deltas_all[case_i,:,:,:])*100, "%")
    println("Minimum discrepancy: ", minimum(R_deltas_all[case_i,:,:,:])*100, "%")
    println("Average discrepancy: ", mean(R_deltas_all[case_i,:,:,:])*100, "%")
    println("Std dev discrepancy: ", std(R_deltas_all[case_i,:,:,:])*100, "%")
end

## Case 1: τ = 0.10, (λ = 530 nm), ρ = 0.0

test_against_6SV1(1,                               # test case no.
                  [180, 90, 0],                    # AZ
                  [23.0739, 53.1301, 78.4630],     # SZA
                  530,                             # λ
                  0.1,                             # τ
                  0.0)                             # ρ

## Case 2: τ = 0.10, (λ = 530 nm), ρ = 0.25

test_against_6SV1(2,                               # test case no.
                  [180, 90, 0],                    # AZ
                  [0.0001, 36.8699, 66.4218],      # SZA
                  530,                             # λ
                  0.1,                             # τ
                  0.25)                            # ρ

## Case 3: τ = 0.25, (λ = 440 nm), ρ = 0.0

test_against_6SV1(3,                               # test case no.
                  [180, 90, 0],                    # AZ
                  [0.0001, 36.8699, 66.4218],      # SZA
                  440,                             # λ
                  0.25,                            # τ
                  0.0)                             # ρ

## Case 4: τ = 0.25, (λ = 440 nm), ρ = 0.25

test_against_6SV1(4,                               # test case no.
                  [180, 90, 0],                    # AZ
                  [23.0739, 53.1301, 78.4630],     # SZA
                  440,                             # λ
                  0.25,                            # τ
                  0.25)                            # ρ

## Case 5: τ = 0.50, (λ = 360 nm), ρ = 0.0

test_against_6SV1(5,                               # test case no.
                  [180, 90, 0],                    # AZ
                  [23.0739, 53.1301, 78.4630],     # SZA
                  360,                             # λ
                  0.50,                            # τ
                  0.0)                             # ρ

## Case 6: τ = 0.50, (λ = 360 nm), ρ = 0.25

test_against_6SV1(6,                               # test case no.
                  [180, 90, 0],                    # AZ
                  [0.0001, 36.8699, 66.4218],      # SZA
                  360,                             # λ
                  0.50,                            # τ
                  0.25)                            # ρ


text:
#=
common code for examples
=#

abstract type ExampleInstance{T <: Real} end

# this is a workaround for randn's lack of support for BigFloat
Random.randn(R::Type{BigFloat}, dims::Vararg{Int, N} where N) = R.(randn(dims...))
Random.randn(R::Type{Complex{BigFloat}}, dims::Vararg{Int, N} where N) =
    R.(randn(ComplexF64, dims...))

# helper for calculating solution violations
relative_residual(residual::Vector{T}, constant::Vector{T}) where {T <: Real} =
    norm(residual, Inf) / (1 + norm(constant, Inf))

# calculate violations for Hypatia certificate equalities
function certificate_violations(
    status::Solvers.Status,
    model::Models.Model{T},
    x::Vector{T},
    y::Vector{T},
    z::Vector{T},
    s::Vector{T},
    ) where {T <: Real}
    (c, A, b, G, h, obj_offset) =
        (model.c, model.A, model.b, model.G, model.h, model.obj_offset)

    if status == Solvers.PrimalInfeasible
        x_res = G' * z + A' * y
        x_viol = relative_residual(x_res, c)
        y_viol = T(NaN)
        z_viol = T(NaN)
    elseif status == Solvers.DualInfeasible
        y_res = A * x
        z_res = G * x + s
        y_viol = relative_residual(y_res, b)
        z_viol = relative_residual(z_res, h)
        x_viol = T(NaN)
    # TODO elseif status == Solvers.IllPosed (primal vs dual ill-posed)
    else
        if status != Solvers.Optimal
            println("status $status not handled, but computing optimality " *
                "certificate violations anyway")
        end
        x_res = G' * z + A' * y + c
        y_res = A * x - b
        z_res = G * x + s - h
        x_viol = relative_residual(x_res, c)
        y_viol = relative_residual(y_res, b)
        z_viol = relative_residual(z_res, h)
    end

    return (x_viol, y_viol, z_viol)
end

# return solve information and certificate violations
function process_result(
    model::Models.Model{T},
    solver::Solvers.Solver{T},
    ) where {T <: Real}
    status = Solvers.get_status(solver)
    solve_time = Solvers.get_solve_time(solver)
    iters = Solvers.get_num_iters(solver)
    primal_obj = Solvers.get_primal_obj(solver)
    dual_obj = Solvers.get_dual_obj(solver)
    rel_obj_diff = (primal_obj - dual_obj) / (1 + abs(dual_obj))

    z = Solvers.get_z(solver)
    s = Solvers.get_s(solver)
    x = Solvers.get_x(solver)
    y = Solvers.get_y(solver)
    compl = dot(s, z)
    (x_viol, y_viol, z_viol) = certificate_violations(status, model, x, y, z, s)

    flush(stdout); flush(stderr)

    solve_stats = (;
        status, solve_time, iters, primal_obj, dual_obj,
        rel_obj_diff, compl, x_viol, y_viol, z_viol,
        :time_rescale => solver.time_rescale,
        :time_initx => solver.time_initx,
        :time_inity => solver.time_inity,
        :time_unproc => solver.time_unproc,
        :time_loadsys => solver.time_loadsys,
        :time_upsys => solver.time_upsys,
        :time_upfact => solver.time_upfact,
        :time_uprhs => solver.time_uprhs,
        :time_getdir => solver.time_getdir,
        :time_search => solver.time_search,
        )
    solution = (; x, y, z, s)

    return (solve_stats, solution)
end

get_model_stats(model::Models.Model) = (;
    :n => model.n, :p => model.p, :q => model.q, :nu => model.nu,
    :cone_types => [string(nameof(c)) for c in unique(typeof.(model.cones))],
    :num_cones => length(model.cones),
    :max_q => maximum(Cones.dimension(c) for c in model.cones),
    )