text:
du[1]=(-4).*(2+P1.^2+(-2).*2.^(1/2).*P2+P2.^2+2.*P1.*((-1).*2.^(1/2)+P2) +(Q1+(-1).*Q2).^2).^(-2).*(2+P1.^2+2.*2.^(1/2).*P2+P2.^2+2.*P1.*(   2.^(1/2)+P2)+(Q1+(-1).*Q2).^2).^(-2).*(P2.^2+Q2.^2).*(2+P1.^2+2.*   P1.*(2.^(1/2)+(-1).*P2)+(-2).*2.^(1/2).*P2+P2.^2+(Q1+Q2).^2).^(-2)   .*(2+P1.^2+2.*2.^(1/2).*P2+P2.^2+(-2).*P1.*(2.^(1/2)+P2)+(Q1+Q2) .^2).^(-2).*((P1.^2+Q1.^2).*(2+(-2).*2.^(1/2).*(P1+P2)+(P1+P2).^2+  (Q1+(-1).*Q2).^2).^(-1).*(2+2.*2.^(1/2).*(P1+P2)+(P1+P2).^2+(Q1+(   -1).*Q2).^2).^(-1).*(P2.^2+Q2.^2).*(2+2.*2.^(1/2).*(P1+(-1).*P2)+(   P1+(-1).*P2).^2+(Q1+Q2).^2).^(-1).*(2+(P1+(-1).*P2).^2+2.*2.^(1/2)  .*((-1).*P1+P2)+(Q1+Q2).^2).^(-1)).^(-1/2).*(3.*P1.^9+24.*P1.^6.*   P2.*Q1.*Q2+4.*P1.^7.*((-4)+(-2).*P2.^2+3.*Q1.^2+2.*Q2.^2)+8.*   P1.^4.*P2.*Q1.*Q2.*(4+(-2).*P2.^2+7.*Q1.^2+2.*Q2.^2)+4.*P1.^3.*   Q1.^2.*(12+3.*P2.^4+3.*Q1.^4+(-4).*Q2.^2+3.*Q2.^4+2.*Q1.^2.*((-2)+   Q2.^2)+(-2).*P2.^2.*((-2)+Q1.^2+Q2.^2))+8.*P2.*Q1.^3.*Q2.*(P2.^4+(   Q1.^2+(-1).*Q2.^2).^2+2.*P2.^2.*(2+Q1.^2+Q2.^2)+(-4).*(3+Q1.^2+   Q2.^2))+2.*P1.^5.*(12+3.*P2.^4+9.*Q1.^4+(-4).*Q2.^2+3.*Q2.^4+8.*   Q1.^2.*((-2)+Q2.^2)+(-2).*P2.^2.*((-2)+4.*Q1.^2+Q2.^2))+(-8).*   P1.^2.*P2.*Q1.*Q2.*(P2.^4+(-5).*Q1.^4+Q2.^4+(-8).*Q1.^2.*(2+Q2.^2)   +(-4).*(3+Q2.^2)+2.*P2.^2.*(2+4.*Q1.^2+Q2.^2))+P1.*((-1).*P2.^8+   3.*Q1.^8+8.*Q1.^2.*((-2)+Q2.^2).*(2+Q2.^2).^2+(-1).*(2+Q2.^2).^4+(   -4).*P2.^6.*((-2)+2.*Q1.^2+Q2.^2)+(-2).*Q1.^4.*(20+36.*Q2.^2+5.*   Q2.^4)+(-2).*P2.^4.*(12+5.*Q1.^4+(-4).*Q2.^2+3.*Q2.^4+4.*Q1.^2.*((   -2)+Q2.^2))+4.*P2.^2.*((-1).*((-2)+Q2.^2).*(2+Q2.^2).^2+Q1.^4.*(   18+23.*Q2.^2)+2.*Q1.^2.*(4+20.*Q2.^2+Q2.^4))));
du[2]=4.*(P1.^2+Q1.^2).*(2+P1.^2+(-2).*2.^(1/2).*P2+P2.^2+2.*P1.*((-1).*  2.^(1/2)+P2)+(Q1+(-1).*Q2).^2).^(-2).*(2+P1.^2+2.*2.^(1/2).*P2+  P2.^2+2.*P1.*(2.^(1/2)+P2)+(Q1+(-1).*Q2).^2).^(-2).*(2+P1.^2+2.*   P1.*(2.^(1/2)+(-1).*P2)+(-2).*2.^(1/2).*P2+P2.^2+(Q1+Q2).^2).^(-2)   .*(2+P1.^2+2.*2.^(1/2).*P2+P2.^2+(-2).*P1.*(2.^(1/2)+P2)+(Q1+Q2)   .^2).^(-2).*((P1.^2+Q1.^2).*(2+(-2).*2.^(1/2).*(P1+P2)+(P1+P2).^2+   (Q1+(-1).*Q2).^2).^(-1).*(2+2.*2.^(1/2).*(P1+P2)+(P1+P2).^2+(Q1+(   -1).*Q2).^2).^(-1).*(P2.^2+Q2.^2).*(2+2.*2.^(1/2).*(P1+(-1).*P2)+(   P1+(-1).*P2).^2+(Q1+Q2).^2).^(-1).*(2+(P1+(-1).*P2).^2+2.*2.^(1/2)   .*((-1).*P1+P2)+(Q1+Q2).^2).^(-1)).^(-1/2).*((-3).*P2.^9+(-24).*   P1.*P2.^6.*Q1.*Q2+8.*P1.*P2.^4.*Q1.*Q2.*((-4)+2.*P1.^2+(-2).*   Q1.^2+(-7).*Q2.^2)+4.*P2.^7.*(4+2.*P1.^2+(-2).*Q1.^2+(-3).*Q2.^2)+  P2.*((P1.^4+2.*P1.^2.*((-2)+Q1.^2)+(2+Q1.^2).^2).^2+8.*(P1.^6+  P1.^4.*((-2)+Q1.^2)+(-1).*((-2)+Q1.^2).*(2+Q1.^2).^2+(-1).*P1.^2.*  (4+20.*Q1.^2+Q1.^4)).*Q2.^2+2.*(20+5.*P1.^4+36.*Q1.^2+5.*Q1.^4+(-2).*P1.^2.*(18+23.*Q1.^2)).*Q2.^4+(-3).*Q2.^8)+(-4).*P2.^3.*  Q2.^2.*(12+3.*P1.^4+(-4).*Q1.^2+3.*Q1.^4+2.*((-2)+Q1.^2).*Q2.^2+   3.*Q2.^4+(-2).*P1.^2.*((-2)+Q1.^2+Q2.^2))+(-8).*P1.*Q1.*Q2.^3.*(   P1.^4+(Q1.^2+(-1).*Q2.^2).^2+2.*P1.^2.*(2+Q1.^2+Q2.^2)+(-4).*(3+   Q1.^2+Q2.^2))+(-2).*P2.^5.*(12+3.*P1.^4+(-4).*Q1.^2+3.*Q1.^4+8.*((  -2)+Q1.^2).*Q2.^2+9.*Q2.^4+(-2).*P1.^2.*((-2)+Q1.^2+4.*Q2.^2))+8.*  P1.*P2.^2.*Q1.*Q2.*((-12)+P1.^4+(-4).*Q1.^2+Q1.^4+(-8).*(2+Q1.^2)   .*Q2.^2+(-5).*Q2.^4+2.*P1.^2.*(2+Q1.^2+4.*Q2.^2)));
du[3]=4.*(2+P1.^2+(-2).*2.^(1/2).*P2+P2.^2+2.*P1.*((-1).*2.^(1/2)+P2)+(Q1+(-1).*Q2).^2).^(-2).*(2+P1.^2+2.*2.^(1/2).*P2+P2.^2+2.*P1.*( 2.^(1/2)+P2)+(Q1+(-1).*Q2).^2).^(-2).*(P2.^2+Q2.^2).*(2+P1.^2+2.* P1.*(2.^(1/2)+(-1).*P2)+(-2).*2.^(1/2).*P2+P2.^2+(Q1+Q2).^2).^(-2) .*(2+P1.^2+2.*2.^(1/2).*P2+P2.^2+(-2).*P1.*(2.^(1/2)+P2)+(Q1+Q2) .^2).^(-2).*((P1.^2+Q1.^2).*(2+(-2).*2.^(1/2).*(P1+P2)+(P1+P2).^2+(Q1+(-1).*Q2).^2).^(-1).*(2+2.*2.^(1/2).*(P1+P2)+(P1+P2).^2+(Q1+( -1).*Q2).^2).^(-1).*(P2.^2+Q2.^2).*(2+2.*2.^(1/2).*(P1+(-1).*P2)+(P1+(-1).*P2).^2+(Q1+Q2).^2).^(-1).*(2+(P1+(-1).*P2).^2+2.*2.^(1/2).*((-1).*P1+P2)+(Q1+Q2).^2).^(-1)).^(-1/2).*(3.*Q1.^9+24.*P1.*P2.*Q1.^6.*Q2+4.*Q1.^7.*(4+3.*P1.^2+2.*P2.^2+(-2).*Q2.^2)+8.*P1.*P2.* Q1.^4.*Q2.*((-4)+7.*P1.^2+2.*P2.^2+(-2).*Q2.^2)+4.*P1.^2.*Q1.^3.*(12+3.*P1.^4+4.*P2.^2+3.*P2.^4+(-2).*(2+P2.^2).*Q2.^2+3.*Q2.^4+2.*P1.^2.*(2+P2.^2+(-1).*Q2.^2))+2.*Q1.^5.*(12+9.*P1.^4+4.*P2.^2+3.*P2.^4+(-2).*(2+P2.^2).*Q2.^2+3.*Q2.^4+8.*P1.^2.*(2+P2.^2+(-1).* Q2.^2))+8.*P1.*P2.*Q1.^2.*Q2.*(5.*P1.^4+8.*P1.^2.*((-2)+P2.^2+(-1) .*Q2.^2)+4.*(3+(-1).*P2.^2+Q2.^2)+(-1).*(P2.^2+Q2.^2).^2)+8.*P1.^3.*P2.*Q2.*(P1.^4+4.*((-3)+P2.^2+(-1).*Q2.^2)+2.*P1.^2.*(2+( -1).*P2.^2+Q2.^2)+(P2.^2+Q2.^2).^2)+Q1.*(3.*P1.^8+(-1).*(P2.^4+2.* P2.^2.*((-2)+Q2.^2)+(2+Q2.^2).^2).^2+(-2).*P1.^4.*(20+5.*P2.^4+ 36.*Q2.^2+5.*Q2.^4+(-2).*P2.^2.*(18+23.*Q2.^2))+8.*P1.^2.*(P2.^6+ P2.^4.*((-2)+Q2.^2)+(-1).*((-2)+Q2.^2).*(2+Q2.^2).^2+(-1).*P2.^2.*(4+20.*Q2.^2+Q2.^4))));
du[4]=(-4).*(P1.^2+Q1.^2).*(2+P1.^2+(-2).*2.^(1/2).*P2+P2.^2+2.*P1.*((-1).*2.^(1/2)+P2)+(Q1+(-1).*Q2).^2).^(-2).*(2+P1.^2+2.*2.^(1/2).*P2+P2.^2+2.*P1.*(2.^(1/2)+P2)+(Q1+(-1).*Q2).^2).^(-2).*((-8).*P1.*P2.^3.*Q1.*(P1.^4+4.*((-3)+P2.^2+(-1).*Q1.^2)+2.*P1.^2.*(2+(-1).*P2.^2+Q1.^2)+(P2.^2+Q1.^2).^2)+(P1.^8+(-3).*P2.^8+8.*P2.^2.*((-2)+Q1.^2).*(2+Q1.^2).^2+(2+Q1.^2).^4+4.*P1.^6.*((-2)+(-2).*P2.^2+Q1.^2)+2.*P1.^4.*(12+8.*P2.^2+5.*P2.^4+(-4).*(1+P2.^2).*Q1.^2+3.*Q1.^4)+2.*P2.^4.*(20+36.*Q1.^2+5.*Q1.^4)+4.*P1.^2.*(((-2)+Q1.^2).*(2+Q1.^2).^2+(-1).*P2.^4.*(18+23.*Q1.^2)+2.*P2.^2.*(4+20.*Q1.^2+Q1.^4))).*Q2+8.*P1.*P2.*Q1.*(P1.^4+(-5).*P2.^4+Q1.^4+8.*P2.^2.*(2+Q1.^2)+(-4).*(3+Q1.^2)+2.*P1.^2.*(2+(-4).*P2.^2+Q1.^2)).*Q2.^2+( -4).*P2.^2.*(12+3.*P1.^4+4.*P2.^2+3.*P2.^4+(-2).*(2+P2.^2).*Q1.^2+3.*Q1.^4+2.*P1.^2.*(2+P2.^2+(-1).*Q1.^2)).*Q2.^3+(-8).*P1.*P2.* Q1.*(2.*P1.^2+7.*P2.^2+(-2).*(2+Q1.^2)).*Q2.^4+(-2).*(12+3.*P1.^4+ 9.*P2.^4+(-4).*Q1.^2+3.*Q1.^4+P1.^2.*(4+8.*P2.^2+(-2).*Q1.^2)+(-8) .*P2.^2.*((-2)+Q1.^2)).*Q2.^5+(-24).*P1.*P2.*Q1.*Q2.^6+(-4).*(4+ 2.*P1.^2+3.*P2.^2+(-2).*Q1.^2).*Q2.^7+(-3).*Q2.^9).*(2+P1.^2+2.* P1.*(2.^(1/2)+(-1).*P2)+(-2).*2.^(1/2).*P2+P2.^2+(Q1+Q2).^2).^(-2) .*(2+P1.^2+2.*2.^(1/2).*P2+P2.^2+(-2).*P1.*(2.^(1/2)+P2)+(Q1+Q2) .^2).^(-2).*((P1.^2+Q1.^2).*(2+(-2).*2.^(1/2).*(P1+P2)+(P1+P2).^2+ (Q1+(-1).*Q2).^2).^(-1).*(2+2.*2.^(1/2).*(P1+P2)+(P1+P2).^2+(Q1+(-1).*Q2).^2).^(-1).*(P2.^2+Q2.^2).*(2+2.*2.^(1/2).*(P1+(-1).*P2)+( P1+(-1).*P2).^2+(Q1+Q2).^2).^(-1).*(2+(P1+(-1).*P2).^2+2.*2.^(1/2) .*((-1).*P1+P2)+(Q1+Q2).^2).^(-1)).^(-1/2);


