text:
using GeometryTypes, StaticArrays, ModernGL, Visualize
import GLAbstraction, GLWindow, ColorVectorSpace
import Transpiler: mix, smoothstep, gli
using Visualize: orthographicprojection

function aastep{T}(threshold1::T, value)
    smoothstep(threshold1 - T(0.05), threshold1 + T(0.05), value)
end

type Uniforms{F}
    projection::Mat4f0
    strokecolor::Vec4f0
    glowcolor::Vec4f0
    distance_func::F
end

immutable Vertex{N, T}
    position::Vec{N, T}
    uvrect::Vec4f0
    color::Vec4f0
    scale::Vec2f0
end

immutable Vertex2Geom
    uvrect::Vec4f0
    color::Vec4f0
    rect::Vec4f0
end

getuv(x) = x.uvrect
getcolor(x) = x.color
getstrokecolor(x) = x.strokecolor
getglowcolor(x) = x.glowcolor
getscale(x) = x.scale
getposition(x) = x.position

function vertex_main(vertex, uniforms)
    p = getposition(vertex)
    scale = getscale(vertex)

    return Vertex2Geom(
        getuv(vertex),
        getcolor(vertex),
        Vec4f0(p[1], p[2], scale[1], scale[2])
    )
end

"""
Emits a vertex with
"""
function emit_vertex(emit!, vertex, uv, offsetted_uv, arg, pos, uniforms)
    datapoint = uniforms.projection * Vec4f0(pos[1], pos[2], 0, 1)
    final_position = uniforms.projection * Vec4f0(vertex[1], vertex[2], 0, 0)
    frag_out = (uv, offsetted_uv, arg.color)
    emit!(datapoint .+ final_position, frag_out)
    return
end

function geometry_main(emit!, vertex_out, uniforms)
    # get arguments from first face
    # (there is only one in there anywas, since primitive type is point)
    # (position, vertex_out)
    arg = vertex_out[1]
    # emit quad as triangle strip
    # v3. ____ . v4
    #    |\   |
    #    | \  |
    #    |  \ |
    #    |___\|
    # v1*      * v2
    pos_scale = arg.rect
    pos = pos_scale[Vec(1, 2)]
    scale = pos_scale[Vec(3, 4)]
    quad = Vec4f0(0f0, 0f0, scale[1], scale[2])
    uv = arg.uvrect
    uvnormed = Vec4f0(-0.5f0, -0.5f0, 0.5f0, 0.5f0)
    emit_vertex(emit!, quad[Vec(1, 2)], uvnormed[Vec(1, 4)], uv[Vec(1, 4)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(1, 4)], uvnormed[Vec(1, 2)], uv[Vec(1, 2)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(3, 2)], uvnormed[Vec(3, 4)], uv[Vec(3, 4)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(3, 4)], uvnormed[Vec(3, 2)], uv[Vec(3, 2)], arg, pos, uniforms)
    return
end
function fragment_main(geom_out, uniforms)
    uv = geom_out[1]; uv_offset = geom_out[2]; color = geom_out[3];
    signed_distance = uniforms.distance_func(uv)
    inside = aastep(0f0, signed_distance)
    bg = Vec4f0(1f0, 1f0, 1f0, 0f0)
    (mix(bg, color, inside),)
end
circle{T}(uv::Vec{2, T}) = T(0.5) - norm(uv)

resolution = (500, 500)
w = GLWindow.create_glcontext(resolution = resolution)

proj = orthographicprojection(SimpleRectangle(0, 0, resolution...), -10_000f0, 10_000f0)
uniforms = Uniforms(
    proj,
    Vec4f0(1, 0, 0, 1),
    Vec4f0(1, 0, 1, 1),
    circle
)

N = 20
vertices = [Vertex(
    Vec2f0((sin(2pi * (i / N)) , cos(2pi * (i / N))) .* 200f0) .+ 250f0,
    Vec4f0(0, 0, 0, 0), Vec4f0(1, i/N, 0, 1), Vec2f0(40, 40)
) for i = 1:N]

vbo = VertexArray(vertices)

uniform_buff = UniformBuffer(uniforms);

draw_particles = GLRasterizer(
    vbo, (uniform_buff,),
    vertex_main, fragment_main;
    geometry_main = geometry_main
)

glDisable(GL_DEPTH_TEST)
glClearColor(1, 1, 1, 0)
GLAbstraction.enabletransparency()

while isopen(w)
    GLWindow.poll_glfw()
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    draw_particles(vbo, (uniform_buff,))
    GLWindow.swapbuffers(w)
end
GLFW.DestroyWindow(w)


text:
using Mass2, Base.Test

# generate a spectra with multiple peaks for calibration testing
dist_energies = [4000,4700,5500,6000,6800,7500,8200,9600,10000]
dist_filt_values = 5*dist_energies.^0.8
dist_use_inds = [1,3,6,7,8,9]
dist_use_filt_values = dist_filt_values[dist_use_inds]
dist_use_energies = dist_energies[dist_use_inds]
filt_value_distribution = MixtureModel([Normal(mu,1) for mu in dist_filt_values])

steps = AbstractStep[]

push!(steps, MockPulsesStep(TupacLikeTwoExponentialPulseGenerator(Int, filt_value_distribution), 100, 5000, [:pulse,:rowstamp]))
push!(steps, PerPulseStep(compute_summary, [:pulse, :pre_samples, :frame_time],
	[:pretrig_mean, :pretrig_rms, :pulse_average, :pulse_rms, :rise_time, :postpeak_deriv, :peak_index, :peak_value, :min_value]))
push!(steps, PerPulseStep(selectfromcriteria, [:pretrig_rms, :pretrig_rms_criteria, :peak_index, :peak_index_criteria, :postpeak_deriv, :postpeak_deriv_criteria], [:selection_good]))
push!(steps, PerPulseStep(filter1lag, [:pulse, :whitenoise_filter], [:filt_value]))
push!(steps, HistogramStep(update_histogram!, [:filt_value_hist, :selection_good, :filt_value], [:filt_value_hist]))
push!(steps, ThresholdStep(calibrate_nofit, [:filt_value_dc_hist,:known_energies],[:calibration],:filt_value_dc_hist, counted, 3000, true))
push!(steps, ThresholdStep(compute_whitenoise_filter, [:pulse, :selection_good], [:whitenoise_filter], :selection_good, sum, 100, true))
push!(steps, ThresholdStep(calc_dc, [:pretrig_mean, :filt_value, :selection_good], [:ptm_dc],:filt_value_hist, counted, 3000, true))
push!(steps, PerPulseStep(applycorrection, [:ptm_dc, :pretrig_mean, :filt_value], [:filt_value_dc]))
push!(steps, HistogramStep(update_histogram!, [:filt_value_dc_hist, :selection_good, :filt_value_dc], [:filt_value_dc_hist]))
push!(steps, PerPulseStep(apply_calibration, [:calibration, :filt_value_dc], [:energy]) )
push!(steps, HistogramStep(update_histogram!, [:energy_hist, :selection_good, :energy], [:energy_hist]))

#push!(steps, ToJLDStep([:filt_value,:pretrig_rms, :energy:
# push!(steps, FreeMemoryStep())

mc=MassChannel()

mc[:pretrig_mean] = RunningVector(Float64)
mc[:pretrig_rms] = RunningVector(Float64)
mc[:pulse_average] = RunningVector(Float64)
mc[:pulse_rms] = RunningVector(Float64)
mc[:rise_time] = RunningVector(Float64)
mc[:postpeak_deriv] = RunningVector(Float64)
mc[:peak_index] = RunningVector(UInt16)
mc[:peak_value] = RunningVector(UInt16)
mc[:min_value] = RunningVector(UInt16)
mc[:filt_value] = RunningVector(Float64)
mc[:filt_value_dc] = RunningVector(Float64)
mc[:pretrig_rms] = RunningVector(Float64)
mc[:selection_good] = RunningSumBitVector()
mc[:energy] = RunningVector(Float64)
mc[:filt_value_hist] = Histogram(0:1:20000)
mc[:filt_value_dc_hist] = Histogram(0:1:20000)
mc[:energy_hist] = Histogram(0:1:20000)
mc[:pulse] = RunningVector(Vector{Int})
mc[:rowstamp] = RunningVector(Int)
mc[:pre_samples] = 100
mc[:frame_time] = 1/100000
mc[:peak_index_criteria] = (170,200)
mc[:pretrig_rms_criteria] = (0.0,10.)
mc[:postpeak_deriv_criteria] = (0.0,50.0)
mc[:known_energies] = dist_use_energies

# metadata
mc[:hdf5_filename] = mktemp()[1]
setsteps!(mc, steps)
preptasks!(mc)

savegraph("graph",graph(steps,mc))

isfile(mc[:hdf5_filename]) && rm(mc[:hdf5_filename])
schedule(mc)
sleep(2) # make sure steps have time do some work before planning to end, test fails without this
# may need better system for this
plantoend(mc)
wait(mc.task.value)

eh = mc[:energy_hist]
@test counted(eh) == sum(counts(eh))

# peakinds returns from lowest counts to highest counts
# lets take the largest peaks and match them up to their known energies from distenergies
peakinds = findpeaks(eh, fwhm=15, n=length(dist_energies))
found_energies = sort(bin_centers(eh)[peakinds])
fvh = mc[:filt_value_hist]
peakinds_filt_value = findpeaks(fvh, fwhm=15, n=length(dist_energies))
found_filt_values = sort(bin_centers(fvh)[peakinds])
cal = mc[:calibration]

# filt value agrees with manual calculation
fvmanual = [dot(p, mc[:whitenoise_filter]) for p in mc[:pulse][1:end]]
@test fvmanual == mc[:filt_value][1:end]

@test all(abs(found_energies-dist_energies).<25) # check that the calibration has the correct peak assignment

ext: jl


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for transformations.jl

@testset "Transformation (testset 1)" begin
    sI    = SMatrix{4, 4, Float32}(I(4))
    m     = Diagonal([1f0, 2f0, 3f0, 4f0])
    sm    = SMatrix{4, 4, Float32}(m)
    invm  = inv(m)
    sinvm = SMatrix{4, 4, Float32}(invm)

    @testset "constructor" begin
        def_t = Transformation()
        @test def_t.m == def_t.invm == sI
        @test isconsistent(def_t)

        t1 = Transformation(m)
        @test t1.m == sm && t1.invm == invm
        @test isconsistent(t1)

        t2 = Transformation(sm)
        @test t2.m == m && t2.invm == invm
        @test isconsistent(t2)

        t3 = Transformation(m, invm)
        @test t3.m == m && t3.invm == invm
        @test isconsistent(t3)

        t4 = Transformation(sm, invm)
        @test t4.m == m && t4.invm == invm
        @test isconsistent(t4)

        t5 = Transformation(m, sinvm)
        @test t5.m == m && t5.invm == invm
        @test isconsistent(t5)

        t6 = Transformation(sm, sinvm)
        @test t6.m == m && t6.invm == invm
        @test isconsistent(t6)

        @test t1 == t2 == t3 == t4 == t5 == t6
    end

    t1 = translation(1f0, 2f0, 3f0)
    t2 = rotationZ(π/4f0)
    vv = [1f0, 2f0, 3f0]
    v = Vec(vv...)
    n = Normal(vv...)
    p = Point(vv...)

    @testset "operations" begin
        prod = t1 * t2

        @test isconsistent(prod)

        @test prod.m    == t1.m * t2.m
        @test prod.invm == t2.invm * t1.invm

        @test prod * v == Vec(prod.m[1:3,1:3] * v)
        @test prod * n == Normal(transpose(prod.invm[1:3,1:3]) * v)
        tmp = prod.m * [p.v..., 1.]
        @test prod * p == Point(@view((tmp/tmp[end])[1:3]))

        @test t1 ≈ t1
    end

    @testset "methods" begin
        t = t1 * t2
        invt = inv(t)
        @test t.m == invt.invm && t.invm == invt.m

        θ = π/4f0
        rot_mat = Dict(
            :X => @SMatrix(Float32[   1       0      0     0;
                                      0     cos(θ) -sin(θ) 0;
                                      0     sin(θ)  cos(θ) 0;
                                      0       0      0     1]),
            :Y => @SMatrix(Float32[ cos(θ)    0    sin(θ)  0;
                                      0       1      0     0;
                                    -sin(θ)   0    cos(θ)  0;
                                      0       0      0     1]),
            :Z => @SMatrix(Float32[ cos(θ) -sin(θ)   0     0;
                                    sin(θ) cos(θ)    0     0;
                                      0      0       1     0;
                                      0      0       0     1])
        )

        @test rotationX(θ) ≈ Transformation(rot_mat[:X], transpose(rot_mat[:X]))
        @test rotationY(θ) ≈ Transformation(rot_mat[:Y], transpose(rot_mat[:Y]))
        @test rotationZ(θ) ≈ Transformation(rot_mat[:Z], transpose(rot_mat[:Z]))

        id = Diagonal(ones(eltype(v), 4)) |> MMatrix{4, 4, Float32}
        id⁻¹ = copy(id)
        id[1:3, end]   =  v
        id⁻¹[1:3, end] = -v
        @test translation(v)    ≈ Transformation(id, id⁻¹)
        @test translation(v...) ≈ Transformation(id, id⁻¹)

        @test scaling(v) ≈ Transformation(Diagonal([v...,true]), Diagonal(true ./ [v..., true]))
        @test scaling(v...) ≈ Transformation(Diagonal([v...,true]), Diagonal(true ./ [v..., true]))
        @test scaling(5) ≈ Transformation(Diagonal([5,5,5,1]), Diagonal(inv.([5,5,5,1])))
    end