text:
module Circuit

import Microwave: AbstractParams, AbstractData
import Base: +, -, *, /, ^, ==, convert, promote_rule, show, getindex

export CircuitParams, Impedance, Admittance, CircuitData, ∥
export capacitor, inductor, resistor, check_frequency_identical

include("CircuitParams.jl")
include("CircuitData.jl")
include("convert.jl")
include("operations.jl")
include("RLC.jl")

end

ext: jl


text:
# ---
# title: 733. Flood Fill
# id: problem733
# author: zhwang
# date: 2022-02-18
# difficulty: Easy
# categories: Depth-first Search
# link: <https://leetcode.com/problems/flood-fill/description/>
# hidden: true
# ---
#
# An `image` is represented by a 2-D array of integers, each integer
# representing the pixel value of the image (from 0 to 65535).
#
# Given a coordinate `(sr, sc)` representing the starting pixel (row and column)
# of the flood fill, and a pixel value `newColor`, "flood fill" the image.
#
# To perform a "flood fill", consider the starting pixel, plus any pixels
# connected 4-directionally to the starting pixel of the same color as the
# starting pixel, plus any pixels connected 4-directionally to those pixels
# (also with the same color as the starting pixel), and so on. Replace the color
# of all of the aforementioned pixels with the newColor.
#
# At the end, return the modified image.
#
# **Example 1:**
#
#
#
#     Input:
#     image = [[1,1,1],[1,1,0],[1,0,1]]
#     sr = 1, sc = 1, newColor = 2
#     Output: [[2,2,2],[2,2,0],[2,0,1]]
#     Explanation:
#     From the center of the image (with position (sr, sc) = (1, 1)), all pixels connected
#     by a path of the same color as the starting pixel are colored with the new color.
#     Note the bottom corner is not colored 2, because it is not 4-directionally connected
#     to the starting pixel.
#
#
# **Note:**
#
# * The length of `image` and `image[0]` will be in the range `[1, 50]`.
# * The given starting pixel will satisfy `0 <= sr < image.length` and `0 <= sc < image[0].length`.
# * The value of each color in `image[i][j]` and `newColor` will be an integer in `[0, 65535]`.
#
#
## @lc code=start
using LeetCode

function flood_fill!(image::Matrix{Int}, sr::Int, sc::Int, new_color::Int)::Matrix{Int}
    ind = CartesianIndex((sr + 1, sc + 1)) ## index correction
    old_color, image[ind] = image[ind], new_color
    old_color == new_color && return image
    inds, queue = CartesianIndices(image), [ind]
    while !isempty(queue)
        ind = popfirst!(queue)
        for I in [(-1, 0), (1, 0), (0, -1), (0, 1)]
            new_I = ind + CartesianIndex(I)
            if new_I ∈ inds && image[new_I] == old_color
                image[new_I] = new_color
                push!(queue, new_I)
            end
        end
    end
    image
end

## @lc code=end

## @lc test=start
@testset "733.flood-fill.jl" begin
    @test flood_fill!([[1,1,1] [1,1,0] [1,0,1]], 1, 1, 1) == [[1,1,1] [1,1,0] [1,0,1]]
    @test flood_fill!([[1,1,1] [1,1,0] [1,0,1]], 1, 1, 2) == [[2,2,2] [2,2,0] [2,0,1]]
    @test flood_fill!([[0,0,0] [0,1,0]], 1, 1, 2) == [[0,0,0] [0,2,0]]
end
## @lc test=end

## @lc info=start
# link: [solution to 733](https://leetcode-cn.com/problems/flood-fill/solution/733-tu-xiang-xuan-ran-by-rex-rs-fod6/)
## @lc info=end


text:
using MLDatasets, Images

x_train, t_train = MNIST.traindata()
x_test, t_test = MNIST.testdata()

println(size(x_train))
println(size(t_train))
println(size(x_test))
println(size(t_test))

Gray.(x_train[:, :, 1]')  # 何故、そこで転置


text:
function elapsed_time_reporting_criterion(report_after::Number)
    previous_time = time()
    while true
        current_time = time()
        if current_time - previous_time >= report_after
            previous_time = time()
            produce(true)
        else
            produce(false)
        end
    end
end

function iterations_reporting_criterion(report_after::Number)
    previous_iteration = 0
    current_iteration  = previous_iteration
    while true
        current_iteration += 1
        if current_iteration - previous_iteration >= report_after
            previous_iteration = current_iteration
            produce(true)
        else
            produce(false)
        end
    end
end


text:
pi = Parameter(3.14, :π)
x = RandomVariable(Normal(0, 1), :x)
y = RandomVariable(Normal(1, 1), :y)
z = RandomVariable(Normal(0, 1), :z)
jd = JointDistribution([x, y], GaussianCopula([1 0; 0 1]))

inputs = [pi, jd, z]

@testset "Inputs" begin
    @testset "sample" begin
        samples = sample(inputs, 10)

        @test size(samples) == (10, 4)
        @test names(samples) == ["π", "x", "y", "z"]
        @test propertynames(samples) == [:π, :x, :y, :z]
    end

    @testset "mean" begin
        means = mean(inputs)
        @test size(means) == (1, 4)
        @test means.π[1] == 3.14
        @test means.x[1] == 0
        @test means.y[1] == 1
        @test means.z[1] == 0
    end

    @testset "names" begin
        @test names(inputs) == [:π, :x, :y, :z]
    end

    @testset "count_rvs" begin
        @test count_rvs(inputs) == 3
    end
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
@kernel function addToLayer!(x, dx, layerIdx)

    i, j = @index(Global, NTuple)
    @inbounds x[i,j,layerIdx] = x[i,j,layerIdx] + dx

end


text:
########################################
## File Name: forward_backward_simulation_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/19
## Description: Test code for src/forward_backward_simulation.jl with
## closed-loop nominal policies.
########################################

using LinearAlgebra
using Random

@testset "Closed-Loop Forward-Backward Simulation Test" begin

# Localization Task
dtc = 0.01;
dto = 0.2;
dtexec = dto*0.4;
Q = Matrix(1.0I, 2, 2);
Cu = Matrix(1.0I, 2, 2);
transPos = TransModel_Pos();
ukfModel = UKFPosRange();
simPosUKF = SimulatePosRangeLocalization2D(transPos,ukfModel,dtc,dto,dtexec,Q,Cu);
gradientPolicy = GradientMultiTargetLocalizationPolicy();
u_param_min_2D = [-2.,-2.];
u_param_max_2D = [2.,2.];

rng = MersenneTwister(1234)
p = PhysPos(0.,zeros(2));
bVec = VecBelMvNormal([BelMvNormal(0.,rand(rng, 2),10*Matrix(1.0I, 2, 2)) for ii = 1:100]);
s = AugState(p,bVec);

UArray = [PosControl(round(dtc*t,digits=5),zeros(2)) for t = 0:Int64(2.0/dtc)-1];

Y,S,S_before_update,UArray = simulateForward(simPosUKF,gradientPolicy,s,u_param_min_2D,u_param_max_2D,UArray,MersenneTwister(1234));

@test length(UArray) == 200
@test isapprox(UArray[1].vel, [-0.704459, -0.709745], rtol=1e-5)
@test isapprox(UArray[end].vel, [-0.526402, -0.850236], rtol=1e-5)

@test length(S_before_update) == 10
@test S_before_update[1].t == 0.19
@test isapprox(S_before_update[1].p.pos, [-0.133847, -0.134851], rtol=1e-5)
@test isapprox(S_before_update[1].b[1].params, [0.590845, 0.766797, 10.0, 0.0, 0.0, 10.0], rtol=1e-5)
@test S_before_update[10].t == 1.99
@test isapprox(S_before_update[10].p.pos, [-1.27097, -1.52526], rtol=1e-5)
@test isapprox(S_before_update[10].b[1].params, [0.0249566, 0.0788592, 10.0039, -2.21546, -2.21546, 9.06684], rtol=1e-5)

CoSArray = simulateBackward(simPosUKF, S, S_before_update, UArray, Y)

@test length(CoSArray) == 201
@test CoSArray[1].t == 0.0
@test isapprox(CoSArray[1].p.pos, [1310.62, 2216.47], rtol=1e-5)
@test isapprox(CoSArray[1].b[1].params, [-19.9427, -30.5194, 10.22, -1.94483, -2.56533, 10.45], rtol=1e-5)

# Manipulation Task
dtc = 0.02;
dto = 0.2;
dtexec = dto*0.8;
QMP = [Matrix(1.0I, 6, 6) zeros(6,5); zeros(5,6) 0.00*Matrix(1.0I, 5, 5)];
RMP = Matrix(1.0I, 9, 9);
CsMP = QMP;
CuMP = Matrix(1.0I, 3, 3);
simMP = SimulateManipulate2D(dtc,dto,dtexec,QMP,RMP,CsMP,CuMP)
pcontrolPolicy = ManipulatePositionControlPolicy();

xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.2);
uMP = MControl2D(0.,6.,7.,8.);
bMP = BelMvNormal(xMP.t,vec(xMP),diagm([0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]));
u_param_min_3D = [-3.,-3.,-3.];
u_param_max_3D = [3.,3.,3.];
posGain = 1.0;
rotGain = 0.2;

UArray = [MControl2D(round(dtc*t,digits=5),0.,0.0,0.0) for t = 0:Int64(2.0/dtc)-1];
Y,S,S_before_update,UArray = simulateForward(simMP,pcontrolPolicy,bMP,u_param_min_3D,u_param_max_3D,posGain,rotGain,UArray,MersenneTwister(123));

@test length(UArray) == 100
@test isapprox(vec(UArray[1]), [-0.1, -0.2, 0.528319], rtol=1e-5)
@test isapprox(vec(UArray[end]), [-0.253098, -0.637423, -0.598584], rtol=1e-5)

@test length(S_before_update) == 10
@test S_before_update[1].t == 0.18
@test isapprox(S_before_update[1].params[1:11], [0.153234, 0.27066, 0.611928, 0.290451, 0.383289, 0.651479, 3.0, 5.0, 10.0, 11.0, 0.2], rtol=1e-5)
@test S_before_update[10].t == 1.98
@test isapprox(S_before_update[10].params[1:11], [0.300082, 0.74635, 0.478632, 0.250788, 0.579737, 3.4435, 2.99207, 5.34697, 10.5676, 11.1221, 0.205542], rtol=1e-5)

CoSArray = simulateBackward(simMP,S,S_before_update,UArray,Y)

@test length(CoSArray) == 101
@test CoSArray[1].t == 0.0
@test isapprox(CoSArray[1].params[1:11], [0.884668, -0.817838, -16.6739, 1.02655, 0.0991195, -6.46987, -0.171831, -12.4607, 4.27364, 3.02749, -0.656432], rtol=1e-5)

end


text:
function boxcox(z, λ)
    reduce(&, z .> 0) || throw(DomainError("boxcox requires strictly positive data"))
	if λ == 0
		return log.(z)
	else
		return (z.^λ .- 1) ./ λ
	end
end

function inv_boxcox(z, λ)
    if λ == 0
        return exp.(z)
    else
        return (λ .* z .+ 1).^(1/λ)
    end
end

function guerrero(z; lower=-1, upper=2, s=2)
    s >= 2 || throw(ArgumentError("guerrero requires seasonality to be higher or equal to 2"))
    N = length(z)
    rz = reshape(view(z, (N % s + 1):N), N ÷ s, s)
    m = mean(rz, dims=1)
    sd = std(rz, mean=m, dims=1)
    function cv(λ)
        w = sd ./ (m.^(1 - λ))
        mw = mean(w)
        return std(w, mean=mw) / mw
    end
    return optimize(cv, lower, upper).minimizer
end

function difference(z::AbstractVector{T}; d::Integer=1,s::Integer=1) where T
    N = length(z)
    differenced = Vector{T}(undef, N)
    @inbounds @simd for i in 1:N
        differenced[i] = z[i]
    end
    @inbounds for j in 1:d
        for i in 1:N - j * s
            differenced[N - i + 1] -= differenced[N - i + 1 - s]
        end
    end
    return differenced[1 + d * s:N]
end

function integrate(z::AbstractVector{T}, z0::AbstractVector; d::Integer=1,s::Integer=1) where T
    N = length(z)
    integrated = Vector{T}(undef, N + d * s)
    N0 = length(z0)
    for i in 1:N + d * s
        integrated[i] = i > d * s ? z[i - d * s] : i > N0 ? zero(T) : z0[i]
    end
    for j in 1:d
        for i in (1 + s):N + d * s
            integrated[i] += integrated[i - s]
        end
    end
    return integrated
end

integrate(z::AbstractVector{T}; d::Integer=1,s::Integer=1) where T = integrate(z, T[], d=d, s=s)


text:
module Benchmark
  using MD5
  using BenchmarkTools
  using PyCall
  using ..Sources
  using ..Converting

  function test_accuracy(gold_data_source::AbstractSource, test_data_source::AbstractSource)
    gold_data_sentences = extract_sentences(gold_data_source)
    test_data_sentences = extract_sentences(test_data_source)

    total_labels = 0
    correct_labels = 0
    correct_vertices = 0

    for gold_data_sentence_pair in gold_data_sentences
      gold_sentence = gold_data_sentence_pair.second

      if !haskey(test_data_sentences, gold_data_sentence_pair.first)
        continue;
      end

      parsed_gold_sentence = split(gold_sentence, "\n")
      parsed_test_sentence =
        split(test_data_sentences[gold_data_sentence_pair.first], "\n")

      for (index, gold_word_data) in enumerate(parsed_gold_sentence)
        total_labels += 1

        parsed_gold_word_data = split(gold_word_data)
        parsed_test_word_data = split(parsed_test_sentence[index])

        if (parsed_gold_word_data[2] == parsed_test_word_data[2])
          correct_vertices += 1

          if lowercase(parsed_gold_word_data[1]) == lowercase(parsed_test_word_data[1])
            correct_labels += 1
          end
        end
      end
    end

    uas = correct_vertices / total_labels
    las = correct_labels / total_labels

    @info "UAS - $(uas)"
    @info "LAS - $(las)"

    [uas, las]
  end

  function extract_sentences(source)
    sentences = Converting.convert(source) |> (text -> split(text, r"SENTENCE_END")) .|> strip
    sentences_dict = Dict()

    foreach(sentences) do sentence
      sentences_dict[sentence_id(sentence)] = sentence
    end

    sentences_dict
  end

  function sentence_id(sentence)
    @info sentence

    split(sentence, "\n") |>
      (words -> map(word -> split(word)[3], words)) |>
      join |>
      MD5.md5 |>
      join
  end

  function time_measure_stanza()
    #Stanza
    b = @benchmarkable stanza_call() samples=10 seconds=260
    run(b)
  end

  function stanza_call()
      run(`python3 /mnt/f/ed_soft/parser_materials/parsers_results/stanza_scrypt.py`)
  end

  function spacy_call()
    run(`python3 /mnt/f/ed_soft/parser_materials/parsers_results/.py`)
  end

  function time_measure_spacy(input_file, model_path)
    #Spacy
    spacy = pyimport("spacy")
    nlp = spacy.load(model_path)
    text = open(f -> read(f, String), input_file)
    b = @benchmarkable nlp_call($nlp, $text) samples=10 seconds=260
    run(b)
  end

  function nlp_call(nlp, text)
      nlp(text)
  end

  function time_measure_coreNLP()
    # CoreNLP
    b = @benchmarkable corenlp_call() samples=10 seconds=260
    run(b)
  end

  function corenlp_call()
    cd("/Users/admin/education/materials/stanford-corenlp-4.4.0")
    run(`java -cp "*" edu.stanford.nlp.parser.nndep.DependencyParser -model edu/stanford/nlp/models/parser/nndep/english_UD.gz -textFile ../test_texts/test_text.txt -outFile corenlp_parsed_result.txt`)
  end
end


text:
mutable struct MockCartesianFunction
    called_grideval::Bool
end
function ApproxTools.grideval(f::MockCartesianFunction, x::NTuple{N,Any}) where {N}
    f.called_grideval = true
    return broadcast(*,ApproxTools.reshape4grideval(x)...)
end

@testset "grideval" begin
    a = [1,2]
    b = [1.,2.]

    @inferred grideval(*,  1 )
    @inferred grideval(*,  a )
    @inferred grideval(*, (a,))
    @test typeof(grideval(*,  1  )) == Int
    @test eltype(grideval(*,  a  )) == Int
    @test eltype(grideval(*, (a,))) == Int
    @test grideval(*,  1  ) == 1
    @test collect(grideval(*,  a  )) == [1,2]
    @test collect(grideval(*, (a,))) == [1,2]

    @inferred grideval(*,  1,1.0 )
    @inferred grideval((x,y)->x*y,  b )
    @inferred grideval(*,  1,b )
    @inferred grideval(*,  a,b )
    @inferred grideval(*, (a,b))
    @test typeof(grideval(*,  1,1.0)) == Float64
    @test eltype(grideval((x,y)->x*y, b)) == Float64
    @test eltype(grideval(*,  1,b )) == Float64
    @test eltype(grideval(*,  a,b )) == Float64
    @test eltype(grideval(*, (a,b))) == Float64
    @test grideval(*,  1,1.0 ) == 1.0
    @test collect(grideval(*,  1,b )) == [1. 2.]
    @test collect(grideval((x,y)->x*y,  b )) == [1. 2.; 2. 4.]
    @test collect(grideval(*,  a,b )) == [1. 2.; 2. 4.]
    @test collect(grideval(*, (a,b))) == [1. 2.; 2. 4.]

    m = MockCartesianFunction(false)
    f = let m = m; @gridfun( x -> x + $m(x) ); end
    @test @inferred(grideval(f, a)) ≈ 2*a
    @test m.called_grideval

    m = MockCartesianFunction(false)
    f = let m = m; @gridfun( (x,y) -> x + $m(x,y) ); end
    @test @inferred(grideval(f, a,b)) ≈ a .+ a .*b'
    @test m.called_grideval
end


text:
struct Entry
    title::String
    source_location::LineNumberNode
    code::Expr
    func::Function
    used_functions::Set{Symbol}
end

function Entry(title, source_location, code, func)
    used_functions = Set{Symbol}()
    MacroTools.postwalk(code) do x
        if @capture(x, f_(xs__))
            push!(used_functions, Symbol(string(f)))
        end
        if @capture(x, f_(xs__) do; body__; end)
            push!(used_functions, Symbol(string(f)))
        end
        return x
    end
    return Entry(title, source_location, code, func, used_functions)
end

function unique_name(entry::Entry)
    loc = entry.source_location
    filename = splitext(basename(string(loc.file)))[1]
    sep = isempty(entry.title) ? "" : "_"
    return replace(string(filename, "_", loc.line, sep, entry.title), " " => "")
end

function nice_title(entry::Entry)
    isempty(entry.title) || return entry.title
    return unique_name(entry)
end

const DATABASE = Dict{String, Entry}()

function cell_expr(name, code, source)
    key = string(source.file, ":", source.line)
    return quote
        closure = () -> $(esc(code))
        entry = ReferenceTests.Entry(
            $(string(name)),
            $(QuoteNode(source)),
            $(QuoteNode(code)),
            closure
        )
        ReferenceTests.DATABASE[$key] = entry
        # display(closure())
    end
end

macro cell(name, code)
    return cell_expr(name, code, __source__)
end

macro cell(code)
    return cell_expr("", code, __source__)
end

"""
    save_result(path, object)

Helper, to more easily save all kind of results from the test database
"""
function save_result(path::String, scene::Makie.FigureLike)
    FileIO.save(path * ".png", scene)
end

function save_result(path::String, stream::VideoStream)
    FileIO.save(path * ".mp4", stream)
end

function save_result(path::String, object)
    FileIO.save(path, object)
end

function load_database()
    empty!(DATABASE)
    include(joinpath(@__DIR__, "tests/text.jl"))
    include(joinpath(@__DIR__, "tests/attributes.jl"))
    include(joinpath(@__DIR__, "tests/examples2d.jl"))
    include(joinpath(@__DIR__, "tests/examples3d.jl"))
    include(joinpath(@__DIR__, "tests/short_tests.jl"))
    return DATABASE
end


text:
#
# BlockSparse storage
#

struct BlockSparse{ElT,VecT,N} <: TensorStorage{ElT}
  data::VecT
  blockoffsets::BlockOffsets{N}  # Block number-offset pairs
  function BlockSparse(
    data::VecT, blockoffsets::BlockOffsets{N}
  ) where {VecT<:AbstractVector{ElT},N} where {ElT}
    return new{ElT,VecT,N}(data, blockoffsets)
  end
end

function BlockSparse(
  ::Type{ElT}, blockoffsets::BlockOffsets, dim::Integer; vargs...
) where {ElT<:Number}
  return BlockSparse(zeros(ElT, dim), blockoffsets; vargs...)
end

function BlockSparse(x::Number, blockoffsets::BlockOffsets, dim::Integer; vargs...)
  return BlockSparse(fill(x, dim), blockoffsets; vargs...)
end

function BlockSparse(
  ::Type{ElT}, ::UndefInitializer, blockoffsets::BlockOffsets, dim::Integer; vargs...
) where {ElT<:Number}
  return BlockSparse(Vector{ElT}(undef, dim), blockoffsets; vargs...)
end

function BlockSparse(blockoffsets::BlockOffsets, dim::Integer; vargs...)
  return BlockSparse(Float64, blockoffsets, dim; vargs...)
end

function BlockSparse(::UndefInitializer, blockoffsets::BlockOffsets, dim::Integer; vargs...)
  return BlockSparse(Float64, undef, blockoffsets, dim; vargs...)
end

copy(D::BlockSparse) = BlockSparse(copy(data(D)), copy(blockoffsets(D)))

setdata(B::BlockSparse, ndata) = BlockSparse(ndata, blockoffsets(B))

#
# Random
#

function randn(
  ::Type{<:BlockSparse{ElT}}, blockoffsets::BlockOffsets, dim::Integer
) where {ElT<:Number}
  return BlockSparse(randn(ElT, dim), blockoffsets)
end

#function BlockSparse{ElR}(data::VecT,offsets) where {ElR,VecT<:AbstractVector{ElT}} where {ElT}
#  ElT == ElR ? BlockSparse(data,offsets) : BlockSparse(ElR.(data),offsets)
#end
#BlockSparse{ElT}() where {ElT} = BlockSparse(ElT[],BlockOffsets())

datatype(::Type{<:BlockSparse{<:Any,DataT}}) where {DataT} = DataT

similar(D::BlockSparse) = setdata(D, similar(data(D)))

# TODO: test this function
similar(D::BlockSparse, ::Type{ElT}) where {ElT} = setdata(D, similar(data(D), ElT))

function similartype(::Type{StoreT}, ::Type{ElT}) where {StoreT<:BlockSparse,ElT}
  return BlockSparse{ElT,similartype(datatype(StoreT), ElT),ndims(StoreT)}
end

# TODO: check the offsets are the same?
function copyto!(D1::BlockSparse, D2::BlockSparse)
  blockoffsets(D1) ≠ blockoffsets(D1) &&
    error("Cannot copy between BlockSparse storages with different offsets")
  copyto!(data(D1), data(D2))
  return D1
end

Base.real(::Type{BlockSparse{T}}) where {T} = BlockSparse{real(T)}

complex(::Type{BlockSparse{T}}) where {T} = BlockSparse{complex(T)}

ndims(::BlockSparse{T,V,N}) where {T,V,N} = N

eltype(::BlockSparse{T}) where {T} = eltype(T)
# This is necessary since for some reason inference doesn't work
# with the more general definition (eltype(Nothing) === Any)
eltype(::BlockSparse{Nothing}) = Nothing
eltype(::Type{BlockSparse{T}}) where {T} = eltype(T)

dense(::Type{<:BlockSparse{ElT,VecT}}) where {ElT,VecT} = Dense{ElT,VecT}

function promote_rule(
  ::Type{<:BlockSparse{ElT1,VecT1,N}}, ::Type{<:BlockSparse{ElT2,VecT2,N}}
) where {ElT1,ElT2,VecT1,VecT2,N}
  return BlockSparse{promote_type(ElT1, ElT2),promote_type(VecT1, VecT2),N}
end

function promote_rule(
  ::Type{<:BlockSparse{ElT1,VecT1,N1}}, ::Type{<:BlockSparse{ElT2,VecT2,N2}}
) where {ElT1,ElT2,VecT1,VecT2,N1,N2}
  return BlockSparse{promote_type(ElT1, ElT2),promote_type(VecT1, VecT2),NR} where {NR}
end

function promote_rule(
  ::Type{<:BlockSparse{ElT1,Vector{ElT1},N1}}, ::Type{ElT2}
) where {ElT1,ElT2<:Number,N1}
  ElR = promote_type(ElT1, ElT2)
  VecR = Vector{ElR}
  return BlockSparse{ElR,VecR,N1}
end

function convert(
  ::Type{<:BlockSparse{ElR,VecR,N}}, D::BlockSparse{ElD,VecD,N}
) where {ElR,VecR,N,ElD,VecD}
  return setdata(D, convert(VecR, data(D)))
end

"""
isblocknz(T::BlockSparse,
          block::Block)

Check if the specified block is non-zero.
"""
function isblocknz(T::BlockSparse{ElT,VecT,N}, block::Block{N}) where {ElT,VecT,N}
  return isassigned(blockoffsets(T), block)
end

# If block is input as Tuple
isblocknz(T::BlockSparse, block) = isblocknz(T, Block(block))

# Given a specified block, return a Dense storage that is a view to the data
# in that block. Return nothing if the block is structurally zero
function blockview(T::BlockSparse, block)
  #error("Block must be structurally non-zero to get a view")
  !isblocknz(T, block) && return nothing
  blockoffsetT = offset(T, block)
  blockdimT = blockdim(T, block)
  dataTslice = @view data(T)[(blockoffsetT + 1):(blockoffsetT + blockdimT)]
  return Dense(dataTslice)
end

# XXX this is not well defined with new Dictionary design
#function (D1::BlockSparse + D2::BlockSparse)
#  # This could be of order nnzblocks, avoid?
#  if blockoffsets(D1) == blockoffsets(D2)
#    return BlockSparse(data(D1)+data(D2),blockoffsets(D1))
#  end
#  blockoffsetsR,nnzR = union(blockoffsets(D1),nnz(D1),
#                             blockoffsets(D2),nnz(D2))
#  R = BlockSparse(undef,blockoffsetsR,nnzR)
#  for (blockR,offsetR) in blockoffsets(R)
#    blockview1 = blockview(D1,blockR)
#    blockview2 = blockview(D2,blockR)
#    blockviewR = blockview(R,blockR)
#    if isnothing(blockview1)
#      copyto!(blockviewR,blockview2)
#    elseif isnothing(blockview2)
#      copyto!(blockviewR,blockview1)
#    else
#      # TODO: is this fast?
#      blockviewR .= blockview1 .+ blockview2
#    end
#  end
#  return R
#end