end

@testset "Transformation (testset 2)" begin
    @testset "constructor" begin
        m1 = Transformation([1f0 2f0 3f0 4f0;
                             5f0 6f0 7f0 8f0;
                             9f0 9f0 8f0 7f0;
                             6f0 5f0 4f0 1f0],
                            [ -3.75f0   2.75f0 -1f0    0f0;
                              4.375f0 -3.875f0  2f0 -0.5f0;
                                0.5f0    0.5f0 -1f0    1f0;
                             -1.375f0  0.875f0  0f0 -0.5f0])
        @test isconsistent(m1)

        m2 = Transformation(m1.m, m1.invm)
        @test m1 ≈ m2

        @test !isconsistent(Transformation([1f0 2f0 3f0 4f0;
                                            5f0 6f0 8f0 8f0;
                                            9f0 9f0 8f0 7f0;
                                            6f0 5f0 4f0 1f0],
                                           [1f0 2f0 3f0 4f0;
                                            5f0 6f0 8f0 8f0;
                                            9f0 9f0 8f0 7f0;
                                            6f0 5f0 4f0 1f0]))

        m3 = Transformation([1f0 2f0 3f0 4f0;
                             5f0 6f0 8f0 8f0;
                             9f0 9f0 8f0 7f0;
                             6f0 5f0 4f0 1f0])
        @test !(m1 ≈ m3)
    end

    @testset "multiplication" begin
        m1 = Transformation([1f0 2f0 3f0 4f0;
                             5f0 6f0 7f0 8f0;
                             9f0 9f0 8f0 7f0;
                             6f0 5f0 4f0 1f0],
                            [ -3.75f0   2.75f0 -1f0    0f0;
                              4.375f0 -3.875f0  2f0 -0.5f0;
                                0.5f0    0.5f0 -1f0    1f0;
                             -1.375f0  0.875f0  0f0 -0.5f0])
        @test isconsistent(m1)

        m2 = Transformation([3f0 5f0 2f0 4f0;
                             4f0 1f0 0f0 5f0;
                             6f0 3f0 2f0 0f0;
                             1f0 4f0 2f0 1f0],
                            [  0.4f0 -0.2f0  0.2f0 -0.6f0;
                               2.9f0 -1.7f0  0.2f0 -3.1f0;
                             -5.55f0 3.15f0 -0.4f0 6.45f0;
                              -0.9f0  0.7f0 -0.2f0  1.1f0])
        @test isconsistent(m2)

        expected = Transformation([ 33f0  32f0 16f0 18f0;
                                    89f0  84f0 40f0 58f0;
                                   118f0 106f0 48f0 88f0;
                                    63f0  51f0 22f0 50f0],
                                  [ -1.45f0    1.45f0  -1.0f0  0.6f0;
                                   -13.95f0   11.95f0  -6.5f0  2.6f0;
                                   25.525f0 -22.025f0 12.25f0 -5.2f0;
                                    4.825f0  -4.325f0   2.5f0 -1.1f0])
        @test isconsistent(expected)

        @test expected ≈ (m1 * m2)
    end

    @testset "vec and point multiplication" begin
        m = Transformation([1f0 2f0 3f0 4f0;
                            5f0 6f0 7f0 8f0;
                            9f0 9f0 8f0 7f0;
                            0f0 0f0 0f0 1f0],
                           [-3.75f0  2.75f0 -1f0  0f0;
                             5.75f0 -4.75f0  2f0  1f0;
                            -2.25f0  2.25f0 -1f0 -2f0;
                                0f0     0f0  0f0  1f0])
        @test isconsistent(m)

        expected_v = Vec(14f0, 38f0, 51f0)
        @test expected_v ≈ (m * Vec(1f0, 2f0, 3f0))

        expected_p = Point(18f0, 46f0, 58f0)
        @test expected_p ≈ (m * Point(1f0, 2f0, 3f0))

        expected_n = Normal(-8.75f0, 7.75f0, -3.0f0)
        @test expected_n ≈ (m * Normal(3.0f0, 2.0f0, 4.0f0))
    end

    @testset "inv" begin
        m1 = Transformation([1f0 2f0 3f0 4f0;
                             5f0 6f0 7f0 8f0;
                             9f0 9f0 8f0 7f0;
                             6f0 5f0 4f0 1f0],
                            [ -3.75f0   2.75f0 -1f0    0f0;
                              4.375f0 -3.875f0  2f0 -0.5f0;
                                0.5f0    0.5f0 -1f0    1f0;
                             -1.375f0  0.875f0  0f0 -0.5f0])

        m2 = inv(m1)
        @test isconsistent(m1)

        prod = m1 * m2
        @test isconsistent(prod)
        @test prod ≈ Transformation()
    end

    @testset "translations" begin
        tr1 = translation(Vec(1f0, 2f0, 3f0))
        @test isconsistent(tr1)

        tr1_2 = translation(Vec(1f0, 2f0, 3f0)...)
        @test isconsistent(tr1_2)
        @test tr1 ≈ tr1_2

        tr2 = translation(Vec(4f0, 6f0, 8f0))
        @test isconsistent(tr1)

        prod = tr1 * tr2
        @test isconsistent(prod)

        expected = translation(Vec(5f0, 8f0, 11f0))
        @test prod ≈ expected
    end

    @testset "rotations" begin
        @test isconsistent(rotationX(0.1f0))
        @test isconsistent(rotationY(0.1f0))
        @test isconsistent(rotationZ(0.1f0))

        @test (rotationZ(π/2) * VEC_X) ≈ VEC_Y
        @test (rotationX(π/2) * VEC_Y) ≈ VEC_Z
        @test (rotationY(π/2) * VEC_Z) ≈ VEC_X
    end

    @testset "scalings" begin
        tr1 = scaling(Vec(2f0, 5f0, 10f0))
        @test isconsistent(tr1)

        tr2 = scaling(Vec(3f0, 2f0, 4f0))
        @test isconsistent(tr2)

        expected = scaling(Vec(6f0, 10f0, 40f0))
        @test expected ≈ (tr1 * tr2)
    end
end


text:
##
##  Java Information Dynamics Toolkit (JIDT)
##  Copyright (C) 2012, Joseph T. Lizier
##
##  This program is free software: you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation, either version 3 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program.  If not, see <http://www.gnu.org/licenses/>.
##

# = Example 4 - Transfer entropy on continuous data using Kraskov estimators =

# Simple transfer entropy (TE) calculation on continuous-valued data using the Kraskov-estimator TE calculator.

# Import the JavaCall package:
using JavaCall;

# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar";
# Start the JVM supplying classpath and heap size
#  (increase memory here if you get crashes due to not enough space)
JavaCall.init(["-Djava.class.path=$(jarLocation)", "-Xmx128M"]);

# Generate some random normalised data.
numObservations = 1000;
covariance=0.4;
sourceArray=randn(numObservations);
destArray = [0, covariance*sourceArray[1:numObservations-1] + (1-covariance)*randn(numObservations - 1)];
sourceArray2=randn(numObservations); # Uncorrelated source

# Create a TE calculator and run it:
teClass = @jimport infodynamics.measures.continuous.kraskov.TransferEntropyCalculatorKraskov;
teCalc=teClass(());
jcall(teCalc, "setProperty", Void, (JString, JString), "k", "4"); # Use Kraskov parameter K=4 for 4 nearest points
jcall(teCalc, "initialise", Void, (jint, ), 1); # Use history length 1 (Schreiber k=1)
# Perform calculation with correlated source:
jcall(teCalc, "setObservations", Void, (Array{jdouble,1}, Array{jdouble,1}),
				     sourceArray, destArray);
result = jcall(teCalc, "computeAverageLocalOfObservations", jdouble, ());
# Note that the calculation is a random variable (because the generated
#  data is a set of random variables) - the result will be of the order
#  of what we expect, but not exactly equal to it; in fact, there will
#  be a large variance around it.
@printf("TE result %.4f nats; expected to be close to %.4f nats for these correlated Gaussians\n",
    result, log(1/(1-covariance^2)));

# Perform calculation with uncorrelated source:
jcall(teCalc, "initialise", Void, ()); # Initialise leaving the parameters the same
jcall(teCalc, "setObservations", Void, (Array{jdouble,1}, Array{jdouble,1}),
				     sourceArray2, destArray);
result2 = jcall(teCalc, "computeAverageLocalOfObservations", jdouble, ());
@printf("TE result %.4f nats; expected to be close to 0 nats for these uncorrelated Gaussians\n", result2);

# We can also compute the local TE values for the time-series samples here:
#  (See more about utility of local TE in the CA demos)
localTE = jcall(teCalc, "computeLocalOfPreviousObservations", Array{jdouble,1}, ());
@printf("Notice that the mean of locals, %.4f nats, equals the previous result\n",
	sum(localTE)/(numObservations-1));

max_issues_repo_name: jmmccracken/ECA


text:
"""

```
coor(atoms; selection)
```

Returns the coordinates of the atoms. The input may be one atom (type `Atom`), a vector of atoms, or a `Residue`.
The coordinates are returned as a vector of static vectors (from `StaticArrays`), more specifically
as a `Vector{SVector{3,Float64}}`.

### Examples

```julia-repl
julia> using PDBTools, StaticArrays

julia> protein = wget("1LBD");

julia> coor(protein[1])
3-element SVector{3, Float64} with indices SOneTo(3):
 45.228
 84.358
 70.638

julia> coor(protein[1],as=SVector{3,Float32})
3-element SVector{3, Float32} with indices SOneTo(3):
 45.228
 84.358
 70.638

julia> coor(protein, "index <= 2")
2-element Vector{SVector{3, Float64}}:
 [45.228, 84.358, 70.638]
 [46.08, 83.165, 70.327]

julia> coor(protein, only = at -> at.resname == "ALA")
110-element Vector{SVector{3, Float64}}:
 [43.94, 81.982, 70.474]
 [43.02, 80.825, 70.455]
 [41.996, 80.878, 69.34]
 ⋮
 [-17.866, 84.088, 51.741]
 [-18.496, 83.942, 52.777]
 [-15.888, 82.583, 51.706]

julia> residues = collect(eachresidue(protein));

julia> coor(residues[1])
6-element Vector{SVector{3, Float64}}:
 [45.228, 84.358, 70.638]
 [46.08, 83.165, 70.327]
 [45.257, 81.872, 70.236]
 [45.823, 80.796, 69.974]
 [47.147, 82.98, 71.413]
 [46.541, 82.639, 72.662]

```

"""
coor(atom::Atom) = SVector{3,Float64}(atom.x, atom.y, atom.z)

function coor(atoms::AbstractVector{Atom}, selection::String)
    query = parse_query(selection)
    return coor(atoms, only = atom -> apply_query(query, atom))
end

function coor(atoms::AbstractVector{Atom}; only = all)
    n = 0
    for atom in atoms
        !only(atom) && continue
        n = n + 1
    end
    x = fill(zero(SVector{3,Float64}),n)
    i = 0
    for atom in atoms
        !only(atom) && continue
        i += 1
        x[i] = coor(atom)
    end
    return x
end

#
# Coordinates of the atoms of a residue/molecule
#
coor(residue::Residue; only = all) = coor(residue.atoms[residue.range], only = only)
coor(residue::Residue, selection::String) = coor(residue.atoms[residue.range], selection)

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
#
# Module to load the files to compute energies of PDB files
#

module PDBEnergy

  include("./src/Atom.jl")
  export Atom

  include("./src/Bond.jl")
  export Bond

  include("./src/Angle.jl")
  export Angle

  include("./src/Dihedral.jl")
  export Dihedral

  include("./src/Improper.jl")
  export Improper

  include("./src/readpsf.jl")
  export readpsf

  include("./src/readprm.jl")
  export readprm!

  include("./src/getpdbcoords.jl")
  export getpdbcoords!

  include("./src/coulomb.jl")
  export coulomb

  include("./src/vdw.jl")
  export vdw

  include("./src/nhbonds.jl")
  export nhbonds

end


text:
export ft_getmodemstatus
export CTS,DSR,RI,DCD
export OE,PE,FE,BI