# Helper function for HDF5 write/read of BlockSparse
function offsets_to_array(boff::BlockOffsets{N}) where {N}
  nblocks = length(boff)
  asize = (N + 1) * nblocks
  n = 1
  a = Vector{Int}(undef, asize)
  for bo in pairs(boff)
    for j in 1:N
      a[n] = bo[1][j]
      n += 1
    end
    a[n] = bo[2]
    n += 1
  end
  return a
end

# Helper function for HDF5 write/read of BlockSparse
function array_to_offsets(a, N::Int)
  asize = length(a)
  nblocks = div(asize, N + 1)
  boff = BlockOffsets{N}()
  j = 0
  for b in 1:nblocks
    insert!(boff, Block(ntuple(i -> (a[j + i]), N)), a[j + N + 1])
    j += (N + 1)
  end
  return boff
end

function HDF5.write(parent::Union{HDF5.File,HDF5.Group}, name::String, B::BlockSparse)
  g = create_group(parent, name)
  attributes(g)["type"] = "BlockSparse{$(eltype(B))}"
  attributes(g)["version"] = 1
  if eltype(B) != Nothing
    write(g, "ndims", ndims(B))
    write(g, "data", data(B))
    off_array = offsets_to_array(blockoffsets(B))
    write(g, "offsets", off_array)
  end
end

function HDF5.read(
  parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{Store}
) where {Store<:BlockSparse}
  g = open_group(parent, name)
  ElT = eltype(Store)
  typestr = "BlockSparse{$ElT}"
  if read(attributes(g)["type"]) != typestr
    error("HDF5 group or file does not contain $typestr data")
  end
  N = read(g, "ndims")
  off_array = read(g, "offsets")
  boff = array_to_offsets(off_array, N)
  # Attribute __complex__ is attached to the "data" dataset
  # by the h5 library used by C++ version of ITensor:
  if haskey(attributes(g["data"]), "__complex__")
    M = read(g, "data")
    nelt = size(M, 1) * size(M, 2)
    data = Vector(reinterpret(ComplexF64, reshape(M, nelt)))
  else
    data = read(g, "data")
  end
  return BlockSparse(data, boff)
end


text:
using Distributions, LinearAlgebra
using UnicodePlots

# Print out the thread count
@info "" Threads.nthreads()

# Utility function
u(X, γ, β) = (1 - exp(-γ * X'β)) / γ

# Single draw of X, gamma, beta, and u
function single_draw(n)
    X = rand(MvNormal(ones(n), 0.25*I))
    gamma = rand(Exponential(0.05))
    beta = rand(MvNormal(ones(n), 0.5*I))
    return u(X, gamma, beta)
end

# Serial simulation
function simulate(n, n_draws)
    draws = zeros(n_draws)

    for i in 1:n_draws
        draws[i] = single_draw(n)
    end

    return draws
end

# Simulation using a lock
function simulate_threads(n, n_draws)
    draws = zeros(n_draws)

    Threads.@threads for i in 1:n_draws
        draw = single_draw(n)
        draws[i] = draw
        # lock(drawlock) do
        #     draws[i] = draw
        # end
    end

    return draws
end

n = 1000
n_draws = 100_000
@time draws1 = simulate(n, n_draws)
@time draws2 = simulate_threads(n, n_draws)

println("Howdy")

histogram(draws1) |> display
histogram(draws2) |> display


text:
# create an index for id6
# building an index is too time consuming
import Base.ht_keyindex
import Base.ht_keyindex2

function buildindex{T}(val::Vector{T})
  index = Dict{T,Vector{Int64}}()
  for (i,v) in enumerate(val)
    dindex = ht_keyindex(index, v)
    if dindex > 0
      push!(index.vals[dindex], i)
    else
      index[v] = [i]
    end
  end
  return index
end

function sumby_index{T,S}(index::Dict{T, Vector{Int64}}, val::Vector{S})
  return Dict{T,S}(kv[1] => sum(val[kv[2]]) for kv in index)
end

@time iid6 = buildindex(id6);  #259,290
@time sumby_index(iid6, v1); # 30.93
@time sumby_index(iid6, v3); #26.7

@time iid4 = buildindex(id4); #9.96
@time sumby_index(iid4, v1); #8.88
@time sumby_index(iid4, v1); #8.63

using SortingAlgorithms
function buildindex2{T}(val::Vector{T})
  index = Dict{T,Vector{Int64}}()
  for (i,v) in enumerate(val)
    dindex = ht_keyindex2(index, v)
    if dindex > 0
      push!(index.vals[dindex], i)
    else
      Base._setindex!(index, [i], v, -dindex)
    end
  end
  return index
end

@time iid42 = buildindex2(id4); #10.94


text:
#########################################################################
#    testing script for distributions and samplers
#########################################################################

using Distributions # used to provide exact cdf of distributions for testing

#########################################################################
#       testing functions and definitions
#########################################################################

N = 10000  # number of steps in MCMC for testing
KSTHRESHOLD = 1.358  #  5% level confidence for Kolmogorov–Smirnov test
KSTHRESHOLD = 10
# TODO : understand why KS is so bad for most samplers * distribs
#         forcing such a high threshold to pass tests

function ksValue(x, distrib)
	global xs = sort(x)
	y = eval( :(cdf($distrib , xs)))
	dn = max( abs([1:length(x)] / length(x) - y) )
	sqrt(length(x))*dn
end

function ksTest(ex::Expr)
	mex = Expr(:block, :(x ~ $ex))
	distrib = Expr(:call,
					Expr(:., :Distributions, Expr(:quote, ex.args[1])),
					ex.args[2:end]...)
	exactMean = eval( :(mean($distrib)) )
	exactMean = isfinite(exactMean) ? exactMean : 1.0  # some distribs (Cauchy) have no defined mean
	exactStd = eval( :(std($distrib)) )
	exactStd = isfinite(exactStd) ? exactStd : 1.0

	for (k,v) in {	"RWM" => RWM(exactStd),
					"HMC" => HMC(2, exactStd/5),
					"MALA" => MALA(exactStd),
					"NUTS" => NUTS()}  # TODO : add other samplers
		print("testing $k sampler on $ex   -")
		srand(1)
		res = run(model(mex, gradient=true, x=exactMean), v, SerialMC((1000:N) ))
		ksv = ksValue(res.samples["x"], distrib)
		println(" KS measure = $ksv")
		@assert ksv < KSTHRESHOLD "correct distrib hyp. rejected"
	end

end

#########################################################################
#  Tests
#########################################################################

ksTest(:(Normal(1, 1)))
ksTest(:(Normal(3, 12)))

ksTest(:(Weibull(1, 1)))
ksTest(:(Weibull(3, 1)))

ksTest(:(Uniform(0, 2)))

ksTest(:(TDist(2.2)))  # very long on NUTS if df <= 2  (due to infinite variance ?)
ksTest(:(TDist(4)))

ksTest(:(Beta(1,2)))
ksTest(:(Beta(3,2)))

ksTest(:(Gamma(1,2)))
ksTest(:(Gamma(3,0.2)))

ksTest(:(Cauchy(0,1)))
ksTest(:(Cauchy(-1,0.2)))

ksTest(:(Exponential(3)))
ksTest(:(Exponential(0.2)))

ksTest(:(LogNormal(-1, 1)))
ksTest(:(LogNormal(2, 0.1)))

# TODO : find a way to test discrete distributions : Bernoulli, Binomial, Poisson


text:
ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)

facts("[sdp] Test simple SDP") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, X[1:3,1:3], SDP)
    @SDconstraint(m, X <= 1/2*eye(3,3))
    @variable(m, Y[1:5,1:5], Symmetric)
    @SDconstraint(m, -ones(5,5) <= Y)
    @SDconstraint(m, Y <= 2*ones(5,5))
    @variable(m, Z[1:4,1:4], Symmetric)
    @SDconstraint(m, ones(4,4) >= Z)

    @constraint(m, trace(X) == 1)
    @constraint(m, trace(Y) == 3)
    @constraint(m, trace(Z) == -1)
    @objective(m, Max, X[1,2] + Y[1,2] + Z[1,2])
    solve(m)

    XX, YY, ZZ, = getvalue(X), getvalue(Y), getvalue(Z)
    Xtrue = [0.25 0.25 0
             0.25 0.25 0
             0    0    0.5]
    Ytrue = fill(0.6, 5, 5)
    Ztrue = [-1.5  3.5 1 1
              3.5 -1.5 1 1
              1    1   1 1
              1    1   1 1]

    @fact norm(XX-Xtrue) --> roughly(0, 1e-2)
    @fact norm(YY-Ytrue) --> roughly(0, 1e-2)
    @fact norm(ZZ-Ztrue) --> roughly(0, 1e-2)
    @fact ispsd(XX) --> true
    @fact ispsd(1/2*eye(3,3)-XX) --> true
    @fact ispsd(YY+ones(5,5)) --> true
    @fact ispsd(2*ones(5,5)-YY) --> true
    @fact ispsd(ones(4,4)-ZZ) --> true
    @fact trace(XX) --> roughly( 1, 1e-4)
    @fact trace(YY) --> roughly( 3, 1e-4)
    @fact trace(ZZ) --> roughly(-1, 1e-4)
    @fact getobjectivevalue(m) --> roughly(4.35, 1e-2)

    @objective(m, Min, X[1,2] + Y[1,2] + Z[1,2])
    solve(m)

    XX, YY, ZZ, = getvalue(X), getvalue(Y), getvalue(Z)
    Xtrue = [ 0.25 -0.25 0
             -0.25  0.25 0
              0     0    0.5]
    Ytrue = fill(0.6, 5, 5)
    Ztrue = [-1.5 -1.5 1 1
             -1.5 -1.5 1 1
              1    1   1 1
              1    1   1 1]

    @fact norm(XX-Xtrue) --> roughly(0, 1e-2)
    @fact norm(YY-Ytrue) --> roughly(0, 1e-2)
    @fact norm(ZZ-Ztrue) --> roughly(0, 1e-2)
    @fact ispsd(XX) --> true
    @fact ispsd(1/2*eye(3,3)-XX) --> true
    @fact ispsd(YY+ones(5,5)) --> true
    @fact ispsd(2*ones(5,5)-YY) --> true
    @fact ispsd(ones(4,4)-ZZ) --> true
    @fact trace(XX) --> roughly( 1, 1e-4)
    @fact trace(YY) --> roughly( 3, 1e-4)
    @fact trace(ZZ) --> roughly(-1, 1e-4)

    # Test SDP constraints
    m = Model(solver=solver)
    @variable(m, X[1:3,1:3], SDP)

    @SDconstraint(m, ones(3,3) <= X)
    @objective(m, Min, trace(ones(3,3)*X))
    stat = solve(m)
    @fact stat --> :Optimal
    XX = getvalue(X)
    @fact ispsd(XX) --> true
    @fact ispsd(XX - ones(3,3)) --> true
    @fact getobjectivevalue(m) --> roughly(9,1e-4)

    # Another test SDP
    m = Model(solver=solver)
    @variable(m, X[1:3,1:3], SDP)
    @variable(m, Y[1:2,1:2], SDP)

    C = eye(3,3)
    A1 = zeros(3,3)
    A1[1,1] = 1.0
    A2 = zeros(3,3)
    A2[2,2] = 1.0
    A3 = zeros(3,3)
    A3[3,3] = 1.0
    D = eye(2,2)
    B1 = ones(2,2)
    B2 = zeros(2,2)
    B2[1,1] = 1
    B3 = zeros(2,2)
    B3[1,2] = B3[2,1] = 2

    @objective(m, Min, trace(C*X)+1+trace(D*Y))
    @constraint(m, trace(A1*X-eye(3,3)/3) == 0)
    @constraint(m, 2*trace(A2*X) == 1)
    @constraint(m, trace(A3*X) >= 2)
    @constraint(m, trace(B1*Y) == 1)
    @constraint(m, trace(B2*Y) == 0)
    @constraint(m, trace(B3*Y) <= 0)
    @constraint(m, trace(A1*X)+trace(B1*Y) >= 1)
    @constraint(m, Y[2,2] == 1)

    stat = solve(m)
    @fact stat --> :Optimal
    XX, YY = getvalue(X), getvalue(Y)
    @fact trace(A1*XX-eye(3,3)/3) --> roughly(0, 1e-5)
    @fact 2*trace(A2*XX) --> roughly(1, 1e-5)
    @fact (trace(A3*XX) >= 2 - 1e-5) --> true
    @fact trace(B1*YY) --> roughly(1, 1e-5)
    @fact trace(B2*YY) --> roughly(0, 1e-5)
    @fact (trace(B3*YY) <= 1e-3) --> true
    @fact (trace(A1*XX)+trace(B1*YY) >= 1) --> true
    @fact YY[2,2] --> roughly(1, 1e-5)
    @fact norm(XX - diagm([1,.5,2])) --> roughly(0, 1e-3)
    @fact norm(YY - [0 0;0 1]) --> roughly(0, 1e-3)
end; end; end

facts("[sdp] Nonsensical SDPs") do
    m = Model()
    @fact_throws @variable(m, unequal[1:5,1:6], SDP)
    # Some of these errors happen at compile time, so we can't use @fact_throws
    @fact macroexpand(:(@variable(m, notone[1:5,2:6], SDP))).head --> :error
    @fact macroexpand(:(@variable(m, oneD[1:5], SDP))).head --> :error
    @fact macroexpand(:(@variable(m, threeD[1:5,1:5,1:5], SDP))).head --> :error
    @fact macroexpand(:(@variable(m, psd[2] <= rand(2,2), SDP))).head --> :error
    @fact macroexpand(:(@variable(m, -ones(3,4) <= foo[1:4,1:4] <= ones(4,4), SDP))).head --> :error
    @fact macroexpand(:(@variable(m, -ones(3,4) <= foo[1:4,1:4] <= ones(4,4), Symmetric))).head --> :error
    @fact macroexpand(:(@variable(m, -ones(4,4) <= foo[1:4,1:4] <= ones(4,5), Symmetric))).head --> :error
    @fact macroexpand(:(@variable(m, -rand(5,5) <= nonsymmetric[1:5,1:5] <= rand(5,5), Symmetric))).head --> :error