#Modem Status (see ft_getmodemstatus)
const CTS = 0x10
const DSR = 0x20
const RI = 0x40
const DCD = 0x80

#Line Status (see ft_getmodemstatus)
const OE = 0x02
const PE = 0x04
const FE = 0x08
const BI = 0x10

function ft_getmodemstatus(ft_handle::Culong)
  modemstatus = Ref{Cuint}()
  ft_status = ccall((:FT_GetModemStatus, d2xx),
                     Cuint,
                     (Culong, Ref{Cuint}),
                     ft_handle, modemstatus)
  checkstatus(ft_status)
  return modemstatus[]
end


text:
@testset "dynamic_regularization" begin
  stats1 = ripqp(
    QuadraticModel(qps1),
    iconf = InputConfig(sp = K2LDLParams(regul = :dynamic), history = true),
    display = false,
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    iconf = InputConfig(sp = K2LDLParams(regul = :dynamic)),
    display = false,
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(
    QuadraticModel(qps3),
    iconf = InputConfig(sp = K2LDLParams(regul = :dynamic)),
    display = false,
  )
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "centrality_corrections" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(kc = -1), display = false) # automatic centrality corrections computation
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(kc = 2), display = false)
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(kc = 2), display = false)
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "refinement" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(refinement = :zoom), display = false) # automatic centrality corrections computation
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(refinement = :ref), display = false) # automatic centrality corrections computation
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    iconf = InputConfig(mode = :multi, scaling = false, refinement = :multizoom),
    display = false,
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(
    QuadraticModel(qps3),
    iconf = InputConfig(mode = :multi, refinement = :multiref),
    display = false,
  )
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "IPF" begin
  stats1 = ripqp(
    QuadraticModel(qps1),
    display = false,
    iconf = InputConfig(solve_method = :IPF, mode = :multi),
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2), display = false, iconf = InputConfig(solve_method = :IPF))
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(solve_method = :IPF, sp = K2_5LDLParams(), refinement = :zoom),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats1 = ripqp(
    QuadraticModel(qps1),
    display = false,
    iconf = InputConfig(solve_method = :IPF, mode = :multi, refinement = :multiref),
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(
      solve_method = :IPF,
      mode = :multi,
      refinement = :multizoom,
      sp = K2_5LDLParams(),
    ),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(
      solve_method = :IPF,
      refinement = :zoom,
      sp = K2LDLParams(regul = :dynamic),
    ),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable
end


text:
"""
    addColRainYday!(df::DataFrame, ColPct::String, NewCol::String, ColNum::Int64)

Adds a column of 0's and 1's indicating whether it has rained or not the previous day.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColPct::String`: The name of the column of `df` containing the total precipitation of a day.
- `NewCol::String`: The name of the new column of `df` that will contain the 0's and 1's.
- `ColNum::Int64`: The column number of `NewCol`.
- `reducedYear::Bool=true`: If true, the first day of the year is given a 0. Default = true.

# Implementation

This function uses the column indicating the total amount of rain for one day to determine if it has rained or not.
One must be aware that it uses a daily partitioning. If reducedYear = true, then the first day of the year is not treated
as the day following the last day of the previous year. For example, if the dataset contains data for May, June and July, the
1st of May will automatically be given a 0 in the new column.

"""

function addColRainYday!(df::DataFrame, ColPct::String, NewCol::String, ColNum::Int64, reducedYear::Bool=true)
    hier = []
    if reducedYear
        for i in 1:nrow(df)
            if i==1
                push!(hier, 0)
            else
                if df[i-1,Symbol(ColPct)] .> 0 && year(df[i,:].Date) .== year(df[i-1,:].Date)
                    push!(hier, 1)
                else
                    push!(hier, 0)
                end
            end
        end
    else
        for i in 1:nrow(df)
            if i==1
                push!(hier, 0)
            else
                if df[i-1,Symbol(ColPct)] .> 0
                    push!(hier, 1)
                else
                    push!(hier, 0)
                end
            end
        end
    end
    insertcols!(df, ColNum, Symbol(NewCol) => hier)
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
@testset "farmer" begin
    NS = 3  # number of scenarios
    probability = ones(3) / 3
    CROPS = 1:3 # set of crops (wheat, corn and sugar beets, resp.)
    PURCH = 1:2 # set of crops to purchase (wheat and corn, resp.)
    SELL = 1:4  # set of crops to sell (wheat, corn, sugar beets under 6K and those over 6K)
    Cost = [150 230 260]    # cost of planting crops
    Budget = 500            # budget capacity
    Purchase = [238 210]    # purchase price
    Sell = [170 150 36 10]  # selling price
    Yield = [3.0 3.6 24.0; 2.5 3.0 20.0; 2.0 2.4 16.0]
    Minreq = [200 240 0]    # minimum crop requirement

    @testset "MILP" begin
        # This creates a Lagrange dual problem for each scenario s.
        function create_scenario_model(s::Int64)
            m = Model(GLPK.Optimizer)
            @variable(m, 0 <= x[i=CROPS] <= 500, Int)
            @variable(m, y[j=PURCH] >= 0)
            @variable(m, w[k=SELL] >= 0)

            @objective(m, Min,
                probability[s] * sum(Cost[i] * x[i] for i=CROPS)
                + probability[s] * sum(Purchase[j] * y[j] for j=PURCH)
                - probability[s] * sum(Sell[k] * w[k] for k=SELL))

            @constraint(m, sum(x[i] for i=CROPS) <= Budget)
            @constraint(m, [j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
            @constraint(m, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
            @constraint(m, w[3] <= 6000)
            return m
        end

        models = Dict{Int,JuMP.Model}(s => create_scenario_model(s) for s in 1:NS)
        coupling_variables = Vector{DD.CouplingVariableRef}()
        for s in 1:NS
            model = models[s]
            xref = model[:x]
            for i in CROPS
                push!(coupling_variables, DD.CouplingVariableRef(s, i, xref[i]))
            end
        end

        @testset "Subgradient Method" begin
            # Create DualDecomposition instance.
            algo = DD.LagrangeDual()

            # Add Lagrange dual problem for each scenario s.
            for s in 1:NS
                DD.add_block_model!(algo, s, models[s])
            end

            # Set nonanticipativity variables as an array of symbols.
            DD.set_coupling_variables!(algo, coupling_variables)

            # Lagrange master method
            LM = DD.SubgradientMaster()
            LM.maxiter = 1000
            LM.step_size = (method) -> (abs(-108390 - method.f) / norm(method.∇f)^2) # Polyak's

            @testset "objective limit" begin
                DD.set_obj_limit!(LM, -110000.)
                DD.run!(algo, LM)
                @show DD.dual_objective_value(algo)
                @test DD.dual_objective_value(algo) > -110000.
            end
            DD.set_obj_limit!(LM, -100000.)

            # Solve the problem with the solver; this solver is for the underlying bundle method.
            DD.run!(algo, LM)

            @show DD.dual_objective_value(algo)
            @show DD.dual_solution(algo)
            @test isapprox(DD.dual_objective_value(algo), -108390, rtol=1e-3)
        end

        @testset "Proximal Method" begin
            # Create DualDecomposition instance.
            algo = DD.LagrangeDual()

            # Add Lagrange dual problem for each scenario s.
            for s in 1:NS
                DD.add_block_model!(algo, s, models[s])
            end

            # Set nonanticipativity variables as an array of symbols.
            DD.set_coupling_variables!(algo, coupling_variables)

            # Lagrange master method
            LM = DD.BundleMaster(BM.ProximalMethod, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))

            @testset "objective limit" begin
                DD.set_obj_limit!(LM, -111000.)
                DD.run!(algo, LM)
                @show DD.dual_objective_value(algo)
                @test DD.dual_objective_value(algo) > -111000.
            end
            DD.set_obj_limit!(LM, -100000.)

            # Solve the problem with the solver; this solver is for the underlying bundle method.
            DD.run!(algo, LM)

            @show DD.dual_objective_value(algo)
            @show DD.dual_solution(algo)
            @test isapprox(DD.dual_objective_value(algo), -108390, rtol=1e-3)
        end

        @testset "Trust Region Method" begin
            # Create DualDecomposition instance.
            algo = DD.LagrangeDual()

            # Add Lagrange dual problem for each scenario s.
            for s in 1:NS
                DD.add_block_model!(algo, s, models[s])
            end

            # Set nonanticipativity variables as an array of symbols.
            DD.set_coupling_variables!(algo, coupling_variables)

            # Lagrange master method
            LM = DD.BundleMaster(BM.TrustRegionMethod, GLPK.Optimizer)

            @testset "objective limit" begin
                DD.set_obj_limit!(LM, -111000.)
                DD.run!(algo, LM)
                @show DD.dual_objective_value(algo)
                @test DD.dual_objective_value(algo) > -111000.
            end
            DD.set_obj_limit!(LM, -100000.)

            # Solve the problem with the solver; this solver is for the underlying bundle method.
            DD.run!(algo, LM)

            @show DD.dual_objective_value(algo)
            @show DD.dual_solution(algo)
            @test isapprox(DD.dual_objective_value(algo), -108390, rtol=1e-3)
        end
    end

    @testset "QP" begin
        function deterministic_model_bound()
            m = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
            @variable(m, 0 <= x[i=CROPS] <= 500)
            @variable(m, y[s=1:NS,j=PURCH] >= 0)
            @variable(m, w[s=1:NS,k=SELL] >= 0)

            @objective(m, Min,
                sum(Cost[i] * x[i]^2 for i=CROPS)
                + sum(probability[s] * Purchase[j] * y[s,j] for s=1:NS, j=PURCH)
                - sum(probability[s] * Sell[k] * w[s,k] for s=1:NS, k=SELL))

            @constraint(m, sum(x[i] for i=CROPS) <= Budget)
            @constraint(m, [s=1:NS,j=PURCH], Yield[s,j] * x[j] + y[s,j] - w[s,j] >= Minreq[j])
            @constraint(m, [s=1:NS], Yield[s,3] * x[3] - w[s,3] - w[s,4] >= Minreq[3])
            @constraint(m, [s=1:NS], w[s,3] <= 6000)
            optimize!(m)
            objval = Inf
            if termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
                objval = objective_value(m)
            end
            return objval
        end

        # This creates a Lagrange dual problem for each scenario s.
        function create_scenario_model(s::Int64)
            m = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
            @variable(m, 0 <= x[i=CROPS] <= 500)
            @variable(m, y[j=PURCH] >= 0)
            @variable(m, w[k=SELL] >= 0)

            @objective(m, Min,
                probability[s] * sum(Cost[i] * x[i]^2 for i=CROPS)
                + probability[s] * sum(Purchase[j] * y[j] for j=PURCH)
                - probability[s] * sum(Sell[k] * w[k] for k=SELL))

            @constraint(m, sum(x[i] for i=CROPS) <= Budget)
            @constraint(m, [j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
            @constraint(m, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
            @constraint(m, w[3] <= 6000)
            return m
        end

        models = Dict{Int,JuMP.Model}(s => create_scenario_model(s) for s in 1:NS)
        coupling_variables = Vector{DD.CouplingVariableRef}()
        for s in 1:NS
            model = models[s]
            xref = model[:x]
            for i in CROPS
                push!(coupling_variables, DD.CouplingVariableRef(s, i, xref[i]))
            end
        end

        objval = deterministic_model_bound()
        @show objval
        @test objval < Inf

        @testset "Proximal Method" begin
            # Create DualDecomposition instance.
            algo = DD.LagrangeDual()

            # Add Lagrange dual problem for each scenario s.
            for s in 1:NS
                DD.add_block_model!(algo, s, models[s])
            end

            # Set nonanticipativity variables as an array of symbols.
            DD.set_coupling_variables!(algo, coupling_variables)

            # Change parameters
            params = BM.Parameters()
            BM.set_parameter(params, "ϵ_s", 1.e-4)
            @test BM.get_parameter(params, "ϵ_s") == 1.e-4

            # Lagrange master method
            LM = DD.BundleMaster(BM.ProximalMethod, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0), params)

            # Solve the problem with the solver; this solver is for the underlying bundle method.
            DD.run!(algo, LM)

            @show DD.dual_objective_value(algo)
            @show DD.dual_solution(algo)
            @test isapprox(DD.dual_objective_value(algo), objval, rtol=1e-3)
        end

        @testset "Trust Region Method" begin
            # Create DualDecomposition instance.
            algo = DD.LagrangeDual()

            # Add Lagrange dual problem for each scenario s.
            for s in 1:NS
                DD.add_block_model!(algo, s, models[s])
            end

            # Set nonanticipativity variables as an array of symbols.
            DD.set_coupling_variables!(algo, coupling_variables)

            # Lagrange master method
            LM = DD.BundleMaster(BM.TrustRegionMethod, GLPK.Optimizer)

            # Solve the problem with the solver; this solver is for the underlying bundle method.
            DD.run!(algo, LM)

            @show DD.dual_objective_value(algo)
            @show DD.dual_solution(algo)
            @test isapprox(DD.dual_objective_value(algo), objval, rtol=1e-3)
        end
    end
end


text:
using Indexing
if VERSION < v"0.7-"
    using Base.Test
else
    using Test
end

@testset "getindices" begin
    d = Dict(:a => "Alice", :b => "Bob", :c => "Charlie")
    @test getindices(d, [:a, :c]) == ["Alice", "Charlie"]
    @test getindices(d, (:a, :c)) == ("Alice", "Charlie")
    @test getindices(d, Dict(:aa => :a, :cc => :c)) == Dict(:aa => "Alice", :cc => "Charlie")
    @static if VERSION > v"0.7-"
        @test getindices(d, (aa = :a, cc = :c)) == (aa = "Alice", cc = "Charlie")
    end
    @test getindices(d, :) == d

    v = [11, 12, 13]
    @test (getindices(v, 2)::Array{Int, 0})[] == 12
    @test getindices(v, [1, 3]) == [11, 13]
    @test getindices(v, Dict(:a => 1, :c => 3)) == Dict(:a => 11, :c => 13)
    @test getindices(v, (1, 3)) === (11, 13)
    @static if VERSION > v"0.7-"
        @test getindices(v, (a = 1, c = 3)) === (a = 11, c = 13)
    end
    @test getindices(v, :) == v

    t = (11, 12, 13)
    @test getindices(t, [1, 3]) == [11, 13]
    @test getindices(t, Dict(:a => 1, :c => 3)) == Dict(:a => 11, :c => 13)
    @test getindices(t, (1, 3)) == (11, 13)
    @static if VERSION > v"0.7-"
        @test getindices(t, (a = 1, c = 3)) === (a = 11, c = 13)
    end
    @test getindices(t, :) == t

    @static if VERSION > v"0.7-"
        nt = (a = 1, b = 2.0, c = "three")
        @test getindices(nt, [:a, :c]) == [1, "three"]
        @test getindices(nt, Dict(:aa => :a, :cc => :c)) == Dict(:aa => 1, :cc => "three")
        @test getindices(nt, (:a, :c)) == (1, "three")
        @test getindices(nt, (aa = :a, cc = :c)) == (aa = 1, cc = "three")
        @test getindices(nt, :) == nt
    end
end

@testset "setindices!" begin
    d = Dict(:a => "Alice", :b => "Bob", :c => "Charlie")
    d2 = copy(d)
    setindices!(d2, "Someone", [:a, :c])
    @test d2 == Dict(:a => "Someone", :b => "Bob", :c => "Someone")

    d3 = copy(d)
    setindices!(d3, "Someone", (:a, :c))
    @test d3 == Dict(:a => "Someone", :b => "Bob", :c => "Someone")

    d4 = copy(d)
    setindices!(d4, "Someone", Dict(:aa => :a, :cc => :c))
    @test d4 == Dict(:a => "Someone", :b => "Bob", :c => "Someone")

    @static if VERSION > v"0.7-"
        d5 = copy(d)
        setindices!(d5, "Someone", (aa = :a, cc = :c))
        @test d5 == Dict(:a => "Someone", :b => "Bob", :c => "Someone")
    end

    d6 = copy(d)
    setindices!(d6, "Someone", :)
    @test d6 == Dict(:a => "Someone", :b => "Someone", :c => "Someone")

    v = [11, 12, 13]
    v2 = copy(v)
    setindices!(v2, 20, [1, 3])
    @test v2 == [20, 12, 20]

    v3 = copy(v)
    setindices!(v3, 20, (1, 3))
    @test v3 == [20, 12, 20]

    v4 = copy(v)
    setindices!(v4, 20, Dict(:a => 1, :c => 3))
    @test v4 == [20, 12, 20]

    @static if VERSION > v"0.7-"
        v5 = copy(v)
        setindices!(v5, 20, (a = 1, c = 3))
        @test v5 == [20, 12, 20]
    end

    v6 = copy(v)
    setindices!(v6, 20, :)
    @test v6 == [20, 20, 20]
end

@testset "view" begin
    d = Dict(:a => "Alice", :b => "Bob", :c => "Charlie")
    d2 = copy(d)

    @test view(d, [:a, :c])::ViewArray == ["Alice", "Charlie"]
    @test view(d, Dict(:aa => :a, :cc => :c))::ViewDict == Dict(:aa => "Alice", :cc => "Charlie")

    av = view(d, [:a, :c])
    @test parent(av) === d
    @test av[1] == "Alice"
    @test Indexing.axes(av) === (Base.OneTo(2),)
    av[1] = "Someone"
    @test d == Dict(:a => "Someone", :b => "Bob", :c => "Charlie")

    dv = view(d2, Dict(:aa => :a, :cc => :c))
    @test parent(dv) === d2
    @test dv[:aa] == "Alice"
    @test keys(dv) ⊆ [:aa, :cc] # Probably will want to change this
    @test length(keys(dv)) == 2
    @test haskey(dv, :aa)
    dv[:aa] = "No-one"
    @test d2 == Dict(:a => "No-one", :b => "Bob", :c => "Charlie")

    v = [11, 12, 13]

    @test view(v, Dict(:a =>1 , :c => 3))::ViewDict == Dict(:a => 11, :c => 13)

    dv2 = view(v, Dict(:a =>1 , :c => 3))
    @test parent(dv2) === v
    @test dv2[:a] == 11
    @test keys(dv2) ⊆ [:a, :c]
    @test length(keys(dv2)) == 2
    @test haskey(dv2, :a)
    @test first(dv2) === (:a => 11) || first(dv2) === (:c => 13)
    dv2[:a] = 21
    @test v == [21, 12, 13]

    @test ViewArray(d, [:b, :c])::ViewArray == ["Bob", "Charlie"]
    @test ViewVector(d, [:b, :c])::ViewVector == ["Bob", "Charlie"]
    @test ViewMatrix(d, [:b :c; :c :b])::ViewMatrix == ["Bob" "Charlie"; "Charlie" "Bob"]
end


text:
using quantumoptics

# System parameters for a pumped decaying cavity
η = 0.9
κ = 1.

# System operators
basis = FockBasis(20)
a = destroy(basis)
aᵀ = create(basis)
n = number(basis)
H = η*(a+aᵀ)
J = [sqrt(κ)*a]

# Initial state
Ψ₀ = fockstate(basis, 10)
ρ₀ = Ψ₀ ⊗ dagger(Ψ₀)

# Integration options
Ntrajectories = 50
T = [0:0.1:10;]

# Master time evolution
tout, ρt_master = timeevolution.master(T, ρ₀, H, J)

# Averaged MCWF time evolution
ρt_average = Operator[ρ₀*0. for i=1:length(T)]
for i=1:Ntrajectories
    tout, Ψt = timeevolution.mcwf(T, Ψ₀, H, J)
    for j=1:length(T)
        ρt_average[j] += (Ψt[j] ⊗ dagger(Ψt[j]))/Ntrajectories
    end
end

# Example trajectories
tout_example1, Ψt_example1 = timeevolution.mcwf(T, Ψ₀, H, J; seed=UInt(1),
                                        display_beforeevent=true,
                                        display_afterevent=true)

tout_example2, Ψt_example2 = timeevolution.mcwf(T, Ψ₀, H, J; seed=UInt(2),
                                        display_beforeevent=true,
                                        display_afterevent=true)

tout_example3, Ψt_example3 = timeevolution.mcwf(T, Ψ₀, H, J; seed=UInt(3),
                                        display_beforeevent=true,
                                        display_afterevent=true)

using PyCall
@pyimport matplotlib.pyplot as plt

plt.subplot(2, 1, 1)
plt.plot(T, real(expect(n, ρt_master)))
plt.plot(T, real(expect(n, ρt_average)))
plt.xlabel("\$\\mathrm\{Time\} \$")
plt.ylabel("\$\\mathrm{Photon number}\$")

plt.subplot(2, 1, 2)
plt.plot(tout_example1, real(expect(n, Ψt_example1)))
plt.plot(tout_example2, real(expect(n, Ψt_example2)))
plt.plot(tout_example3, real(expect(n, Ψt_example3)))
plt.xlabel("\$\\mathrm\{Time\} \$")
plt.ylabel("\$\\mathrm{Photon number}\$")

plt.show()


text:
using RealPolynomialRoots
using Polynomials
using Test

@testset "Real roots" begin

    # produces a very close pair of roots
    # like 2^(-n*τ/2)
    function mignotte(n,τ)
        x = variable(SparsePolynomial)
        p = x^n - ((2^(τ ÷ 2)-1)*x - 1)^2
        coeffs(convert(Polynomial, p))
    end

    p = mignotte(7, 14)
    st = ANewDsc(p)
    @test length(st) == 3
    @test length(st.Unresolved) == 0

    p = mignotte(15, 14)
    st = ANewDsc(p)
    @test length(st) == 3
    @test length(st.Unresolved) == 0

    p = mignotte(63, 14)    # gets slow ~ 2secs
    st = ANewDsc(p)
    @test length(st) == 3
    @test length(st.Unresolved) == 0

    𝐩 = Polynomial(mignotte(7, 14))
    𝐱 = variable(Polynomial{BigInt})

    p = coeffs(prod(𝐩(𝐱+i) for i in -3:3))
    st = ANewDsc(p)
    @test length(st) == 7*3
    @test length(st.Unresolved) == 0

    𝐪 = 𝐩 * prod(𝐱^2 + i for i in 1:10)
    p = coeffs(𝐪)
    st = ANewDsc(p)
    @test length(st) == 3
    @test length(st.Unresolved) == 0

    N = 50
    ps = [big(i)/N for i in 1:N]
    𝐩 = fromroots(Polynomial, ps);
    p = coeffs(𝐩)
    st = ANewDsc(p)
    @test length(st) == 50
    @test length(st.Unresolved) == 0

    # These two tests show that polynomials may not have the roots
    # that they theoretically might expect to have due to round
    # off.
    # This passes with ^30
    𝐩 = setprecision(256) do
        delta = 1/big(10)^30
        fromroots(Polynomial,[1.0, 1+delta, 1 + 2delta])
    end
    p = coeffs(𝐩)
    st = ANewDsc(p)
    @test length(st) == 3
    @test length(st.Unresolved) == 0

    # This shows mathematical expectations may fail due to precision issues
    delta = 1e-5
    p = coeffs(fromroots(Polynomial,[1.0, 1+delta, 1 + 2delta]))
    st = ANewDsc(p)
    @test length(st) == 3 # works

    delta = 1e-8
    p = coeffs(fromroots(Polynomial,[1.0, 1+delta, 1 + 2delta]))
    st = ANewDsc(p)
    @test_broken length(st) == 3 # doesn't to find 3

    delta = 1e-8
    p = coeffs(fromroots(Polynomial, BigFloat[1.0, 1+delta, 1 + 2delta]))
    st = ANewDsc(p)
    @test length(st) == 3 # works

    ## test real_roots
    rroots = (refine_roots ∘ ANewDsc ∘ coeffs)

    𝐱 = variable(Polynomial)

    # d = 0
    𝐩 = zero(𝐱)
    rts = rroots(𝐩)
    @test isempty(rts)

    𝐩 = one(𝐱)
    rts = rroots(𝐩)
    @test isempty(rts)

    # d = 1

    𝐩 = 𝐱
    rts = rroots(𝐩)
    @test length(rts) == 1 && abs(rts[1]) <= eps()

    𝐩 = 𝐱 + 1
    rts = rroots(𝐩)
    @test Float64.(rts) ≈ [-1]

    for i ∈ (1e2, 1e5, 1e8) # bigger values are just close
        𝐩 = 𝐱 - i
        rts = rroots(𝐩)
        @test Float64.(rts) ≈ [i]
    end

    # d > 1, nreal = 1
    for i ∈ -2000:117:2000
        𝐩 = (𝐱^2+1)*(𝐱-i)
        rts = rroots(𝐩)
        @test Float64.(rts) ≈ [i]
    end

    # d = 2
    𝐩 = 𝐱*(𝐱-1)
    rts = rroots(𝐩)
    @test sort(Float64.(rts)) ≈ [0,1]

    # larger d
    xs = -10:2:10
    𝐩 = fromroots(Polynomial, xs)
    𝐱 = variable(𝐩)
    𝐪 = 𝐩 * (𝐱^2 + 1)*(𝐱^4 + 1)
    rts = rroots(𝐪)
    @test length(rts) == length(xs)
    @test maximum(abs, sort(rts) - xs) <= eps()

end


text:
import NetCDF
import NetCDF: NcDim,NcVar
"""
This function save a DataCube object `d` generated by genDataCube to a NetCDF file `filename` in the directory `savedir`. The resulting NetCDF will have 3
variables, `datacube`, which is the actual multivariate dataset to be analysed, `components` will store the independent components used
to generate the cube and `events` will contain a map of events that should be detected by the algorithm. All settings will also be written
 as to the NetCDF file as attributes.
"""
function save_cube(d::DataCube, savedir = joinpath(Pkg.dir(),"SurrogateCube","data"),filename="")
	# Test if output directory exists
	isdir(savedir) || mkdir(savedir)
	# Generate new file name if none is given
	if filename==""
		i=1
		while isfile(joinpath(savedir,"file$(i).nc"))
			i=i+1
		end
		filename = "file$(i).nc"
	else
		isfile(joinpath(savedir,filename)) && rm(joinpath(savedir,filename))
	end
	fullname=joinpath(savedir,filename)
	# Now create essential variables
	Ntime,Nlat,Nlon,Nvar = size(d.variables)
	Ncomp = size(d.components,4)

	#Create variables
	#First define variable attributes
	varatts = Dict()
	compatts= Dict()
	objectToAtt(varatts,d.varnoise,"VarNoise")
	objectToAtt(compatts,d.compnoise,"CompNoise")
	objectToAtt(compatts,d.baselines,"Baseline")
	objectToAtt(compatts,d.events,"Event")

	dvar = NcDim("variable", Nvar)
	dcomp= NcDim("component", Ncomp)
	dlon = NcDim("longitude", Nlon)
	dlat = NcDim("latitude", Nlat)
	dtim = NcDim("time", Ntime)

	vcube = NcVar("datacube",[dtim,dlat,dlon,dvar],atts=varatts)
	vcomp = NcVar("components",[dtim,dlat,dlon,dcomp],atts=compatts)
	veven = NcVar("events",[dtim,dlat,dlon,dcomp])

	nc = NetCDF.create(fullname,NcVar[vcube,vcomp,veven])

	NetCDF.putvar(nc,"datacube",d.variables)
	NetCDF.putvar(nc,"components",d.components)
	NetCDF.putvar(nc,"events",d.eventmap)

	NetCDF.close(nc)
end
export save_cube

function objectToAtt(atts::Dict,t,prestring)

	T=string(typeof(t))
	n=fieldnames(t)
	atts[prestring]=T
	for i=1:length(n)
		atts[@sprintf("%s_%s",prestring,string(n[i]))]=t.(n[i])
	end
end

# Define a special function for Trend Events
function objectToAtt(atts::Dict,t::TrendEvent,prestring)

	T=string(typeof(t))
	n=fieldnames(t.ev)
	atts[prestring]=T
	for i=1:length(n)
		atts[@sprintf("%s_%s",prestring,string(n[i]))]=t.ev.(n[i])
	end
end

function objectToAtt(atts::Dict,t::Vector,prestring)

	for i=1:length(t)
		objectToAtt(atts,t[i],@sprintf("%s_%d",prestring,i))
	end
end

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
### Define Structs
abstract type AbstractAMR end
abstract type AMRGrid <: VisClaw.AbstractAMR end

###################################
"""
Struct:
 storm data
"""
mutable struct Storm <: VisClaw.AMRGrid
    gridnumber::Integer
    AMRlevel::Integer
    mx::Integer
    my::Integer
    xlow::Float64
    ylow::Float64
    dx::Float64
    dy::Float64
    u::AbstractArray{Float64,2}
    v::AbstractArray{Float64,2}
    slp::AbstractArray{Float64,2}
    # Constructor
    VisClaw.Storm(gridnumber, AMRlevel, mx, my, xlow, ylow, dx, dy, u, v, slp) =
    new(gridnumber, AMRlevel, mx, my, xlow, ylow, dx, dy, u, v, slp)
end
###################################

###################################
"""
Struct: Water veloccity
"""
mutable struct Velocity <: VisClaw.AMRGrid
    gridnumber::Integer
    AMRlevel::Integer
    mx::Integer
    my::Integer
    xlow::Float64
    ylow::Float64
    dx::Float64
    dy::Float64
    u :: AbstractArray{Float64,2}
    v :: AbstractArray{Float64,2}
    vel :: AbstractArray{Float64,2}
    # Constructor
    VisClaw.Velocity(gridnumber, AMRlevel, mx, my, xlow, ylow, dx, dy, u, v, vel) =
    new(gridnumber, AMRlevel, mx, my, xlow, ylow, dx, dy, u, v, vel)
end
###################################

###################################
"""
Struct: Sea Surface Height
"""
mutable struct SurfaceHeight <: VisClaw.AMRGrid
    gridnumber::Integer
    AMRlevel::Integer
    mx::Integer
    my::Integer
    xlow::Float64
    ylow::Float64
    dx::Float64
    dy::Float64
    eta::AbstractArray{Float64,2}
    # Constructor
    VisClaw.SurfaceHeight(gridnumber, AMRlevel, mx, my, xlow, ylow, dx, dy, eta) =
    new(gridnumber, AMRlevel, mx, my, xlow, ylow, dx, dy, eta)
end
###################################

###################################
"""
Struct:
 time-seies of AMR data
"""
mutable struct AMR <: VisClaw.AbstractAMR
    nstep::Integer
    timelap::AbstractVector
    amr :: AbstractVector{Vector{VisClaw.AMRGrid}}
    unittime :: Symbol
    # Constructor
    VisClaw.AMR(nstep, timelap, amr) = new(nstep, timelap, amr, :second)
end
###################################

abstract type AbstractTopo end
###################################
"""
Struct:
 topography and bathymetry
"""
struct Topo <: AbstractTopo
    ncols :: Integer
    nrows :: Integer
    x :: AbstractVector
    y :: AbstractVector
    dx :: Float64
    dy :: Float64
    elevation :: AbstractArray
    # Constructor
    VisClaw.Topo(ncols, nrows, x, y, dx, dy, elevation) =
             new(ncols, nrows, x, y, dx, dy, elevation)
end
###################################

##########################################################
"""
Struct:
 seafloor deformation for tsunami computation
"""
struct DTopo <: AbstractTopo
    mx :: Integer
    my :: Integer
    x :: AbstractVector
    y :: AbstractVector
    dx :: Float64
    dy :: Float64
    mt :: Integer
    t0 :: Float64
    dt :: Float64
    deform :: AbstractArray{Float64}
    # Constructor
    VisClaw.DTopo(mx,my,x,y,dx,dy,mt,t0,dt,deform) = new(mx,my,x,y,dx,dy,mt,t0,dt,deform)
end
##########################################################

##########################################################
"""
Struct:
 region
"""
struct Region
    level :: AbstractVector
    tlims :: Tuple
    xlims :: Tuple
    ylims :: Tuple
    # Constructor
    VisClaw.Region(level, tlims, xlims, ylims) = new(level, tlims, xlims, ylims)
end
##########################################################

########################################
"""
Struct: parameters in geoclaw.data
"""
struct GeoParam
    cs :: Integer # coordinate system
    p0:: Float64 # ambient pressure
    R :: Float64 # earth radious
    eta0 :: Float64 # sea level
    n ::Float64 # manning coafficient
    dmin :: Float64 # dry tolerance
    # Constructor
    VisClaw.GeoParam() = new(2,101300.0,6367500.0,0.0,0.025,0.001)
    VisClaw.GeoParam(cs,p0,R,eta0,n,dmin) = new(cs,p0,R,eta0,n,dmin)
end
########################################

########################################
"""
Struct: parameters in amr.data
"""
struct AMRParam
    maxlevel::Integer
    # Constructor
    VisClaw.AMRParam(maxlevel) = new(maxlevel)
end
########################################

########################################
"""
Struct: parameters in surge.data
"""
struct SurgeParam
    windindex::Integer
    slpindex::Integer
    stormtype::Integer
    # Constructor
    VisClaw.SurgeParam() = new(5,7,1)
    VisClaw.SurgeParam(windindex,slpindex,stormtype) = new(windindex,slpindex,stormtype)
end
########################################

########################################
"""
Struct: gauge data
"""
mutable struct Gauge
    label :: AbstractString # Name
    id :: Integer # gauge id
    nt :: Integer # number of time step
    loc :: AbstractVector{Float64} # gauge location
    AMRlevel :: AbstractVector{Integer}
    time :: AbstractVector # time
    eta :: AbstractVector{Float64} # surface
    u :: AbstractVector{Float64} # u
    v :: AbstractVector{Float64} # v
    unittime :: Symbol
    # Constructor
    VisClaw.Gauge(label,id,nt,loc,AMRlevel,time,eta) = new(label,id,nt,loc,AMRlevel,time,eta,[],[], :second)
    VisClaw.Gauge(label,id,nt,loc,AMRlevel,time,eta,u,v) = new(label,id,nt,loc,AMRlevel,time,eta,u,v, :second)
end
########################################

########################################
"""
Struct: max values at a gauge
"""
mutable struct Gaugemax
    label :: AbstractString # Name
    id :: Integer # gauge id
    loc :: AbstractVector{Float64} # gauge location
    AMRlevel :: Integer
    eta :: Float64
    vel :: Float64
    t_eta
    t_vel
    unittime :: Symbol
    # Constructor
    VisClaw.Gaugemax(label,id,loc,AMRlevel,eta,vel,t_eta,t_vel,unittime) = new(label,id,loc,AMRlevel,eta,vel,t_eta,t_vel,unittime)
end
########################################

########################################
"""
Struct: Fixed grid
"""
struct FixedGrid
    id :: Integer
    style :: Integer
    nval :: Integer
    ## point_style = 2, 4
    nx :: Integer
    ny :: Integer
    xlims :: Tuple
    ylims :: Tuple
    ## point_style = 0, 1, 4
    npts :: Integer
    x :: AbstractVector{Float64}
    y :: AbstractVector{Float64}
    ## point_style = 4
    flag :: AbstractArray

    # Constructor
    ## point_style = 0, 1
    VisClaw.FixedGrid(id,style,nval,npts,x,y) = new(id,style,nval,0,0,(NaN,NaN),(NaN,NaN),npts,x,y,[])
    ## point_style = 2
    VisClaw.FixedGrid(id,style,nval,nx,ny,xlims,ylims) = new(id,style,nval,nx,ny,xlims,ylims,0,[NaN],[NaN],[])
    ## point_style = 3
    VisClaw.FixedGrid(id,style,nval,nx,ny,xlims,ylims,npts,x,y) = new(id,style,nval,nx,ny,xlims,ylims,npts,x,y,[])
    ## point_style = 4
    VisClaw.FixedGrid(id,style,nval,nx,ny,xlims,ylims,npts,x,y,flag) = new(id,style,nval,nx,ny,xlims,ylims,npts,x,y,flag)

end
########################################

########################################
"""
Struct: fgmax values
"""
mutable struct FGmax
    topo :: AbstractArray{Float64}
    D :: AbstractArray{Float64}
    v :: AbstractArray{Float64}
    M :: AbstractArray{Float64}
    Mflux :: AbstractArray{Float64}
    Dmin :: AbstractArray{Float64}
    tD :: AbstractArray
    tv :: AbstractArray
    tM :: AbstractArray
    tMflux :: AbstractArray
    tDmin :: AbstractArray
    tarrival :: AbstractArray
    unittime :: Symbol

    # Constructor
    VisClaw.FGmax(topo,D,tD,tarrival) =
              new(topo,D,emptyF,emptyF,emptyF,emptyF,tD,emptyF,emptyF,emptyF,emptyF,tarrival,:second)
    VisClaw.FGmax(topo,D,v,tD,tv,tarrival) =
              new(topo,D,v,emptyF,emptyF,emptyF,tD,tv,emptyF,emptyF,emptyF,tarrival,:second)
    VisClaw.FGmax(topo,D,v,M,Mflux,Dmin,tD,tv,tM,tMflux,tDmin,tarrival) =
              new(topo,D,v,M,Mflux,Dmin,tD,tv,tM,tMflux,tDmin,tarrival,:second)
end
########################################

########################################
"""
Struct: track data container
"""
mutable struct Track
    timelap :: AbstractVector
    lon :: AbstractVector{Float64}
    lat :: AbstractVector{Float64}
    direction :: AbstractVector{Float64}
    unittime :: Symbol

    # Constructor
    VisClaw.Track(lon,lat) = new(empty([0.0]),lon,lat,empty([0.0]),:second)
    VisClaw.Track(timelap,lon,lat,direction) = new(timelap,lon,lat,direction,:second)
end
########################################


text:
using NonconvexMMA, LinearAlgebra, Test, Zygote, FiniteDifferences
const FDM = FiniteDifferences

f(x::AbstractVector) = sqrt(x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

_m = Model(f)
addvar!(_m, [0.0, 0.0], [10.0, 10.0])
add_ineq_constraint!(_m, x -> g(x, 2, 0))
add_ineq_constraint!(_m, x -> g(x, -1, 1))
m = NonconvexCore.tovecmodel(_m)[1]

x = fill(0.5, 2)

@testset "Gradient and hessian of MMA approximation" begin
    exactf = NonconvexMMA.getobjective(m)
    approxf = NonconvexMMA.MMAApprox(exactf, x)
    val1, grad1 = NonconvexCore.value_gradient(exactf, x)
    val2, grad2 = exactf(x), Zygote.gradient(exactf, x)[1]
    val3, grad3 = NonconvexCore.value_gradient(approxf, x)
    val4, grad4 = approxf(x), Zygote.gradient(approxf, x)[1]
    grad5 = FDM.grad(central_fdm(5, 1), approxf, x)[1]
    @test val1 ≈ val2 ≈ val3 ≈ val4
    @test grad1 ≈ grad2 ≈ grad3 ≈ grad4 ≈ grad5
    for i in 1:100
        y = rand(2)*10
        H = Zygote.hessian(approxf, y)
        for i in 1:2, j in 1:2
            if i == j
                @test H[i, j] >= 0
            else
                @test H[i, j] == 0
            end
        end
    end

    exactf = NonconvexMMA.getobjective(m)
    approxf = NonconvexMMA.MMAApprox(exactf, x)
    val1, grad1 = NonconvexCore.value_gradient(exactf, x)
    val2, grad2 = exactf(x), Zygote.gradient(exactf, x)[1]
    val3, grad3 = NonconvexCore.value_gradient(approxf, x)
    val4, grad4 = approxf(x), Zygote.gradient(approxf, x)[1]
    grad5 = FDM.grad(central_fdm(5, 1), approxf, x)[1]
    @test val1 ≈ val2 ≈ val3 ≈ val4
    @test grad1 ≈ grad2 ≈ grad3 ≈ grad4 ≈ grad5
    for i in 1:100
        y = rand(2)*10
        H = Zygote.hessian(approxf, y)
        for i in 1:2, j in 1:2
            if i == j
                @test H[i, j] >= 0
            else
                @test H[i, j] == 0
            end
        end
    end

    exactf = NonconvexMMA.getineqconstraints(m)
    approxf = NonconvexMMA.MMAApprox(exactf, x)
    val1, jac1 = NonconvexCore.value_jacobian(exactf, x)
    val2 = exactf(x)
    jac2 = vcat([Zygote.gradient(x -> exactf(x)[i], x)[1]' for i in 1:length(val2)]...)
    val3, jac3 = NonconvexCore.value_jacobian(approxf, x)
    val4 = approxf(x)
    jac4 = vcat([Zygote.gradient(x -> approxf(x)[i], x)[1]' for i in 1:length(val2)]...)
    jac5 = vcat([FDM.grad(central_fdm(5, 1), x -> approxf(x)[i], x)[1]' for i in 1:length(val2)]...)

    @test val1 ≈ val2 ≈ val3 ≈ val4
    @test jac1 ≈ jac2 ≈ jac3 ≈ jac4 ≈ jac5
end
@testset "Jacobian of MMA approximation" begin
    approx = NonconvexMMA.MMAApproxModel(m, x)
    exactfs = approx.objective_ineq_constraints
    approxfs = approx.approx_objective_ineq_constraints
    @test exactfs(x) ≈ approxfs(x)

    val1, jac1 = NonconvexCore.value_jacobian(exactfs, x)
    val2, jac2 = NonconvexCore.value_jacobian(approxfs, x)
    @test val1 ≈ val2
    @test jac1 ≈ jac2
end

@testset "Dual objective" begin
    approx = NonconvexMMA.MMAApproxModel(m, x)
    dualmodel = NonconvexMMA.MMADualModel(approx)
    dualobj = dualmodel.obj
    λ = [2.0, 3.0]

    @testset "Value and gradient wrt λ" begin
        val1, grad1 = NonconvexCore.value_gradient(dualobj, λ)
        optimalx = copy(NonconvexMMA.optimizeprimal!(dualobj, λ))
        fg = approx.approx_objective_ineq_constraints
        val2 = fg(optimalx)' * [1.0; λ]
        grad2 = FDM.grad(central_fdm(5, 1), dualobj, λ)[1]
        @test val1 ≈ val2
        @test grad1 ≈ grad2

        grad1 = Zygote.gradient(λ -> dualobj(optimalx, λ), λ)[1]
        grad2 = FDM.grad(central_fdm(5, 1), λ -> dualobj(optimalx, λ), λ)[1]
        @test grad1 ≈ grad2 ≈ fg(optimalx)[2:end]
    end

    @testset "Value and gradient wrt x" begin
        optimalx = copy(NonconvexMMA.optimizeprimal!(dualobj, λ))
        grad1 = Zygote.gradient(x -> dualobj(x, λ), optimalx)[1]
        grad2 = FDM.grad(central_fdm(5, 1), x -> dualobj(x, λ), optimalx)[1]
        @test grad1 ≈ grad2
        @test all(1:length(x)) do j
            abs(grad1[j]) <= 1e-6 ||
                optimalx[j] == 0 && grad1[j] > 0 ||
                optimalx[j] == 10 && grad1[j] < 0
        end

        approxfs = approx.approx_objective_ineq_constraints
        val, jac = NonconvexCore.value_jacobian(approxfs, optimalx)
        @test grad1 ≈ jac' * [1; λ]
    end
end


text:
using Test, GraphGrinder, Flux

@testset matchvertex2edges begin
    x = rand(3,3)
    Δ = Flux.data(Flux.Tracker.gradient(x -> sum(sin.(matchvertex2edges(x, edges(g)))), x)[1])
    nΔ = Flux.Tracker.ngradient(x -> sum(sin.(matchvertex2edges(x, edges(g)))), x)[1]
    @test isapprox(Δ, nΔ, atol = 1e-4)
end


text:
using Metatheory
using Metatheory.EGraphs
using Metatheory.Library

Context = commutative_monoid(:∧, :(:I)) ∪ @theory begin
    (A.Ob ⊢ id(A)) == (A.Ob ∧ id(A).(A→A) ⊢ id(A))
end

println(Context)

using Test
@testset "Context" begin
    @testset "Commutative Monoid in the Context" begin
        @test @areequal Context (A ∧ B⊢ f) (B ∧ A ⊢ f)
        @test @areequal Context (A ∧ (B ∧ C) ⊢ f) ((A∧B) ∧ C ⊢ f)
        @test @areequal Context (A ∧ (B ∧ C) ⊢ f) ((B∧A) ∧ C ⊢ f)
        @test @areequal Context (A ∧ (B ∧ C) ⊢ f) ((C∧B) ∧ A ⊢ f)
        @test @areequal Context ((C ∧ A) ∧ B ⊢ f) ((C ∧ B) ∧ A ⊢ f)
        @test @areequal Context ((A ∧ A) ∧ B ⊢ f) (B ∧ A ⊢ f)
        @test @areequal Context (B∧(A ∧ A) ∧ B ⊢ f) (B ∧ A ⊢ f)
    end
    @testset "Not equal things" begin
        @test ! @areequal Context (A∧B ⊢ A)  (A ⊢ A)
        @test ! @areequal Context (A ⊢ A)  (A∧B ⊢ A)
        @test ! @areequal Context (A ⊢ A)  (A ⊢ A∧B)
        @test ! @areequal Context (A ⊢ A)  (A ⊢ A ∨ B)
    end
    @testset "Unitality" begin
        @test @areequal Context A.Ob ⊢ id(A)        A.Ob∧ id(A).(A→A) ⊢ id(A)
        @test @areequal Context ((A.Ob ∧ A.Ob) ⊢ id(A))        ((A.Ob) ⊢ id(A))
        @test @areequal Context (A.Ob ∧ A.Ob) ⊢ id(A)        A.Ob ∧id(A).(A→A)⊢ id(A)
        @test @areequal Context A.Ob ∧ B.Ob ⊢ id(A)        B.Ob ∧ A.Ob ∧ id(A).(A→A) ⊢ id(A)
        @test @areequal Context (A.Ob ∧ B.Ob ⊢ id(A))        A.Ob ∧ id(A).(A→A) ⊢ id(A)
        @test ! @areequal Context A.Ob ∧ B.Ob ⊢ id(A)  B.Ob ∧ id(B).(A→A)⊢ id(B)
    end
end


text:
#!/usr/bin/env julia

using Cooler

coolName = "MB286.mcool"
startCoord=1000
endCoord = 10000000
resolution = "resolutions/1000"

matrix = Cooler.matrix(coolName,startCoord,endCoord,resolution)
println(length(matrix))

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(sin(e+f*x)^m*(1+m+-1*sin(e+f*x)^2*(2+m)), x) == :(f^-1*sin(e+f*x)^(1+m)*cos(e+f*x))
@test integrate(sin(e+f*x)^5*(6+-7*sin(e+f*x)^2), x) == :(f^-1*sin(e+f*x)^6*cos(e+f*x))
@test integrate(sin(e+f*x)^4*(5+-6*sin(e+f*x)^2), x) == :(f^-1*sin(e+f*x)^5*cos(e+f*x))
@test integrate(sin(e+f*x)^3*(4+-5*sin(e+f*x)^2), x) == :(f^-1*sin(e+f*x)^4*cos(e+f*x))
@test integrate(sin(e+f*x)^2*(3+-4*sin(e+f*x)^2), x) == :(f^-1*sin(e+f*x)^3*cos(e+f*x))
@test integrate((2+-3*sin(e+f*x)^2)*sin(e+f*x), x) == :(f^-1*sin(e+f*x)^2*cos(e+f*x))
@test integrate(1+-2*sin(e+f*x)^2, x) == :(f^-1*cos(e+f*x)*sin(e+f*x))
@test integrate(-1*sin(e+f*x)^2*csc(e+f*x), x) == :(f^-1*cos(e+f*x))
@test integrate(-1*csc(e+f*x)^2, x) == :(f^-1*cot(e+f*x))
@test integrate(csc(e+f*x)^3*(-2+sin(e+f*x)^2), x) == :(f^-1*cot(e+f*x)*csc(e+f*x))
@test integrate(csc(e+f*x)^4*(-3+2*sin(e+f*x)^2), x) == :(f^-1*csc(e+f*x)^2*cot(e+f*x))
@test integrate(csc(e+f*x)^5*(-4+3*sin(e+f*x)^2), x) == :(f^-1*csc(e+f*x)^3*cot(e+f*x))
@test integrate(sin(e+f*x)^5*(A+C*sin(e+f*x)^2), x) == :(-1*f^-1*(A+C)*cos(e+f*x)+-1//5*f^-1*cos(e+f*x)^5*(A+3C)+1//3*f^-1*cos(e+f*x)^3*(2A+3C)+1//7*C*f^-1*cos(e+f*x)^7)

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.4.1 (a+b sin)^m (A+B sin+C sin^2).jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractDensityTransformTarget

Abstract type for probability density transformation targets.
"""
abstract type AbstractDensityTransformTarget end
export AbstractDensityTransformTarget

"""
    abstract type AbstractTransformToUnitspace <: AbstractDensityTransformTarget

Abstract type for density transformation targets that specify a
transformation into the unit hypercube.
"""
abstract type AbstractTransformToUnitspace <: AbstractDensityTransformTarget end
export AbstractTransformToUnitspace

"""
    abstract type AbstractTransformToInfinite <: AbstractDensityTransformTarget

Abstract type for density transformation targets that specify are
transformation into unbounded space.

The density that results from such a transformation must be unbounded in all
dimensions and that should taper off to zero somewhat smoothly in any
direction.
"""
abstract type AbstractTransformToInfinite <: AbstractDensityTransformTarget end
export AbstractTransformToInfinite

"""
    abstract type TransformAlgorithm

Abstract type for density transformation algorithms.
"""
abstract type TransformAlgorithm end
export TransformAlgorithm

"""
    bat_transform(
        target::AbstractDensityTransformTarget,
        density::AnyDensityLike,
        [algorithm::TransformAlgorithm]
    )::AbstractDensity

Transform `density` to another variate space defined/implied by `target`.

Returns a NamedTuple of the shape

```julia
(result = newdensity::AbstractDensity, trafo = vartrafo::Function, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_transform`, add methods to
    `bat_transform_impl` instead.
"""
function bat_transform end
export bat_transform

function bat_transform_impl end

function bat_transform(
    target::AbstractDensityTransformTarget,
    density::AnyDensityLike,
    algorithm::TransformAlgorithm = bat_default_withinfo(bat_transform, Val(:algorithm), target, density)
)
    r = bat_transform_impl(target, density, algorithm)
    result_with_args(r, (algorithm = algorithm,))
end

function argchoice_msg(::typeof(bat_transform), ::Val{:algorithm}, x::TransformAlgorithm)
    "Using transform algorithm $x"
end

"""
    struct NoDensityTransform <: AbstractDensityTransformTarget

The identity density transformation target, specifies that densities
should not be transformed.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct NoDensityTransform <: AbstractDensityTransformTarget end
export NoDensityTransform

"""
    struct DensityIdentityTransform <: TransformAlgorithm

A no-op density transform algorithm that leaves any density unchanged.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct DensityIdentityTransform <: TransformAlgorithm end
export DensityIdentityTransform

function bat_transform_impl(target::NoDensityTransform, density::AnyDensityLike, algorithm::DensityIdentityTransform)
    (result = convert(AbstractDensity, density), trafo = identity)
end

# ToDo: Merge PriorToUniform and PriorToGaussian into PriorTo{Uniform|Normal}.

"""
    struct PriorToUniform <: AbstractTransformToUnitspace

Specifies that posterior densities should be transformed in a way that makes
their pior equivalent to a uniform distribution over the unit hypercube.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct PriorToUniform <: AbstractTransformToUnitspace end
export PriorToUniform

_distribution_density_trafo(target::PriorToUniform, density::DistributionDensity) = DistributionTransform(Uniform, parent(density))

function bat_transform_impl(target::PriorToUniform, density::DistributionDensity{<:StandardUniformDist}, algorithm::DensityIdentityTransform)
    (result = density, trafo = identity)
end

"""
    struct PriorToGaussian <: AbstractTransformToInfinite

Specifies that posterior densities should be transformed in a way that makes
their pior equivalent to a standard multivariate normal distribution with an
identity covariance matrix.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct PriorToGaussian <: AbstractTransformToInfinite end
export PriorToGaussian

_distribution_density_trafo(target::PriorToGaussian, density::DistributionDensity) = DistributionTransform(Normal, parent(density))

function bat_transform_impl(target::PriorToGaussian, density::DistributionDensity{<:StandardNormalDist}, algorithm::DensityIdentityTransform)
    (result = density, trafo = identity)
end

"""
    struct FullDensityTransform <: TransformAlgorithm

Transform the density as a whole a given specified target space. Operations
that use the gradient of the density will require to the `log(abs(jacobian))`
of the transformation to be auto-differentiable.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct FullDensityTransform <: TransformAlgorithm end
export FullDensityTransform

_get_deep_prior_for_trafo(density::DistributionDensity) = density
_get_deep_prior_for_trafo(density::AbstractPosteriorDensity) = _get_deep_prior_for_trafo(getprior(density))
_get_deep_prior_for_trafo(density::RenormalizedDensity) = _get_deep_prior_for_trafo(parent(density))

function bat_transform_impl(target::Union{PriorToUniform,PriorToGaussian}, density::AbstractPosteriorDensity, algorithm::FullDensityTransform)
    orig_prior = _get_deep_prior_for_trafo(density)
    trafo = _distribution_density_trafo(target, orig_prior)
    (result = TransformedDensity(density, trafo, TDLADJCorr()), trafo = trafo)
end

function bat_transform_impl(target::Union{PriorToUniform,PriorToGaussian}, density::DistributionDensity, algorithm::FullDensityTransform)
    trafo = _distribution_density_trafo(target, density)
    (result = TransformedDensity(density, trafo, TDLADJCorr()), trafo = trafo)
end

"""
    struct PriorSubstitution <: TransformAlgorithm

Substitute the prior by a given distribution and transform the
likelihood accordingly. The `log(abs(jacobian))` of the transformation does
not need to be auto-differentiable even for operations that use the
gradient of the posterior.

Constructors:

* ```$(FUNCTIONNAME)()```
"""
struct PriorSubstitution <: TransformAlgorithm end
export PriorSubstitution

function bat_transform_impl(target::Union{PriorToUniform,PriorToGaussian}, density::DistributionDensity, algorithm::PriorSubstitution)
    trafo = _distribution_density_trafo(target, density)
    transformed_density = DistributionDensity(trafo.target_dist)
    (result = transformed_density, trafo = trafo)
end

function bat_transform_impl(target::Union{PriorToUniform,PriorToGaussian}, density::AbstractPosteriorDensity, algorithm::PriorSubstitution)
    orig_prior = getprior(density)
    orig_likelihood = getlikelihood(density)
    new_prior, trafo = bat_transform_impl(target, orig_prior, algorithm)
    new_likelihood = TransformedDensity(orig_likelihood, trafo, TDNoCorr())
    (result = PosteriorDensity(new_likelihood, new_prior), trafo = trafo)
end

function bat_transform_impl(target::Union{PriorToUniform,PriorToGaussian}, density::RenormalizedDensity, algorithm::PriorSubstitution)
    new_parent_density, trafo = bat_transform_impl(target, parent(density), algorithm)
    (result = RenormalizedDensity(new_parent_density, density.logrenormf), trafo = trafo)
end

# ToDo: Support bat_transform for vectors of variates and DensitySampleVector?

unshaping_trafo(::ArrayShape{Real, 1}) = identity
unshaping_trafo(vs::AbstractValueShape) = inverse(vs)

function transform_and_unshape(trafotarget::AbstractDensityTransformTarget, object::Any)
    transform_and_unshape(bat_transform, trafotarget, object)
end

function transform_and_unshape(bat_trafofunc::Function, trafotarget::AbstractDensityTransformTarget, object::Any)
    orig_denstiy = convert(AbstractDensity, object)
    trafoalg = bat_default(bat_trafofunc, Val(:algorithm), trafotarget, orig_denstiy)
    transformed_density, initial_trafo = bat_trafofunc(trafotarget, orig_denstiy, trafoalg)
    us_trafo = unshaping_trafo(varshape(transformed_density))
    result_density = us_trafo(transformed_density)
    result_trafo = us_trafo ∘ initial_trafo
    return result_density, result_trafo
end


text:
module SDFResults

export read_simulation, EPOCHSimulation, SDFFile,
    get_time, get_parameter, get_npart,
    domain_length, domain_volume, cell_volume,
    compute_L

using Unitful
using PhysicalConstants.CODATA2018: c_0, ε_0, μ_0, m_e, e
using Statistics, StatsBase
using BangBang
using Transducers, ThreadsX
using SDFReader
using SDFReader: get_units
using PICDataStructures
using EPOCHInput
using StaticArrays: MVector

include("files/sdffile.jl")
include("simulation.jl")
include("traits.jl")
include("read.jl")
include("utils.jl")

end


text:
#TODO: more tests
@testset "Model" begin
  adata, mdata, models = runmodel(param_file)
  @test size(mdata, 1) == 15
  @test maximum(mdata.step) == 14
end

@testset "replicates" begin
  replicates = 3
  adata, mdata, models = runmodel(param_file, replicates = replicates)
  size(mdata)
  @test size(mdata, 1) == 15*replicates
  @test maximum(mdata.step) == 14
  @test maximum(mdata.ensemble) == 3
end


text:
function a_ideal(model::SAFTgammaMie, conditions)
    0
end
function a_mono(model::SAFTgammaMie, conditions)
    0
end
function a_chain(model::SAFTgammaMie, conditions)
    0
end
function a_assoc(model::SAFTgammaMie, conditions)
    0
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function prodfactor(pv::PV, latitude::Real, longitude::Real; timeframe="hourly")

    if !(ismissing(pv.prod_factor_series_kw))
        return pv.prod_factor_series_kw
    end

    url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
        "&lat=", latitude , "&lon=", longitude, "&tilt=", pv.tilt,
        "&system_capacity=1", "&azimuth=", pv.azimuth, "&module_type=", pv.module_type,
        "&array_type=", pv.array_type, "&losses=", round(pv.losses*100, digits=3), "&dc_ac_ratio=", pv.dc_ac_ratio,
        "&gcr=", pv.gcr, "&inv_eff=", pv.inv_eff*100, "&timeframe=", timeframe, "&dataset=nsrdb",
        "&radius=", pv.radius
    )

    try
        @info "Querying PVWatts for prodfactor with " pv.name
        r = HTTP.get(url)
        response = JSON.parse(String(r.body))
        if r.status != 200
            error("Bad response from PVWatts: $(response["errors"])")
        end
        @info "PVWatts success."
        watts = collect(get(response["outputs"], "ac", []) / 1000)  # scale to 1 kW system (* 1 kW / 1000 W)
        if length(watts) != 8760
            @error "PVWatts did not return a valid production factor. Got $watts"
        end
        return watts
    catch e
        @error "Error occurred when calling PVWatts: $e"
    end
end

function prodfactor(g::AbstractGenerator; ts_per_hour::Int=1)
    return ones(8760 * ts_per_hour)
end

"""
    prodfactor(wind::Wind, latitude::Real, longitude::Real)

If the user does not provide their own prod_factor_series_kw for the Wind turbine, then this method creates
a production factor time-series using resource data and the System Advisor Model Wind module.
If the user does not provide the resource data, the latitude and longitude are used to get the resource data from the
Wind Toolkit.
"""
function prodfactor(wind::Wind, latitude::Real, longitude::Real, time_steps_per_hour::Int)

    if !(ismissing(wind.prod_factor_series_kw))
        return wind.prod_factor_series_kw
    end

    resources = []
    heights_for_sam = [wind.hub_height]

    if all(length(a) > 0 for a in [wind.temperature_celsius, wind.pressure_atmospheres, wind.wind_direction_degrees,
                                   wind.wind_meters_per_sec])
        push!(resources, [wind.temperature_celsius, wind.pressure_atmospheres, wind.wind_direction_degrees,
                          wind.wind_meters_per_sec])
    else  # download resource data from WindToolKit

        # Allowed hub heights in meters for the Wind Toolkit
        windtoolkit_hub_heights = [10, 40, 60, 80, 100, 120, 140, 160, 200]
        """
        SAM can interpolate the wind power if the wind.hub_height is not one of the windtoolkit_hub_heights.
            If we do need to interpolate then we need to provide SAM with the resources for two hub heights.
        """
        if !(wind.hub_height in windtoolkit_hub_heights)
            if wind.hub_height < minimum(windtoolkit_hub_heights)
                heights_for_sam = [windtoolkit_hub_heights[1]]
            elseif wind.hub_height > maximum(windtoolkit_hub_heights)
                heights_for_sam = [windtoolkit_hub_heights[end]]
            else
                upper_index = findfirst(x -> x > wind.hub_height, windtoolkit_hub_heights)
                heights_for_sam = [windtoolkit_hub_heights[upper_index-1], windtoolkit_hub_heights[upper_index]]
            end
        end
        # TODO validate against API with different hub heights (not in windtoolkit_hub_heights)

        for height in heights_for_sam
            url = string("https://developer.nrel.gov/api/wind-toolkit/v2/wind/wtk-srw-download",
                "?api_key=", nrel_developer_key,
                "&lat=", latitude , "&lon=", longitude,
                "&hubheight=", Int(height), "&year=", 2012
            )
            resource = []
            try
                @info "Querying Wind Toolkit for resource data ..."
                r = HTTP.get(url; retries=5)
                if r.status != 200
                    error("Bad response from Wind Toolkit: $(response["errors"])")
                end
                @info "Wind Toolkit success."

                resource = readdlm(IOBuffer(String(r.body)), ',', Float64, '\n'; skipstart=5);
                # columns: Temperature, Pressure, Speed, Direction (C, atm, m/s, Degrees)
                if size(resource) != (8760, 4)
                    @error "Wind Toolkit did not return valid resource data. Got an array with size $(size(resource))"
                end
            catch e
                @error "Error occurred when calling Wind Toolkit: $e"
            end
            push!(resources, resource)
        end
        resources = hcat(resources...)
    end  # done filling in resources (can contain one Vector of Vectors or multiple Vectors of Vectors)

    # Initialize SAM inputs
    global hdl = nothing
    sam_prodfactor = []

    # Corresponding size in kW for generic reference turbines sizes
    system_capacity_lookup = Dict(
        "large"=> 2000,
        "medium" => 250,
        "commercial"=> 100,
        "residential"=> 2.5
    )
    system_capacity = system_capacity_lookup[wind.size_class]

    # Corresponding rotor diameter in meters for generic reference turbines sizes
    rotor_diameter_lookup = Dict(
        "large" => 55*2,
        "medium" => 21.9*2,
        "commercial" => 13.8*2,
        "residential" => 1.85*2
    )
    wind_turbine_powercurve_lookup = Dict(
        "large" => [0, 0, 0, 70.119, 166.208, 324.625, 560.952, 890.771, 1329.664,
                    1893.213, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000,
                    2000, 2000, 2000, 2000, 2000, 2000],
        "medium"=> [0, 0, 0, 8.764875, 20.776, 40.578125, 70.119, 111.346375, 166.208,
                    236.651625, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250,
                    250, 250, 250, 250, 250],
        "commercial"=> [0, 0, 0, 3.50595, 8.3104, 16.23125, 28.0476, 44.53855, 66.4832,
                        94.66065, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100,
                        100, 100, 100, 100, 100],
        "residential"=> [0, 0, 0, 0.070542773, 0.1672125, 0.326586914, 0.564342188,
                        0.896154492, 1.3377, 1.904654883, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5,
                        2.5, 2.5, 2.5, 0, 0, 0, 0, 0, 0, 0]
    )

    try
        if Sys.isapple()
            libfile = "libssc.dylib"
        elseif Sys.islinux()
            libfile = "libssc.so"
        elseif Sys.iswindows()
            libfile = "ssc.dll"
        else
            @error """Unsupported platform for using the SAM Wind module.
                      You can alternatively provide the Wind.prod_factor_series_kw"""
        end

        global hdl = joinpath(dirname(@__FILE__), "..", "sam", libfile)
        wind_module = @ccall hdl.ssc_module_create("windpower"::Cstring)::Ptr{Cvoid}
        wind_resource = @ccall hdl.ssc_data_create()::Ptr{Cvoid}  # data pointer
        @ccall hdl.ssc_module_exec_set_print(0::Cint)::Cvoid

        @ccall hdl.ssc_data_set_number(wind_resource::Ptr{Cvoid}, "latitude"::Cstring, latitude::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(wind_resource::Ptr{Cvoid}, "longitude"::Cstring, longitude::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(wind_resource::Ptr{Cvoid}, "elevation"::Cstring, 0::Cdouble)::Cvoid  # not used in SAM
        @ccall hdl.ssc_data_set_number(wind_resource::Ptr{Cvoid}, "year"::Cstring, 2012::Cdouble)::Cvoid

        heights_array = []  # have to repeat heights for each resource column
        for h in heights_for_sam
            append!(heights_array, repeat([h], 4))
        end
        heights_array = convert(Array{Float64}, heights_array)
        @ccall hdl.ssc_data_set_array(wind_resource::Ptr{Cvoid}, "heights"::Cstring,
           heights_array::Ptr{Cdouble}, length(heights_array)::Cint)::Cvoid

        # setup column data types: temperature=1, pressure=2, degree=3, speed=4
        fields = collect(repeat(range(1, stop=4), length(heights_for_sam)))
        fields = convert(Array{Float64}, fields)
        @ccall hdl.ssc_data_set_array(wind_resource::Ptr{Cvoid}, "fields"::Cstring,
            fields::Ptr{Cdouble}, length(fields)::Cint)::Cvoid

        nrows, ncols = size(resources)
        t = [row for row in eachrow(resources)];
        t2 = reduce(vcat, t);
        # the values in python api are sent to SAM as vector (35040) with rows concatenated
        c_resources = [convert(Float64, t2[i]) for i in eachindex(t2)]
        @ccall hdl.ssc_data_set_matrix(wind_resource::Ptr{Cvoid}, "data"::Cstring, c_resources::Ptr{Cdouble},
            Cint(nrows)::Cint, Cint(ncols)::Cint)::Cvoid

        data = @ccall hdl.ssc_data_create()::Ptr{Cvoid}  # data pointer
        @ccall hdl.ssc_data_set_table(data::Ptr{Cvoid}, "wind_resource_data"::Cstring, wind_resource::Ptr{Cvoid})::Cvoid
        @ccall hdl.ssc_data_free(wind_resource::Ptr{Cvoid})::Cvoid

        # # can get the same values back with:
        # a = @ccall hdl.ssc_data_get_matrix(wind_resource::Ptr{Cvoid}, "data"::Cstring, Ref(Cint(4))::Ptr{Cint}, Ref(Cint(4))::Ptr{Cint})::Ptr{Cdouble}
        # unsafe_load(a, 1)
        # unsafe_load(a, 35040)

        # Scaler inputs
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_resource_shear"::Cstring, 0.14000000059604645::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_resource_turbulence_coeff"::Cstring,
            0.10000000149011612::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "system_capacity"::Cstring,
            system_capacity::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_resource_model_choice"::Cstring, 0::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "weibull_reference_height"::Cstring, 50::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "weibull_k_factor"::Cstring, 2::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "weibull_wind_speed"::Cstring, 7.25::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_turbine_rotor_diameter"::Cstring,
            rotor_diameter_lookup[wind.size_class]::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_turbine_hub_ht"::Cstring, wind.hub_height::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_turbine_max_cp"::Cstring, 0.44999998807907104::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_farm_losses_percent"::Cstring, 0::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "wind_farm_wake_model"::Cstring, 0::Cdouble)::Cvoid
        @ccall hdl.ssc_data_set_number(data::Ptr{Cvoid}, "adjust:constant"::Cstring, 0::Cdouble)::Cvoid

        # Array inputs
        speeds = convert(Array{Float64},
            [0., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25])
        @ccall hdl.ssc_data_set_array(data::Ptr{Cvoid}, "wind_turbine_powercurve_windspeeds"::Cstring,
            speeds::Ptr{Cdouble}, length(speeds)::Cint)::Cvoid

        powercurve = convert(Array{Float64}, wind_turbine_powercurve_lookup[wind.size_class])
        @ccall hdl.ssc_data_set_array(data::Ptr{Cvoid}, "wind_turbine_powercurve_powerout"::Cstring,
            powercurve::Ptr{Cdouble}, length(powercurve)::Cint)::Cvoid

        wind_farm_xCoordinates = [Float64(0)]
        @ccall hdl.ssc_data_set_array(data::Ptr{Cvoid}, "wind_farm_xCoordinates"::Cstring,
            wind_farm_xCoordinates::Ptr{Cdouble}, 1::Cint)::Cvoid

        wind_farm_yCoordinates = [Float64(0)]
        @ccall hdl.ssc_data_set_array(data::Ptr{Cvoid}, "wind_farm_yCoordinates"::Cstring,
            wind_farm_yCoordinates::Ptr{Cdouble}, 1::Cint)::Cvoid

        # example of getting a number:
        # val = convert(Cdouble, 0.0)
        # ref = Ref(val)
        # @ccall hdl.ssc_data_get_number(data::Ptr{Cvoid}, "wind_resource_shear"::Cstring, ref::Ptr{Cdouble})::Cvoid
        # Float64(ref[])

        if !Bool(@ccall hdl.ssc_module_exec(wind_module::Ptr{Cvoid}, data::Ptr{Cvoid})::Cint)
            log_type = 0
            log_type_ref = Ref(log_type)
            log_time = 0
            log_time_ref = Ref(log_time)
            msg_ptr = @ccall hdl.ssc_module_log(wind_module::Ptr{Cvoid}, 0::Cint, log_type_ref::Ptr{Cvoid},
                                            log_time_ref::Ptr{Cvoid})::Cstring
            msg = "no message from ssc_module_log."
            try
                msg = unsafe_string(msg_ptr)
            finally
                @error("SAM Wind simulation error: $msg")
            end
        end

        len = 0
        len_ref = Ref(len)
        a = @ccall hdl.ssc_data_get_array(data::Ptr{Cvoid}, "gen"::Cstring, len_ref::Ptr{Cvoid})::Ptr{Float64}
        for i in range(1, stop=8760)
            push!(sam_prodfactor, unsafe_load(a, i))
        end
        @ccall hdl.ssc_module_free(wind_module::Ptr{Cvoid})::Cvoid
        @ccall hdl.ssc_data_free(data::Ptr{Cvoid})::Cvoid

    catch e
        @error "Problem calling SAM C library!"
        showerror(stdout, e)
    end

    if !(length(sam_prodfactor) == 8760)
        @error "Wind production factor from SAM has length $(length(sam_prodfactor)) (should be 8760)."
    end

    @assert !(nothing in sam_prodfactor) "Did not get complete Wind production factor from SAM."

    # normalize by system_capacity
    normalized_prod_factor = sam_prodfactor ./ system_capacity
    # upsample if time steps per hour > 1
    if time_steps_per_hour > 1
        normalized_prod_factor = repeat(normalized_prod_factor, inner=time_steps_per_hour)
    end
    return normalized_prod_factor
end

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite


text:
"""
Copyright Marcus Greiff 2020

Tests of the attitude conversions and distances pertaining to SO(3)
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

@verbose(1, "Testing the SO(3) maps...")
@testset "Attitude maps" begin
    @testset "Test the SO(3) hat map" begin
        @testset "Good inputs" begin
            # Fixture
            A, Ac, Am = get_matrix_and_copy(rand(3,3))
            u, uc, um = get_matrix_and_copy(rand(3,1))
            # Test
            status = ccall((:cont_SO3_hat, AerialVehicleControl.CONT_LIB_PATH),
                Cint,
                (Ref{matrix_double_t}, Ref{matrix_double_t}),
                um,
                Am)
            # Assert
            @test isapprox(status,1);
            @test isapprox(u, uc);
            @test !isapprox(A, Ac);
            @test isapprox(A, skew(u));
        end
        @testset "Bad inputs" begin
            # Fixture
            A, Ac, Am = get_matrix_and_copy(rand(3,4))
            u, uc, um = get_matrix_and_copy(rand(3,1))
            # Test
            status = ccall((:cont_SO3_hat, AerialVehicleControl.CONT_LIB_PATH),
                Cint,
                (Ref{matrix_double_t}, Ref{matrix_double_t}),
                um,
                Am)
            # Assert
            @test isapprox(status, 0);
        end
    end
    @testset "Test the SO(3) vee map" begin
        @testset "Good inputs" begin
            # Fixture
            u               = rand(3,1);
            S, Sc, Sm       = get_matrix_and_copy(skew(u))
            out, outc, outm = get_matrix_and_copy(rand(3,1))
            # Test
            status = ccall((:cont_SO3_vee, AerialVehicleControl.CONT_LIB_PATH),
                Cint,
                (Ref{matrix_double_t}, Ref{matrix_double_t}),
                Sm,
                outm)
            # Assert
            @test isapprox(status,1);
            @test isapprox(S, Sc);
            @test !isapprox(out, outc);
            @test isapprox(u, out);
        end
        @testset "Bad inputs" begin
            # Fixture
            u               = rand(3,1);
            S, Sc, Sm       = get_matrix_and_copy(skew(u))
            out, outc, outm = get_matrix_and_copy(rand(3,2))
            # Test
            status = ccall((:cont_SO3_vee, AerialVehicleControl.CONT_LIB_PATH),
                Cint,
                (Ref{matrix_double_t}, Ref{matrix_double_t}),
                Sm,
                outm)
            # Assert
            @test isapprox(status, 0);
        end
    end
    @testset "Test the SO(3) Psi distance" begin
        @testset "Good inputs" begin
            for i = 1:10
                # Fixture
                R1, R1c, R1m  = get_matrix_and_copy(exp(skew(rand(3,1))))
                R2, R2c, R2m  = get_matrix_and_copy(exp(skew(rand(3,1))))
                # Test
                value = ccall((:cont_SO3_distance, AerialVehicleControl.CONT_LIB_PATH),
                    Cdouble,
                    (Ref{matrix_double_t}, Ref{matrix_double_t}),
                    R1m,
                    R2m)
                # Assert
                @test isapprox(value, tr(I-R1'R2)/2.0);
                @test isapprox(R1, R1c)
                @test isapprox(R2, R2c)
            end
        end
        # TODO write bad input tests
    end
    @testset "Test the SO(3) exponential map" begin
        @testset "Good inputs" begin
            for i = 1:5
                # Fixture
                u, uc, um = get_matrix_and_copy(rand(3,1));
                R, Rc, Rm = get_matrix_and_copy(rand(3,3));
                # Test
                status = ccall((:cont_SO3_Exp, AerialVehicleControl.CONT_LIB_PATH),
                    Cint,
                    (Ref{matrix_double_t}, Ref{matrix_double_t}),
                    um,
                    Rm)
                # Assert
                @test isapprox(status,1);
                @test isapprox(u, uc);
                @test !isapprox(R, Rc);
                @test isapprox(R, exp(skew(uc)));
            end
        end
        # TODO write bad input tests
    end
    @testset "Test the SO(3) logarithmic map" begin
        @testset "Good inputs" begin
            for i = 1:5
                # Fixture
                u, uc, um = get_matrix_and_copy(rand(3,1));
                R, Rc, Rm = get_matrix_and_copy(exp(skew(u)));
                out, outc, outm = get_matrix_and_copy(rand(3,1));
                # Test
                status = ccall((:cont_SO3_Log, AerialVehicleControl.CONT_LIB_PATH),
                    Cint,
                    (Ref{matrix_double_t}, Ref{matrix_double_t}),
                    Rm,
                    outm)

                theta = acos((tr(R) - 1)/2)
                # Assert
                @test isapprox(status,1);
                # We should have changed out, but not R or u
                @test isapprox(R, Rc);
                @test isapprox(u, uc);
                @test !isapprox(out, outc);
                # We chould be able to recover u from R
                @test isapprox(u, out);
            end
        end
        # TODO write bad input tests
    end
end


text:
module TOML

using Dates

include("buffer.jl")
include("value.jl")
include("token.jl")
include("tokenizer.jl")
include("parser.jl")

end # module