end

facts("[sdp] SDP with quadratics") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, X[1:2,1:2], SDP)
    @variable(m, y[0:2])
    @constraint(m, y[0] >= 0)
    @constraint(m, y[1]^2 + y[2]^2 <= y[0]^2)
    @SDconstraint(m, X <= eye(2))
    @constraint(m, X[1,1] + X[1,2] == y[1] + y[2])
    @objective(m, Max, trace(X) - y[0])
    stat = solve(m)

    @fact stat --> :Optimal
    XX, yy = getvalue(X), getvalue(y)
    @fact ispsd(XX) --> true
    @fact (yy[0] >= 0) --> true
    @fact (yy[1]^2 + yy[2]^2 <= yy[0]^2 + 1e-4) --> true
    @fact ispsd(eye(2)-XX) --> true
    @fact ((XX[1,1] + XX[1,2]) - (yy[1] + yy[2])) --> roughly(0,1e-4)
    @fact norm(XX - eye(2)) --> roughly(0, 1e-4)
    @fact norm(yy[:] - [1/sqrt(2), 0.5, 0.5]) --> roughly(0, 1e-4)
    @fact getobjectivevalue(m) --> roughly(1.293, 1e-2)
end; end; end

facts("[sdp] Trivial symmetry constraints are removed (#766)") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    q = 2
    m = 3
    angles1 = linspace(3*pi/4, pi, m)
    angles2 = linspace(0, -pi/2, m)
    V = [3.*cos(angles1)' 1.5.*cos(angles2)';
         3.*sin(angles1)' 1.5.*sin(angles2)']
    p = 2*m
    n = 100

    mod = Model(solver=solver)
    @variable(mod, x[j=1:p] >= 1, Int)
    @variable(mod, u[i=1:q] >= 1)
    @objective(mod, Min, sum(u))
    @constraint(mod, sum(x) <= n)
    for i=1:q
        @SDconstraint(mod, [V*diagm(x./n)*V' eye(q)[:,i] ; eye(q)[i:i,:] u[i]] >= 0)
    end
    f, A, b, var_cones, con_cones = JuMP.conicdata(mod)
    @fact length(f) --> 8
    @fact size(A) --> (23,8)
    @fact length(b) --> 23
    @fact var_cones --> [(:Free,[1,2,3,4,5,6,7,8])]
    @fact con_cones --> [(:NonNeg,[1]),(:NonPos,[2,3,4,5,6,7,8,9]),(:SDP,10:15),(:Zero,16:16),(:SDP,17:22),(:Zero,23:23)]
end; end; end

# Adapt SDP atom tests from Convex.jl:
#   https://github.com/JuliaOpt/Convex.jl/blob/master/test/test_sdp.jl
# facts("[sdp] Test problem #1") do
# for solver in sdp_solvers
# context("With solver $(typeof(solver))") do
#     m = Model(solver=solver)
#     @variable(m, X[1:2,1:2], SDP)
#     @objective(m, Max, X[1,1])
#     stat = solve(m)
#     @fact stat --> :Unbounded

#     setobjectivesense(m, :Min)
#     stat = solve(m)
#     @fact stat --> :Optimal
#     @fact getobjectivevalue(m) --> roughly(0, 1e-4)

#     @constraint(m, X[1,1] == 1)
#     stat = solve(m)
#     @fact stat --> :Optimal
#     @fact getobjectivevalue(m) --> roughly(1, 1e-5)
#     @fact getvalue(X[1,1]) --> roughly(1, 1e-5)

#     @objective(m, Min, sum(diag(X)))
#     stat = solve(m)
#     @fact stat --> :Optimal
#     @fact getobjectivevalue(m) --> roughly(1, 1e-5)
# end; end; end

facts("[sdp] Test problem #2") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, Y[1:3,1:3], SDP)
    @constraint(m, Y[2,1] <= 4)
    @constraint(m, Y[2,2] >= 3)
    @objective(m, Min, trace(Y))
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getobjectivevalue(m) --> roughly(3, 1e-5)
end; end; end

facts("[sdp] Test problem #3") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, x >= 0)
    @variable(m, Y[1:3,1:3], SDP)
    @constraint(m, Y[2,1] == 1)
    @objective(m, Min, Y[1,2])
    stat = solve(m)
    # @fact stat --> :Optimal
    @fact getobjectivevalue(m) --> roughly(1, 1e-4)
end; end; end

facts("[sdp] Test problem #4") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, x >= 0)
    @variable(m, Y[1:3,1:3], SDP)
    @constraint(m, x >= 1)
    @constraint(m, Y[2,1] == 1)
    @objective(m, Min, x + Y[1,1])
    stat = solve(m)
    # @fact stat --> :Optimal # TODO: remove this once SCS starts behaving
    @fact getobjectivevalue(m) --> roughly(1, 1e-3)
end; end; end

function nuclear_norm(model, A)
    m, n = size(A,1), size(A,2)
    @variable(model, U[1:m,1:m])
    @variable(model, V[1:n,1:n])
    @SDconstraint(model, 0 ⪯ [U A; A' V])
    return 0.5(trace(U) + trace(V'))
end

facts("[sdp] Test problem #5") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, Y[1:3,1:3], SDP)
    @constraint(m, Y[2,1] <= 4)
    @constraint(m, Y[2,2] >= 3)
    @constraint(m, Y[3,3] <= 2)
    @objective(m, Min, nuclear_norm(m, Y))
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getobjectivevalue(m) --> roughly(3, 1e-5)
end; end; end

function operator_norm(model, A)
    m, n = size(A,1), size(A,2)
    @variable(model, t >= 0)
    @SDconstraint(model, [t*eye(n) A; A' eye(n)*t] >= 0)
    return t
end

facts("[sdp] Test problem #6") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, Y[1:3,1:3])
    @constraint(m, Y[2,1] <= 4)
    @constraint(m, Y[2,2] >= 3)
    @constraint(m, sum(Y) >= 12)
    @objective(m, Min, operator_norm(m, Y))
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getobjectivevalue(m) --> roughly(4, 1e-5)
end; end; end

function lambda_max(model, A)
    m, n = size(A,1), size(A,2)
    @variable(model, t)
    @SDconstraint(model, speye(n)*t - A ⪰ 0)
    @SDconstraint(model, A >= 0)
    return t
end

facts("[sdp] Test problem #7") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, Y[1:3,1:3])
    @constraint(m, Y[1,1] >= 4)
    @objective(m, Min, lambda_max(m, Y))
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getobjectivevalue(m) --> roughly(4, 1e-5)
end; end; end

function lambda_min(model, A)
    m, n = size(A,1), size(A,2)
    @variable(model, t)
    @SDconstraint(model, A - eye(n)*t >= 0)
    @SDconstraint(model, A >= 0)
    return t
end

facts("[sdp] Test problem #8") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, Y[1:3,1:3], SDP)
    @constraint(m, trace(Y) <= 6)
    @objective(m, Max, lambda_min(m, Y))
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getobjectivevalue(m) --> roughly(2, 1e-5)
end; end; end

function matrix_frac(model, x, P)
    n = size(P,1)
    @variable(model, t)
    @variable(model, M[1:(n+1),1:(n+1)], SDP)
    @constraint(model, M[1:n,1:n] .== P)
    @constraint(model, M[1:n,n+1] .== x)
    @constraint(model, M[n+1,n+1] == t)
    return t
end

facts("[sdp] Test problem #9") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    n = 3
    x = [1,2,3]
    lb = 0.5eye(n)
    ub = 2eye(n)
    @variable(m, lb[i,j] <= P[i=1:n,j=1:n] <= ub[i,j])
    @objective(m, Min, matrix_frac(m, x, P))
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getobjectivevalue(m) --> roughly(7, 1e-5)
end; end; end

facts("[sdp] Correlation example") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)

    @variable(m, X[1:3,1:3], SDP)

    # Diagonal is 1s
    @constraint(m, X[1,1] == 1)
    @constraint(m, X[2,2] == 1)
    @constraint(m, X[3,3] == 1)

    # Bounds on the known correlations
    @constraint(m, X[1,2] >= -0.2)
    @constraint(m, X[1,2] <= -0.1)
    @constraint(m, X[2,3] >=  0.4)
    @constraint(m, X[2,3] <=  0.5)

    # Find upper bound
    @objective(m, Max, X[1,3])
    stat = solve(m)
    @fact stat --> :Optimal
    @fact (+0.8719 <= getvalue(X)[1,3] <= +0.8720) --> true

    # Find lower bound
    @objective(m, Min, X[1,3])
    stat = solve(m)
    @fact stat --> :Optimal
    @fact (-0.9779 >= getvalue(X)[1,3] >= -0.9799) --> true
end; end; end

facts("[sdp] Robust uncertainty example") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    include(joinpath("data","robust_uncertainty.jl"))
    R = 1
    d = 3
    𝛿 = 0.05
    ɛ = 0.05
    N = ceil((2+2log(2/𝛿))^2) + 1

    Γ1(𝛿,N) = (R/sqrt(N))*(2+sqrt(2*log(1/𝛿)))
    Γ2(𝛿,N) = (2R^2/sqrt(N))*(2+sqrt(2*log(2/𝛿)))

    for d in [3,5,8]; context("d = $d") do

        μhat = μhats[d]
        M = Ms[d]
        Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat')

        m = Model(solver=solver)

        @variable(m, Σ[1:d,1:d], SDP)
        @variable(m, u[1:d])
        @variable(m, μ[1:d])

        @variable(m, t1 >= 0)
        @variable(m, L1[1:d])
        @constraint(m, L1 .== (μ-μhat))
        @constraint(m, sum{L1[i]^2, i=1:d} <= t1^2)
        @constraint(m, t1 <= Γ1(𝛿/2,N))

        @variable(m, t2 >= 0)
        @variable(m, L2[1:d,1:d])
        @constraint(m, L2 .== (Σ-Σhat))
        @constraint(m, sum{L2[i,j]^2, i=1:d, j=1:d} <= t2^2)
        @constraint(m, t2 <= Γ2(𝛿/2,N))

        A = [(1-ɛ)/ɛ (u-μ)';
             (u-μ)     Σ   ]
        @SDconstraint(m, A >= 0)

        c = cs[d]
        @objective(m, Max, dot(c,u))

        stat = solve(m)

        object = getobjectivevalue(m)
        exact = dot(μhat,c) + Γ1(𝛿/2,N)*norm(c) + sqrt((1-ɛ)/ɛ)*sqrt(dot(c,(Σhat+Γ2(𝛿/2,N)*eye(d,d))*c))
        @fact stat --> :Optimal
        @fact abs(object - exact) --> roughly(0, 1e-5)
    end; end
end; end; end

facts("[sdp] Robust uncertainty example (with norms)") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    include(joinpath("data","robust_uncertainty.jl"))
    R = 1
    d = 3
    𝛿 = 0.05
    ɛ = 0.05
    N = ceil((2+2log(2/𝛿))^2) + 1

    Γ1(𝛿,N) = (R/sqrt(N))*(2+sqrt(2*log(1/𝛿)))
    Γ2(𝛿,N) = (2R^2/sqrt(N))*(2+sqrt(2*log(2/𝛿)))

    for d in [3,5,8]; context("d = $d") do

        μhat = μhats[d]
        M = Ms[d]
        Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat')

        m = Model(solver=solver)

        @variable(m, Σ[1:d,1:d], SDP)
        @variable(m, u[1:d])
        @variable(m, μ[1:d])

        @constraint(m, norm(μ-μhat) <= Γ1(𝛿/2,N))
        @constraint(m, vecnorm(Σ-Σhat) <= Γ2(𝛿/2,N))

        A = [(1-ɛ)/ɛ (u-μ)';
             (u-μ)     Σ   ]
        @SDconstraint(m, A >= 0)

        c = cs[d]
        @objective(m, Max, dot(c,u))

        stat = solve(m)

        object = getobjectivevalue(m)
        exact = dot(μhat,c) + Γ1(𝛿/2,N)*norm(c) + sqrt((1-ɛ)/ɛ)*sqrt(dot(c,(Σhat+Γ2(𝛿/2,N)*eye(d,d))*c))
        @fact stat --> :Optimal
        @fact abs(object - exact) --> roughly(0, 1e-5)
    end; end
end; end; end

facts("[sdp] Can't mix SDP and QP (#665)") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    model = Model(solver=solver)
    @variable(model, Q[1:2, 1:2], SDP)
    @constraint(model, Q[1,1] - 1 == Q[2,2])
    @objective(model, Min, Q[1,1] * Q[1,1])
    @fact_throws solve(model)
end; end; end

facts("[sdp] Just another SDP") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    model = Model(solver=solver)
    @variable(model, Q[1:2, 1:2], SDP)
    @constraint(model, Q[1,1] - 1 == Q[2,2])
    @variable(model, objective)
    T = [1 Q[1,1]; Q[1,1] objective]
    @SDconstraint(model, T ⪰ 0)
    @objective(model, Min, objective)
    @fact solve(model) --> :Optimal
    @fact getvalue(Q) --> roughly([1 0;0 0], 1e-3)
    @fact getobjectivevalue(model) --> roughly(1, TOL)
end; end; end


text:
function UpdateThermo!(sys::System,comm::MPI.Comm,root::Int64)
    UpdateTemperature!(sys,comm,root)
    UpdateMomentum!(sys,comm,root)
end

function kinetic_energy(vel::Array{Float64,1},masses::Float64)
    k=vel[1]*vel[1]*0.5*masses+vel[2]*vel[2]*0.5*masses+vel[3]*vel[3]*0.5*masses
    return k
end

function UpdateTemperature!(sys::System,comm::MPI.Comm,root::Int64)
    Ktot=0.0
    Ktot_local = sum(kinetic_energy.(sys.atoms.vels,sys.atoms.masses))
    MPI.Barrier(comm)
    Ktot=MPI.Allreduce(Ktot_local, +, comm)
    MPI.Barrier(comm)
    temp=Ktot*2/(sys.atoms.total_N[1]*sys.pbc.dim-sys.pbc.dim)
    sys.current_temp[1]=temp
end

function UpdateMomentum!(sys::System,comm::MPI.Comm,root::Int64)
    momentum_local = sum(sys.atoms.vels .* sys.atoms.masses)
    momentum_tot=MPI.Allreduce(momentum_local, +, comm)
    momentum_particle=momentum_tot/(sys.atoms.total_N[1]*sys.pbc.dim-sys.pbc.dim)
    sys.current_momentum[1]=momentum_particle[1]+momentum_particle[2]+momentum_particle[3]
end


text:
using Plots
using Random
Random.seed!(1)

include_ddp()

# Model
include_model("double_integrator")

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2]; u[1]]
end

function fd(model::DoubleIntegratorContinuous{Midpoint, FixedTime}, x, u, w, h, t)
	if t == 1
		# h = max(0.001, u[2])
		h = u[2]^2.0
    	return [x + h * f(model, x + 0.5 * h * f(model, x, view(u, 1:1), w), view(u, 1:1), w); h]
	else
		h = x[3]^2.0
		return [view(x, 1:2) + h * f(model, view(x, 1:2) + 0.5 * h * f(model, view(x, 1:2), u, w), u, w); h]
	end
end

# Time
T = 11
h = 0.1

model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 1)
n = [t == 1 ? model.n : model.n + 1 for t = 1:T]
m = [t == 1 ? model.m + 1 : model.m for t = 1:T]

# Initial conditions, controls, disturbances
x1 = [1.0; 0.0]
xT = [0.0; 0.0]
ū = [t == 1 ? [1.0 * randn(1); h] : 1.0 * randn(1) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# Objective
obj = StageCosts([NonlinearCost() for t = 1:T], T)

Q = [t < T ? Diagonal([1.0; 1.0]) : Diagonal([1.0; 1.0]) for t = 1:T]
q = [-2.0 * Q[t] * xT for t = 1:T]
R = [Diagonal(1.0 * ones(1)) for t = 1:T-1]
r = [[0.0] for t = 1:T-1]

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		t == 1 ? (h = u[2]) : (h = x[3])
        return h * (view(x, 1:2)' * Q[t] * view(x, 1:2) + q[t]' * view(x, 1:2) + view(u, 1:1)' * R[t] * view(u, 1:1) + r[t]' * view(u, 1:1)) + 1.0e-5 * u' * u
    elseif t == T
        return view(x, 1:2)' * Q[t] * view(x, 1:2) + q[t]' * view(x, 1:2)
    else
        return 0.0
    end
end

# Constraints
p = [t == 1 ? 4 : (t == T ? 2 : 2) for t = 1:T]
info_1 = Dict(:ul => [-5.0; 0.01], :uu => [5.0; 1.0], :inequality => (1:4))
info_t = Dict(:ul => [-5.0], :uu => [5.0], :inequality => (1:2))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? (t == 1 ? info_1 : info_t) : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	if t == 1
		c[1:2] .= view(u, 1:2) - cons.con[t].info[:uu]
		c[3:4] = cons.con[t].info[:ul] - view(u, 1:2)
	end
	if t > 1 && t < T
		c[1:1] .= view(u, 1:1) - cons.con[t].info[:uu]
		c[2:2] = cons.con[t].info[:ul] - view(u, 1:1)
	end
	if t == T
		c[1:2] .= view(x, 1:2) - cons.con[t].info[:xT]
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T, n = n, m = m)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 6,
	ρ_init = 1.0, ρ_scale = 10.0)

# # Solve
# @time ddp_solve!(prob,
#     max_iter = 100, verbose = true)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

@show u[1][2] * (T - 1)
@show ū[1][2] * (T - 1)

# Visualize
using Plots
plot(hcat([xt[1:2] for xt in x]...)[:, 1:T]', color = :black, label = "")
plot(hcat([ut[1:1] for ut in u]...)[:, 1:T-1]', color = :black, label = "", linetype = :steppost)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
    project(M::AbstractManifold, p)

Project point `p` from the ambient space of the [`AbstractManifold`](@ref) `M` to `M`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given. Additionally, the projection includes changing data representation, if applicable,
i.e. if the points on `M` are not represented in the same array data, the data is changed
accordingly.

See also: [`EmbeddedManifold`](@ref), [`embed`](@ref embed(M::AbstractManifold, p))
"""
function project(M::AbstractManifold, p)
    q = allocate_result(M, project, p)
    project!(M, q, p)
    return q
end

"""
    project!(M::AbstractManifold, q, p)

Project point `p` from the ambient space onto the [`AbstractManifold`](@ref) `M`. The result is
storedin `q`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given. Additionally, the projection includes changing data representation, if applicable,
i.e. if the points on `M` are not represented in the same array data, the data is changed
accordingly.

See also: [`EmbeddedManifold`](@ref), [`embed!`](@ref embed!(M::AbstractManifold, q, p))
"""
function project!(M::AbstractManifold, q, p)
    return error(manifold_function_not_implemented_message(M, project!, q, p))
end

"""
    project(M::AbstractManifold, p, X)

Project ambient space representation of a vector `X` to a tangent vector at point `p` on
the [`AbstractManifold`](@ref) `M`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given.
Additionally, `project` includes changing data representation, if applicable, i.e.
if the tangents on `M` are not represented in the same way as points on the embedding,
the representation is changed accordingly. This is the case for example for Lie groups,
when tangent vectors are represented in the Lie algebra. after projection the change to the
Lie algebra is perfomed, too.

See also: [`EmbeddedManifold`](@ref), [`embed`](@ref embed(M::AbstractManifold, p, X))
"""
function project(M::AbstractManifold, p, X)
    # Note that the order is switched,
    # since the allocation by default takes the type of the first.
    Y = allocate_result(M, project, X, p)
    project!(M, Y, p, X)
    return Y
end

"""
    project!(M::AbstractManifold, Y, p, X)

Project ambient space representation of a vector `X` to a tangent vector at point `p` on
the [`AbstractManifold`](@ref) `M`. The result is saved in vector `Y`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given.
Additionally, `project!` includes changing data representation, if applicable, i.e.
if the tangents on `M` are not represented in the same way as points on the embedding,
the representation is changed accordingly. This is the case for example for Lie groups,
when tangent vectors are represented in the Lie algebra. after projection the change to the
Lie algebra is perfomed, too.

See also: [`EmbeddedManifold`](@ref), [`embed!`](@ref embed!(M::AbstractManifold, Y, p, X))
"""
function project!(M::AbstractManifold, Y, p, X)
    return error(manifold_function_not_implemented_message(M, project!, Y, p, X))
end


text:
using Makie, VideoIO, Dates, Observables

include("testvideo.jl")

# create the test video if it doesn't exist
if !isfile(testvideo)
    createtestvideo(testvideo)
end

# open the video
f = VideoIO.openvideo(testvideo)

# read one image to get its type <-- this might be unecessary since we might know that in advance for all videos
img = Observable(read(f))
timestamp = Observable(Time(0))

# reset the pointer to the begining of the file <-- there might be a more elegant way to achieve this
f = VideoIO.openvideo(testvideo)

# create a buffer of 2 seconds
nframesbuffer = round(Int, 2f.framerate)
buff = Channel{Pair{Time, typeof(img[])}}(nframesbuffer)
# fill up the buffer with the next 15 seconds worth of frames
task = @async begin
    while !eof(f)
        t = Time(0) + Millisecond(round(Int, 1000gettime(f)))
        put!(buff, Pair(t, read(f)))
        yield() # needed
    end
end

# all the stuff needed for display
pixelaspectratio = VideoIO.aspect_ratio(f)
h = f.height
w = round(typeof(h), f.width * pixelaspectratio)
AbstractPlotting.set_theme!(backgroundcolor = :black)
scene = Scene(resolution = (w, h))
Makie.image!(scene, 1:h, 1:w, img, show_axis = false, scale_plot = false)
Makie.rotate!(scene, -0.5π)
update_limits!(scene)
update_cam!(scene)
update!(scene)

# everything starts at this observable, the current frame
old_slider_pos = Observable(1)

# get an estimate for the number of frames, this could be very innacurate for VFR videos
nframes = round(Int, VideoIO.get_duration(f.avin.io)*f.framerate)

# the slider, updating currentframe and getting its string from timestamps
slider_h = slider(1:nframes, start = 1, valueprinter = i -> "no, this will be wrong", textcolor = :white)

lift(slider_h[end][:value]) do frame
    while old_slider_pos[] ≤ frame
        timestamp[], img[] = take!(buff)
        old_slider_pos[] += 1
    end
    0
    # currentframe[] = frame
end

# this is for a separate display of the time stamp, due to the fact that we can not calculate with the time stamp is from the frame number (which is what the slider has)
timestamp_h = text(lift(string, timestamp), color = :white, raw = true, camera = campixel!)

# I have to have a forward button becuase apparently the stuff in the lift gets evaluated, so to have the movie start in frame #1 I need to go forward once and the backwards...
fwdbutton = button(">", textcolor = :white)
lift(fwdbutton[end][:clicks]) do _
    # if currentframe[] < nframes # check we are not in the last frame
    timestamp[], img[] = take!(buff)
    slider_h[end][:value][] += 1
        # so now the movie is in frame #2
    # end
    0
end
#=bckbutton = button!(rsc(), "<", textcolor = :white)
lift(bckbutton[end][:clicks]) do _
    if currentframe[] > 1 # check we're not in frame #1
        currentframe[] -= 1
        slider_h[end][:value][] -= 1
        # and now we're back in frame #1
    end
    0
end=#

# done!
hbox(vbox(slider_h, fwdbutton, timestamp_h), scene)

# played = Observables.async_latest(play, play_bubtton[end][:clicks])

# hbox(vbox(play_button, next_button, slider_h), scene)


text:
function solve_perspective_relaxation(U, mu, lambda, k_sparse, k_rank;
    zero_indices=[], one_indices=[], solver_output=0)
    """
    This function exactly solves the perspective reformulation relaxation of the
    problem of interest.

    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
                   (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).
    :param zero_indices: List of 2-tuples of Int64 where each entry consists of
                         an index (i, j) such that Y_ij is constrained to take
                         value 0.
    :param one_indices: List of 2-tuples of Int64 where each entry consists of
                        an index (i, j) such that S_ij is constrained to take
                        value 1 where S is the binary matrix denoting the
                        sparsity pattern of Y.
    :param solver_output: Solver_output param to be passed to Mosek (Int64).

    :return: This function returns two values. The first value is a tuple of 5
             n-by-n arrays that correspond to the solution of the perspective
             relaxation (the first element in the tuple is the matrix X, the
             second element is the matrix Y, the third element is the matrix Z,
             the fourth element is the matrix P and the fifth element is the
             matrix Theta). The second value is the optimal objective value of
             the relaxation (Float64).
    """

    n = size(U)[1]

    num_zeros = size(zero_indices)[1]
    num_ones = size(one_indices)[1]

    @assert num_ones <= k_sparse
    @assert num_zeros <= n * n - k_sparse

    # Build perspective formulation relaxation using JuMP
    m = Model(Mosek.Optimizer)
    set_optimizer_attribute(m, "MSK_IPAR_LOG", solver_output)

    @variable(m, X[i=1:n, j=1:n])
    @variable(m, Y[i=1:n, j=1:n])
    @variable(m, Z[i=1:n, j=1:n] >= 0)
    @variable(m, P[i=1:n, j=1:n])
    @variable(m, P_hat[i=1:n, j=1:n])
    @variable(m, alpha[i=1:n, j=1:n])
    @variable(m, theta[i=1:n, j=1:n] >= 0)

    @constraint(m, [i=1:n, j=1:n], Z[i, j] <= 1)
    @SDconstraint(m, [i=1:n, j=1:n],
        [theta[i, j] Y[i, j]; Y[i, j] Z[i, j]] >= 0)
    for (i, j) in zero_indices
        @constraint(m, Z[i, j] == 0)
        @constraint(m, Y[i, j] == 0)
        @constraint(m, theta[i, j] == 0)
    end
    for (i, j) in one_indices
        @constraint(m, Z[i, j] == 1)
        @constraint(m, Y[i, j] == (U[i, j] - X[i, j]) / (1 + mu))
    end
    @constraint(m, sum(Z[i, j] for i=1:n, j=1:n) <= k_sparse)

    @constraint(m, sum(P[i, i] for i = 1:n) <= k_rank)
    @constraint(m, 1*Matrix(I, n, n) - P in PSDCone())
    @constraint(m, [alpha X; X' P] in PSDCone())

    @objective(m, Min, sum((U[i, j]-X[i, j]-Y[i, j])^2 for i=1:n for j=1:n) +
        lambda * sum(alpha[i, i] for i=1:n) +
        mu * sum(theta[i, j] for i=1:n for j=1:n))

    # Solve perspective relaxation
    optimize!(m)

    X_sol = value.(X)
    Y_sol = value.(Y)

    P_sol = value.(P)
    Z_sol = value.(Z)
    alpha_sol = value.(alpha)
    theta_sol = value.(theta)

    return (X_sol, Y_sol, Z_sol, P_sol, alpha_sol), objective_value(m)

end;


text:
using SafeTestsets
import Base.Threads.@spawn

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
const is_TRAVIS = haskey(ENV,"TRAVIS")

# Start Test Script

@time begin

    # if GROUP == "All" || GROUP == "Components"
    #    #@time @safetestset "Test for regression against original code" begin include("regression_test.jl") end
    #    @time @safetestset "MOLFiniteDifference Utils" begin include("utils_test.jl") end
    #    @time @safetestset "Discretization of space and grid types" begin include("components/DiscreteSpace.jl") end
    #    @time @safetestset "Variable PDE mapping and interior construction" begin include("components/interiormap_test.jl") end
    #    @time @safetestset "Fornberg" begin include("components/MOLfornberg_weights.jl") end
    #    #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end
    # end

    # if GROUP == "All" || GROUP == "Interface"
    #     @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest.jl") end
    # end

    # if GROUP == "All" || GROUP == "Diffusion"
    #     @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion" begin include("pde_systems/MOL_1D_Linear_Diffusion.jl") end
    # end

    # if GROUP == "All" || GROUP == "Nonlinear_Diffusion"
    #     @time @safetestset "MOLFiniteDifference Interface: 1D Non-Linear Diffusion" begin include("pde_systems/MOL_1D_NonLinear_Diffusion.jl") end
    # end

    # if GROUP == "All" || GROUP == "2D_Diffusion"
    #     @time @safetestset "MOLFiniteDifference Interface: 2D Diffusion" begin include("pde_systems/MOL_2D_Diffusion.jl") end
    # end

    # if GROUP == "All" || GROUP == "Convection"
    #     @time @safetestset "MOLFiniteDifference Interface: Linear Convection" begin include("pde_systems/MOL_1D_Linear_Convection.jl") end
    # end

    # if GROUP == "All" || GROUP == "Higher_Order"
    #     @time @safetestset "MOLFiniteDifference Interface: 1D HigherOrder" begin include("pde_systems/MOL_1D_HigherOrder.jl") end
    # end

    # if GROUP == "All" || GROUP == "Stationary"
    #     @time @safetestset "MOLFiniteDifference Interface: Stationary Nonlinear Problems" begin include("pde_systems/MOL_NonlinearProblem.jl") end
    # end

    # if GROUP == "All" || GROUP == "DAE"
    #     @time @safetestset "MOLFiniteDifference Interface: 1D Partial DAE" begin include("pde_systems/MOL_1D_PDAE.jl") end
    # end

    # if GROUP == "All" || GROUP == "Burgers"
    #     @time @safetestset "MOLFiniteDifference Interface: 2D Burger's Equation" begin include("pde_systems/burgers_eq.jl") end
    # end

    if GROUP == "All" || GROUP == "Brusselator"
        @time @safetestset "MOLFiniteDifference Interface: 2D Brusselator Equation" begin include("pde_systems/brusselator_eq.jl") end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

code = create_basis(
    :Tri3,
    "3 node linear triangle element",
    (
     (0.0, 0.0), # N1
     (1.0, 0.0), # N2
     (0.0, 1.0), # N3
    ),
    "1 + u + v",
   )
eval(code)

code = create_basis(
    :Tri6,
    "6 node quadratic triangle element",
    (
     (0.0, 0.0), # N1
     (1.0, 0.0), # N2
     (0.0, 1.0), # N3
     (0.5, 0.0), # N4
     (0.5, 0.5), # N5
     (0.0, 0.5), # N6
    ),
    "1 + u + v + u^2 + u*v + v^2",
   )
eval(code)

code = create_basis(
    :Tri7,
    "7 node quadratic triangle element (has middle node)",
    (
     (0.0, 0.0), # N1
     (1.0, 0.0), # N2
     (0.0, 1.0), # N3
     (0.5, 0.0), # N4
     (0.5, 0.5), # N5
     (0.0, 0.5), # N6
     (1/3, 1/3), # N7
    ),
    "1 + u + v + u^2 + u*v + v^2 + u^2*v^2",
   )
eval(code)


text:
using Whitespaces
using Test

@testset "Whitespaces.jl" begin
    # Write your own tests here.
    hi = "   \t  \t   \n\t\n     \t\t \t  \t\n\t\n  \n\n\n"
    hi_ws = Whitespace(hi)
    compile!(hi_ws)
    res = sprint(io->execute!(io, hi_ws))
    @test res == "Hi"

    seven = Whitespace("   \t\t\t\n   \t\n\t   \t\n \t\n\n\n")
    compile!(seven)
    res = sprint(io->execute!(io, seven))
    @test res == "8"
end


text:
# Application Product Test
module APT
include("../src/sdistance.jl")
include("../src/environments.jl")
import .Sdistance: computing_bench, signedDistance2D
import Plots
using Plots
function plots_contours(Ns::Array, _csv_datafiles, circle_n::Union{String,Nothing} = nothing)
    plots = []
    for (i, N) in enumerate(Ns)
        _x = [i for i = -L:2*L/N:L] # len:N+1
        _y = [i for i = -L:2*L/N:L] # len:N+1
        _phi = signedDistance2D(_csv_datafiles, N, circle_n)
        push!(plots, plot(_x, _y, _phi, title = "N=$N", st = :contour))
    end
    println(plots)
    plot(plots..., size = (2300, 2500), layout = (3, 2))
    if circle_n == "multi"
        savefig("test/image/NsPlots_floodfill.png")
    else
        savefig("test/image/NsPlots_isinside.png")
    end
end
function plots_wireframe(Ns::Array, _csv_datafiles, circle_n::Union{String,Nothing} = nothing)
    plots = []
    for (i, N) in enumerate(Ns)
        _x = [i for i = -L:2*L/N:L] # len:N+1
        _y = [i for i = -L:2*L/N:L] # len:N+1
        _phi = signedDistance2D(_csv_datafiles, N, circle_n)
        push!(plots, plot(_x, _y, _phi, title = "N=$N", st = :wireframe))
    end
    println(plots)
    plot(plots..., size = (2300, 2500), layout = (3, 2))
    if circle_n == "multi"
        savefig("test/image/NsPlots_floodfill_wireframe.png")
    else
        savefig("test/image/NsPlots_isinside_wireframe.png")
    end
end
export plots_contours, plots_wireframe
end

max_stars_repo_name: jabelic/senior-thesis-julia

max_issues_repo_name: jabelic/senior-thesis-julia

max_forks_repo_name: jabelic/senior-thesis-julia


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct RoleType <: AbstractRoleType
    label::Label
    is_root::Bool
end

# Porting note: Java client calls into RequstBuilder but it really
# has nothing to do with requests... I think it's probably better
# migrating the function here.
function proto(t::AbstractRoleType)
    @assert t.label.scope !== nothing
    Proto._Type(
        label = t.label.name,
        scope = t.label.scope,
        encoding = encoding(t)
    )
end

function get_players(transaction::AbstractCoreTransaction, label::Label)
    play_req = RoleTypeRequestBuilder.get_players_req(label)
    player_res = stream(transaction, play_req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.role_type_get_players_res_part.thing_types for r in player_res)))
end


text:
# from: https://github.com/BBN-Q/Instruments.jl

# Helper functions to find instruments
# resmgr = viOpenDefaultRM()
# viFindRsrc(resmgr, expr)

function check_status(status)
	status_signed = reinterpret(Int32, UInt32(status))
	if status_signed != VI_SUCCESS
		errMsg = codes[status_signed]
		error("VISA C call failed with status $(errMsg[1]): $(errMsg[2])")
	end
	status_signed
end

#- Resource Manager Functions and Operations -------------------------------#
function viOpenDefaultRM()
	rm = ViPSession(0)
	check_status(ccall((:viOpenDefaultRM, libvisa), ViStatus, (ViPSession,), rm))
	rm.x
end

function viFindRsrc(sesn::ViSession, expr::AbstractString)
	returnCount = ViPUInt32(0)
	findList = ViPFindList(0)
	desc = zeros(ViByte, VI_FIND_BUFLEN)
	descp = pointer(desc)
	check_status(ccall((:viFindRsrc, libvisa), ViStatus,
						(ViSession, ViString, ViPFindList, ViPUInt32, ViPByte),
						sesn, expr, findList, returnCount, descp))

	#Create the array of instrument strings and push them on
	instrStrs = Vector{String}()
	if returnCount.x > 0
		push!(instrStrs, unsafe_string(descp))
	end
	for i=1:returnCount.x-1
		check_status(ccall((:viFindNext, libvisa), ViStatus,
						(ViFindList, ViPByte), findList.x, descp))
		push!(instrStrs, unsafe_string(descp))

	end

	instrStrs
end

function viClose(viObj::ViObject)
	check_status(ccall((:viClose, libvisa), ViStatus, (ViObject,), viObj))
end

# #- Resource Template Operations --------------------------------------------*/

function viSetAttribute(viObj::ViObject, attrName::ViAttr, attrValue::ViAttrState)
	check_status(ccall((:viSetAttribute, libvisa), ViStatus,
						(ViObject, ViAttr, ViAttrState),
						viObj, attrName, attrValue))
end

function viGetAttribute(viObj::ViObject, attrName::ViAttr)
	value = ViAttrState[0]
	check_status( ccall((:viGetAttribute, libvisa), ViStatus,
						(ViObject, ViAttr, Ptr{Cvoid}),
						viObj, attrName, value))
	value[]
end

function viEnableEvent(instrHandle::ViSession, eventType::Integer,
					   mechanism::Integer)
	check_status(ccall((:viEnableEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt16, ViEventFilter),
						 instrHandle, eventType, mechanism, 0))
end

function viDisableEvent(instrHandle::ViSession, eventType::Integer,
					   mechanism::Integer)
	check_status(ccall((:viEnableEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt16),
						 instrHandle, eventType, mechanism))
end

function viDiscardEvents(instrHandle::ViSession, eventType::ViEventType,
					   mechanism::UInt16)
	check_status(ccall((:viEnableEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt16),
						 instrHandle, eventType, mechanism))
end

function viWaitOnEvent(instrHandle::ViSession, eventType::ViEventType, timeout::UInt32 = VI_TMO_INFINITE)
	outType = Array(ViEventType)
	outEvent = Array(ViEvent)
	check_status(ccall((:viWaitOnEvent,libvisa), ViStatus,
						(ViSession, ViEventType, UInt32, Ptr{ViEventType}, Ptr{ViEvent}),
						 instrHandle, eventType, timeout, outType, outEvent))
	(outType[], outEvent[])
end

#- Basic I/O Operations ----------------------------------------------------#

function viWrite(instrHandle::ViSession, data::Union{String, Vector{UInt8}})
	bytesWritten = ViUInt32[0]
	check_status(ccall((:viWrite, libvisa), ViStatus,
						(ViSession, ViBuf, ViUInt32, ViPUInt32),
						instrHandle, pointer(data), length(data), bytesWritten))
	bytesWritten[1]
end

function viRead!(instrHandle::ViSession, buffer::Array{UInt8})
	bytesRead = ViUInt32[0]
	status = check_status(ccall((:viRead, libvisa), ViStatus,
						(ViSession, ViBuf, ViUInt32, ViPUInt32),
						instrHandle, buffer, sizeof(buffer), bytesRead))
	return (status != VI_SUCCESS_MAX_CNT, bytesRead[])
end

function viRead(instrHandle::ViSession; bufSize::UInt32=0x00000400)
	buf = zeros(UInt8, bufSize)
	(done, bytesRead) = viRead!(instrHandle, buf)
	unsafe_string(pointer(buf))
end

function readavailable(instrHandle::ViSession)
	ret = IOBuffer()
	buf = Array(UInt8, 0x400)
	while true
		(done, bytesRead) = viRead!(instrHandle, buf)
		write(ret,buf[1:bytesRead])
		if done
			break
		end
	end
	take!(ret)
end

write(obj::T, msg::AbstractString) where {T<:Instrument} = viWrite(obj.handle, msg)

read(obj::T) where {T<:Instrument} = rstrip(viRead(obj.handle; bufSize=obj.bufSize), ['\r', '\n'])

readavailable(obj::T) where {T<:Instrument} = readavailable(obj.handle)

# connect,disconnect, query
function connect!(sesn, obj, mode=VI_NO_LOCK, timeout=VI_TMO_IMMEDIATE)
	if !obj.initialized
		#Pointer for the instrument handle
		vi = ViPSession(0)
		status = viOpen(sesn, obj.address, mode, timeout, vi)
		println("status: ", status)
		check_status(status)
		obj.handle = vi.x
		obj.initialized = true
	end
end

function disconnect!(obj::T) where {T<:Instrument}
	if obj.initialized
		viClose(obj.handle)
		obj.initialized = false
	end
end

function query(obj::T, msg::AbstractString; delay::Real=0) where {T<:Instrument}
	write(obj, msg)
	sleep(delay)
	read(obj)
end


text:
baremodule Orihime

macro plot end

function cleanup end
function readablestring end

module Internal

import ..Orihime: @plot
using ..Orihime: Orihime

import VegaLite
using DataFrames: DataFrames, DataFrame

include("utils.jl")
include("cleanup.jl")
include("plot.jl")

end  # module Internal

end  # baremodule Orihime


text:
# %% CCM illustration
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

ds = Systems.lorenz()
tr = trajectory(ds, 100; Ttr = 100)
x, y, z = columns(tr)

τx = estimate_delay(x, "mi_min")
τy = estimate_delay(y, "mi_min")
τz = estimate_delay(z, "mi_min")

X = embed(x, 3, τx)
Y = embed(y, 3, τy)
Z = embed(z, 3, τz)

using3D()
fig = figure()
axts = fig.add_subplot(1, 3, 1)
axx = fig.add_subplot(1,3,2; projection="3d")
axy = fig.add_subplot(1,3,3; projection="3d")

axts.plot(x .- 10)
axts.plot(y .+ 20; color = "C2")
axts.set_xlim(0, 1000)
axts.set_xticklabels([])
axts.set_yticklabels([])

axx.plot3D(columns(X)...; lw = 1)
axy.plot3D(columns(Y)...; lw = 1, color = "C2")
for s in (:x, :y, :z)
    f = Symbol(:set_, s, :ticklabels)
    @eval axx.$(f)([])
    @eval axy.$(f)([])
    g = Symbol(:set_, s, :lim)
    @eval axx.$(g)(-15, 15)
    @eval axy.$(g)(-20, 20)
end

# Axis pretty-fication
axx.dist = 8
axy.dist = 8
axx.elev = 20
axy.elev = 20
axx.text3D(-10, -10, 15, "\$M_x\$", size = 40)
axy.text3D(-10, -10, 20, "\$M_y\$", size = 40)
axts.set_ylabel("timeseries")

axx.text(15, 12, 21, "b"; bbox = bbox, zorder = 99, va = "top")
axy.text(17, 17, 27.5, "c"; bbox = bbox, zorder = 99, va = "top")
# axy.text(45, 40, 54, "c"; bbox = bbox, zorder = 99, va = "top")

add_identifiers!(fig)
fig.tight_layout(pad=0.35)
fig.subplots_adjust(wspace = 0.2)
wsave(plotsdir("7", "ccm"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Brownian
using Base.Test

@test convert(FGN, FBM(0:0.5:10, 0.4)) == FGN(0.757858283255199,0.4)

p = FBM(0:1/2^10:1, 0.4)

rand(p)
rand(p, fbm=false)
rand([p, p])
rand(p, rtype=:chol)
rand(p, fbm=false, rtype=:chol)
rand([p, p], rtype=:chol)


text:
## References
#==========================================================================================#

abstract type AbstractReferences end

"""
    References(V, ρ, μ, S, b, c, r)
    References(; speed, density, viscosity,
                 sound_speed, area, span,
                 chord, location)
Define reference values with speed ``V``, density ``ρ``, _dynamic_ viscosity ``μ``, area ``S``, span ``b``, chord ``c``, location ``r`` for a vortex lattice analysis.

# Arguments
- `speed       :: Real         = 1.`: Speed (m/s)
- `density     :: Real         = 1.225`: Density (m)
- `viscosity   :: Real         = 1.5e-5`: Dynamic viscosity (m)
- `sound_speed :: Real         = 330.`: Speed of sound (m/s)
- `span        :: Real         = 1.`: Span length (m)
- `area        :: Real         = 1.`: Area (m²)
- `chord       :: Real         = 1.`: Chord length (m)
- `location    :: Vector{Real} = [0,0,0]`: Position
"""
struct References{T} <: AbstractReferences
    speed       :: T
    density     :: T
    viscosity   :: T
    sound_speed :: T
    area        :: T
    span        :: T
    chord       :: T
    location    :: SVector{3,T}
end

function References(V, ρ, μ, a, S, b, c, r)
    T = promote_type(eltype(V), eltype(ρ), eltype(μ), eltype(a), eltype(S), eltype(b), eltype(c), eltype(r))
    References{T}(V, ρ, μ, a, S, b, c, r)
end

# References(; V, rho, mu, a, S, b, c, r) = References(V, rho, mu, a, S, b, c, r)

References(; speed = 1., density = 1.225, viscosity = 1.5e-5, sound_speed = 330., area = 1., span = 1., chord = 1., location = zeros(3)) = References(speed, density, viscosity, sound_speed, area, span, chord, location)

Base.broadcastable(refs :: References) = Ref(refs)

function Base.show(io :: IO, refs :: References)
    println(io, "References: ")
    for fname in fieldnames(typeof(refs))
        println(io, "    ", fname, " = ", getfield(refs, fname))
    end
end

kinematic_viscosity(refs :: References) = refs.viscosity / refs.density
mach_number(refs :: References) = refs.speed / refs.sound_speed

force_coefficient(force, refs :: References) = force_coefficient(force, dynamic_pressure(refs.density, refs.speed), refs.area)
moment_coefficient(moment, refs :: References) = moment_coefficient(moment, dynamic_pressure(refs.density, refs.speed), refs.area, refs.span, refs.chord)

rate_coefficient(fs :: Freestream, refs :: References) = rate_coefficient(fs.omega, refs.speed, refs.span, refs.chord)

## System
#==========================================================================================#

abstract type AbstractVortexLatticeSystem end

"""
    VortexLatticeSystem

A system consisting of the relevant variables for a vortex lattice analysis for post-processing.

# Arguments
The accessible fields are:
- `vortices`:  The array of vortices (currently `Horseshoe`).
- `circulations`: The circulation strengths of the vortices obtained by solving the linear system.
- `influence_matrix`: The influence matrix of the linear system.
- `boundary_vector`: The boundary condition corresponding to the right-hand-side of the linear system.
- `freestream :: Freestream`: The freestream conditions.
- `reference :: References`: The reference values.
"""
struct VortexLatticeSystem{M,N,R,S,P <: AbstractFreestream, Q <: AbstractReferences} <: AbstractVortexLatticeSystem
    vortices          :: M
    circulations      :: N
    influence_matrix  :: R
    boundary_vector   :: S
    freestream        :: P
    reference         :: Q
end

# Made out of annoyance and boredom
function Base.show(io :: IO, sys :: VortexLatticeSystem)
    println(io, "VortexLatticeSystem -")
    println(io, length(sys.vortices), " ", eltype(sys.vortices), " Elements\n")
    show(io, sys.freestream)
    println(io, "")
    show(io, sys.reference)
end

# Miscellaneous
rate_coefficient(system :: VortexLatticeSystem) = rate_coefficient(system.freestream, system.reference)

## THINK ABOUT USING ONLY WIND AXES FOR PG-TRANSFORMATION AND MAPPING BACK

# Velocities
"""
    surface_velocities(system :: VortexLatticeSystem;
                       axes   :: AbstractAxisSystem = Geometry())

Compute the induced velocities for all components of the `VortexLatticeSystem` in the reference axis system.
"""
surface_velocities(system :: VortexLatticeSystem; axes = Geometry()) = surface_velocities(system, axes)

surface_velocities(system :: VortexLatticeSystem, ::Geometry) = surface_velocities(system.vortices, system.vortices, system.circulations, system.reference.speed * body_frame_velocity(system.freestream), system.freestream.omega)

surface_velocities(system :: VortexLatticeSystem, ::Body) = geometry_to_body_axes.(surface_velocities(system, Geometry()), system.freestream.alpha, system.freestream.beta)

surface_velocities(system :: VortexLatticeSystem, ::Stability) = geometry_to_stability_axes.(surface_velocities(system, Geometry()), system.freestream.alpha)

surface_velocities(system :: VortexLatticeSystem, ::Wind) = geometry_to_wind_axes.(surface_velocities(system, Geometry()), system.freestream.alpha, system.freestream.beta)

## Forces
"""
    surface_forces(system :: VortexLatticeSystem;
                   axes   :: AbstractAxisSystem = Geometry())

Compute the forces for all components of the `VortexLatticeSystem` in the reference axis system.
"""
surface_forces(system; axes :: AbstractAxisSystem = Geometry()) = surface_forces(system, axes)

surface_forces(system :: VortexLatticeSystem, ::Geometry) = surface_forces(system.vortices, system.circulations, system.reference.speed * body_frame_velocity(system.freestream), system.freestream.omega, system.reference.density)

surface_forces(system :: VortexLatticeSystem, ::Body) = geometry_to_body_axes.(surface_forces(system, Geometry()), system.freestream.alpha, system.freestream.beta)

surface_forces(system :: VortexLatticeSystem, ::Stability) = geometry_to_stability_axes.(surface_forces(system, Geometry()), system.freestream.alpha)

surface_forces(system :: VortexLatticeSystem, ::Wind) = geometry_to_wind_axes.(surface_forces(system, Geometry()), system.freestream.alpha, system.freestream.beta)

## Moments
"""
    surface_moments(system :: VortexLatticeSystem;
                    axes   :: AbstractAxisSystem = Geometry())

Compute the moments for all components of the `VortexLatticeSystem` in the reference axis system.
"""
surface_moments(system; axes :: AbstractAxisSystem = Geometry()) = surface_moments(system, axes)

surface_moments(system :: VortexLatticeSystem, ::Geometry) = surface_moments(system.vortices, surface_forces(system, Geometry()), system.reference.location)

surface_moments(system :: VortexLatticeSystem, ::Body) = surface_moments(system.vortices, surface_forces(system, Body()), system.reference.location)

surface_moments(system :: VortexLatticeSystem, ::Stability) = surface_moments(system.vortices, flip_xz.(surface_forces(system, Stability())), system.reference.location)

surface_moments(system :: VortexLatticeSystem, ::Wind) =  surface_moments(system.vortices, flip_xz.(surface_forces(system, Wind())), system.reference.location)

## Dynamics
"""
    surface_dynamics(system :: VortexLatticeSystem;
                     axes   :: AbstractAxisSystem = Geometry())

Compute the forces and moments for all components of the `VortexLatticeSystem` in the reference axis system.
"""
function surface_dynamics(system :: VortexLatticeSystem)
    # Compute nearfield forces and moments
    surf_forces  = surface_forces(system)
    surf_moments = surface_moments(system.vortices, surf_forces, system.reference.location)

    surf_forces, surf_moments
end

surface_dynamics(system :: VortexLatticeSystem, ::Geometry) = surface_dynamics(system)

function surface_dynamics(system :: VortexLatticeSystem, ::Body)
    # Compute surface forces and moments in geometry axes
    surface_forces, surface_moments = surface_dynamics(system)

    # Transform to body axes
    stability_forces  = @. geometry_to_body_axes(surface_forces)
    stability_moments = @. geometry_to_body_axes(surface_moments)

    stability_forces, stability_moments
end

function surface_dynamics(system :: VortexLatticeSystem, ::Stability)
    # Get angle of attack
    α = system.freestream.alpha

    # Compute surface forces and moments
    surface_forces, surface_moments = surface_dynamics(system)

    # Transform to stability axes
    stability_forces  = @. geometry_to_stability_axes(surface_forces, α)
    stability_moments = @. geometry_to_stability_axes(flip_xz(surface_moments), α)

    stability_forces, stability_moments
end

function surface_dynamics(system :: VortexLatticeSystem, ::Wind)
    # Get angles of attack and sideslip
    α = system.freestream.alpha
    β = system.freestream.beta

    # Compute nearfield forces and moments
    surface_forces, surface_moments = surface_dynamics(system)

    # Transform to wind axes
    wind_forces  = @. geometry_to_wind_axes(surface_forces, α, β)
    wind_moments = @. geometry_to_wind_axes(flip_xz(surface_moments), α, β)

    wind_forces, wind_moments
end

surface_dynamics(system; axes :: AbstractAxisSystem = Wind()) = surface_dynamics(system, axes)

"""
    surface_coefficients(system :: VortexLatticeSystem;
                         axes   :: AbstractAxisSystem = Wind())

Compute the force and moment coefficients on the surface given the `VortexLatticeSystem` in the reference axis system.
"""
function surface_coefficients(system :: VortexLatticeSystem; axes :: AbstractAxisSystem = Wind())
    # Compute surface forces in whichever axes
    forces, moments = surface_dynamics(system, axes)

    # Compute coefficients
    CFs = @. force_coefficient(forces, system.reference)
    CMs = @. moment_coefficient(moments, system.reference)

    CFs, CMs
end

"""
    nearfield_coefficients(system :: VortexLatticeSystem)

Compute the force and moment coefficients in **wind axes** for all components of the `VortexLatticeSystem`.
"""
function nearfield_coefficients(system :: VortexLatticeSystem)
    CFs, CMs = surface_coefficients(system; axes = Wind())
    NamedTuple(key => [sum(CFs[key]); sum(CMs[key])] for key in keys(CFs))
end

"""
    nearfield(system :: VortexLatticeSystem)

Compute the **total** force and moment coefficients in **wind axes** for all components of the `VortexLatticeSystem`.
"""
function nearfield(system :: VortexLatticeSystem)
    CDi, CY, CL, Cl, Cm, Cn = mapreduce(sum, vcat, surface_coefficients(system; axes = Wind()))

    (CDi = CDi, CY = CY, CL = CL, Cl = Cl, Cm = Cm, Cn = Cn)
end

"""
    farfield_forces(system :: VortexLatticeSystem)

Compute the **farfield** forces in **wind axes** for all components of the `VortexLatticeSystem`.
"""
function farfield_forces(system :: VortexLatticeSystem)
    hs = system.vortices
    Γs = system.circulations
    α  = system.freestream.alpha
    β  = system.freestream.beta
    V  = system.reference.speed
    ρ  = system.reference.density

    NamedTuple(key => farfield_forces(Γs[key], hs[key], V, α, β, ρ) for key in keys(hs))
end

"""
    farfield_coefficients(system :: VortexLatticeSystem)

Compute the **total farfield** force coefficients in **wind axes** for all components of the `VortexLatticeSystem`.
"""
farfield_coefficients(system :: VortexLatticeSystem) = map(ff -> force_coefficient(ff, dynamic_pressure(system.reference.density, system.reference.speed), system.reference.area), farfield_forces(system))

"""
    farfield_coefficients(system :: VortexLatticeSystem)

Compute the **total farfield** force in **wind axes** for all components of the `VortexLatticeSystem`.
"""
function farfield(system :: VortexLatticeSystem)
    ffs = farfield_coefficients(system)
    # Massive hack
    if length(ffs) == 1
        coeffs = ffs[1]
    else
        coeffs = vec(sum(reduce(hcat, ffs), dims = 2))
    end

    return (CDi_ff = coeffs[1], CY_ff = coeffs[2], CL_ff = coeffs[3])
end


text:
module Qiskit

using PyCall

const qiskit = PyNULL()

function init()
    try
	copy!(qiskit, pyimport("qiskit"))
    catch err
        error("The Python Qiskit package could not be imported. You must install Python Qiskit before using this package. The error was $err")
    end
end

include("Terra.jl")

function execute(circuit, backend)
    qiskit.execute(circuit, backend)
end

function execute(circuit, backend, shots)
    qiskit.execute(circuit, backend, shots=1024)
end

include("BasicAer.jl")

include("Aer.jl")

include("Visualization.jl")

include("Version.jl")

end # module


text:
f(x) = x^3 - 3*x^2 + 2*x
fp(x) = 3*x^2-6*x+2

x_s, count = newton(f,fp,1.00)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode