text:
################################################################################
#
#   module LatticePhysics_LatticeModification
#   -> LatticePhysics_Base
#   -> LinearAlgebra
#
#   --> MODIFICATION OF LATTICES COMPONENTS
#       - Adding / Removing Connections
#       - Adding / Removing sites
#
#   --> RELABELING BONDS / SITES
#       - relabeling bonds
#       - relabeling sites
#
#   --> SPATIAL MODIFICATIONS
#       - Scaling  in space
#       - Rotating in space
#       - Shifting in space
#       - Lattice vector optimization (relabeling unitcells)
#
################################################################################

# Start of module
module LatPhysLatticeModification

# Using libraries
using LatPhysBase

# ADDING / REMOVING SITES
include("site_add_remove.jl")

# ADDING / REMOVING BONDS
include("bond_add_remove.jl")

# TRANSFORMATIONS IN REAL SPACE
include("real_space_transformation.jl")

# RELABELING

# setting labels of sites to bipartite
include("labeling_bipartite.jl")

# relabel all bonds to string
include("labeling_to_string.jl")

# relabel all sties to int
include("labeling_to_int.jl")

# End of module
end


text:
function cstrt_bezdek_ngon(n)
if mod(log2(n),1)==0 && n>=16
    # initialization
    v(u) = 4*pi/n-3*u;
    d(u) = (2*cos(2*u+v(u))+1)/(2*cos(u)+cos(3*u+v(u)));
    F(u) = 1+d(u)^2-2*d(u)*cos(u) - (2-2*cos(v(u)));
    u0 = find_zero(F,(pi/n, 4*pi/(3*n)));
    v0 = 4*pi/n-3*u0;
    d0 = (2*cos(2*u0+v0)+1)/(2*cos(u0)+cos(3*u0+v0));
    r = ones(Int(3*n/8));
    for i = 1:Int(3*n/16)
        r[2*i] = d0;
    end
    x = zeros(Int(3*n/8));
    for i = 1:Int(3*n/8)
        if mod(i,3) == 2
            x[i] = u0+v0;
        else
            x[i] = u0;
        end
    end
    p = zeros(Int(3*n/8-1)); q = zeros(Int(3*n/8-1));
    p[1] = r[1]*sin(x[1]); q[1] = r[1]*cos(x[1]);
    for i = 2:Int(3*n/8-1)
        p[i] = p[i-1] - (-1)^i*r[i]*sin(sum(x[1:i]));
        q[i] = q[i-1] - (-1)^i*r[i]*cos(sum(x[1:i]));
    end
    # construction
    a = zeros(n-1); b = zeros(n-1);
    a[Int(n/2-1)] = p[1]; a[Int(n/2+1)] = -a[Int(n/2-1)];
    b[Int(n/2-1)] = q[1]; b[Int(n/2+1)] = b[Int(n/2-1)];
    a[n-1] = a[Int(n/2-1)]-sin(x[1]+v0); a[1] = -a[n-1];
    b[n-1] = b[Int(n/2-1)]-cos(x[1]+v0); b[1] = b[n-1];
    a[n-2] = p[2]; a[2] = -a[n-2];
    b[n-2] = q[2]; b[2] = b[n-2];
    for i = 1:Int(n/8-1)
        a[Int(mod(n/2-1+i*(3*n/2-2),n))] = p[3*i+1];
        a[Int(n-mod(n/2-1+i*(3*n/2-2),n))] = -a[Int(mod(n/2-1+i*(3*n/2-2),n))];
        b[Int(mod(n/2-1+i*(3*n/2-2),n))] = q[3*i+1];
        b[Int(n-mod(n/2-1+i*(3*n/2-2),n))] = b[Int(mod(n/2-1+i*(3*n/2-2),n))];
        a[Int(mod(n-1+i*(3*n/2-2),n))] = a[Int(mod(n/2-1+i*(3*n/2-2),n))]-(-1)^i*sin(sum(x[1:3*i+1])+u0);
        a[Int(n-mod(n-1+i*(3*n/2-2),n))] = -a[Int(mod(n-1+i*(3*n/2-2),n))];
        b[Int(mod(n-1+i*(3*n/2-2),n))] = b[Int(mod(n/2-1+i*(3*n/2-2),n))]-(-1)^i*cos(sum(x[1:3*i+1])+u0);
        b[Int(n-mod(n-1+i*(3*n/2-2),n))] = b[Int(mod(n-1+i*(3*n/2-2),n))];
        a[Int(mod(n-2+i*(3*n/2-2),n))] = p[3*i+2];
        a[Int(n-mod(n-2+i*(3*n/2-2),n))] = -a[Int(mod(n-2+i*(3*n/2-2),n))];
        b[Int(mod(n-2+i*(3*n/2-2),n))] = q[3*i+2];
        b[Int(n-mod(n-2+i*(3*n/2-2),n))] = b[Int(mod(n-2+i*(3*n/2-2),n))];
        a[Int(mod(n+i*(3*n/2-2),n))] = p[3*i];
        a[Int(n-mod(n+i*(3*n/2-2),n))] = -a[Int(mod(n+i*(3*n/2-2),n))];
        b[Int(mod(n+i*(3*n/2-2),n))] = q[3*i];
        b[Int(n-mod(n+i*(3*n/2-2),n))] = b[Int(mod(n+i*(3*n/2-2),n))];
    end
    a[Int(n/2)] = 0; b[Int(n/2)] = d0;
    return(a,b)
end
end

ext: jl

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
using Mocking

accept_handshake_1_patch = @patch function accept_handshake(a...; b...)
    return Dict{String, Any}("Handshake" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "20170214T1215Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION"), Dict{String, Any}("Id" => "juan@example.com", "Type" => "EMAIL")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@amazon.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Org Master Account", "Type" => "MASTER_NAME"), Dict{String, Any}("Value" => "ALL", "Type" => "ORGANIZATION_FEATURE_SET")]), Dict{String, Any}("Value" => "222222222222", "Type" => "ACCOUNT")], "ExpirationTimestamp" => "20170228T1215Z", "Action" => "INVITE", "State" => "ACCEPTED"))
end

cancel_handshake_1_patch = @patch function cancel_handshake(a...; b...)
    return Dict{String, Any}("Handshake" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "20170214T1215Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION"), Dict{String, Any}("Id" => "susan@example.com", "Type" => "EMAIL")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@example.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Master Account", "Type" => "MASTER_NAME"), Dict{String, Any}("Value" => "CONSOLIDATED_BILLING", "Type" => "ORGANIZATION_FEATURE_SET")]), Dict{String, Any}("Value" => "222222222222", "Type" => "ACCOUNT"), Dict{String, Any}("Value" => "This is a request for Susan's account to join Bob's organization.", "Type" => "NOTES")], "ExpirationTimestamp" => "20170228T1215Z", "Action" => "INVITE", "State" => "CANCELED"))
end

create_account_1_patch = @patch function create_account(a...; b...)
    return Dict{String, Any}("CreateAccountStatus" => Dict{String, Any}("Id" => "car-examplecreateaccountrequestid111", "State" => "IN_PROGRESS"))
end

create_organization_1_patch = @patch function create_organization(a...; b...)
    return Dict{String, Any}("Organization" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:organization/o-exampleorgid", "Id" => "o-exampleorgid", "MasterAccountId" => "111111111111", "AvailablePolicyTypes" => Any[Dict{String, Any}("Status" => "ENABLED", "Type" => "SERVICE_CONTROL_POLICY")], "MasterAccountEmail" => "bill@example.com", "FeatureSet" => "ALL", "MasterAccountArn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/111111111111"))
end

create_organization_2_patch = @patch function create_organization(a...; b...)
    return Dict{String, Any}("Organization" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:organization/o-exampleorgid", "Id" => "o-exampleorgid", "MasterAccountId" => "111111111111", "AvailablePolicyTypes" => Any[], "MasterAccountEmail" => "bill@example.com", "FeatureSet" => "CONSOLIDATED_BILLING", "MasterAccountArn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/111111111111"))
end

create_organizational_unit_1_patch = @patch function create_organizational_unit(a...; b...)
    return Dict{String, Any}("OrganizationalUnit" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:ou/o-exampleorgid/ou-examplerootid111-exampleouid111", "Id" => "ou-examplerootid111-exampleouid111", "Name" => "AccountingOU"))
end

create_policy_1_patch = @patch function create_policy(a...; b...)
    return Dict{String, Any}("Policy" => Dict{String, Any}("Content" => "{\"Version\":\"2012-10-17\",\"Statement\":{\"Effect\":\"Allow\",\"Action\":\"s3:*\"}}", "PolicySummary" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:policy/o-exampleorgid/service_control_policy/p-examplepolicyid111", "Description" => "Allows delegation of all S3 actions", "Type" => "SERVICE_CONTROL_POLICY", "Name" => "AllowAllS3Actions")))
end

decline_handshake_1_patch = @patch function decline_handshake(a...; b...)
    return Dict{String, Any}("Handshake" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "2016-11-30T19:27:58Z", "Parties" => Any[Dict{String, Any}("Id" => "222222222222", "Type" => "ACCOUNT"), Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@example.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Master Account", "Type" => "MASTER_NAME")]), Dict{String, Any}("Value" => "222222222222", "Type" => "ACCOUNT"), Dict{String, Any}("Value" => "This is an invitation to Susan's account to join the Bill's organization.", "Type" => "NOTES")], "ExpirationTimestamp" => "2016-12-15T19:27:58Z", "Action" => "INVITE", "State" => "DECLINED"))
end

describe_account_1_patch = @patch function describe_account(a...; b...)
    return Dict{String, Any}("Account" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/555555555555", "Id" => "555555555555", "Email" => "anika@example.com", "Name" => "Beta Account"))
end

describe_create_account_status_1_patch = @patch function describe_create_account_status(a...; b...)
    return Dict{String, Any}("CreateAccountStatus" => Dict{String, Any}("Id" => "car-exampleaccountcreationrequestid", "AccountId" => "333333333333", "State" => "SUCCEEDED"))
end

describe_handshake_1_patch = @patch function describe_handshake(a...; b...)
    return Dict{String, Any}("Handshake" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "2016-11-30T17:24:58.046Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION"), Dict{String, Any}("Id" => "333333333333", "Type" => "ACCOUNT")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@example.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Master Account", "Type" => "MASTER_NAME")]), Dict{String, Any}("Value" => "333333333333", "Type" => "ACCOUNT")], "ExpirationTimestamp" => "2016-11-30T17:24:58.046Z", "Action" => "INVITE", "State" => "OPEN"))
end

describe_organization_1_patch = @patch function describe_organization(a...; b...)
    return Dict{String, Any}("Organization" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:organization/o-exampleorgid", "Id" => "o-exampleorgid", "AvailablePolicyTypes" => Any[Dict{String, Any}("Status" => "ENABLED", "Type" => "SERVICE_CONTROL_POLICY")], "MasterAccountEmail" => "bill@example.com", "FeatureSet" => "ALL", "MasterAccountArn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/111111111111"))
end

describe_organizational_unit_1_patch = @patch function describe_organizational_unit(a...; b...)
    return Dict{String, Any}("OrganizationalUnit" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:ou/o-exampleorgid/ou-examplerootid111-exampleouid111", "Id" => "ou-examplerootid111-exampleouid111", "Name" => "Accounting Group"))
end

describe_policy_1_patch = @patch function describe_policy(a...; b...)
    return Dict{String, Any}("Policy" => Dict{String, Any}("Content" => "{\\n  \\\"Version\\\": \\\"2012-10-17\\\",\\n  \\\"Statement\\\": [\\n    {\\n      \\\"Effect\\\": \\\"Allow\\\",\\n      \\\"Action\\\": \\\"*\\\",\\n      \\\"Resource\\\": \\\"*\\\"\\n    }\\n  ]\\n}", "PolicySummary" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:policy/o-exampleorgid/service_control_policy/p-examplepolicyid111", "Id" => "p-examplepolicyid111", "Description" => "Enables admins to delegate S3 permissions", "AwsManaged" => false, "Type" => "SERVICE_CONTROL_POLICY", "Name" => "AllowAllS3Actions")))
end

disable_policy_type_1_patch = @patch function disable_policy_type(a...; b...)
    return Dict{String, Any}("Root" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:root/o-exampleorgid/r-examplerootid111", "Id" => "r-examplerootid111", "PolicyTypes" => Any[], "Name" => "Root"))
end

enable_all_features_1_patch = @patch function enable_all_features(a...; b...)
    return Dict{String, Any}("Handshake" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/enable_all_features/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "2017-02-13T09:35:40.05Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION")], "ExpirationTimestamp" => "2017-02-28T09:35:40.05Z", "Action" => "ENABLE_ALL_FEATURES", "State" => "REQUESTED"))
end

enable_policy_type_1_patch = @patch function enable_policy_type(a...; b...)
    return Dict{String, Any}("Root" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:root/o-exampleorgid/r-examplerootid111", "Id" => "r-examplerootid111", "PolicyTypes" => Any[Dict{String, Any}("Status" => "ENABLED", "Type" => "SERVICE_CONTROL_POLICY")], "Name" => "Root"))
end

invite_account_to_organization_1_patch = @patch function invite_account_to_organization(a...; b...)
    return Dict{String, Any}("Handshake" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "2017-02-01T09:36:05.02Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION"), Dict{String, Any}("Id" => "juan@example.com", "Type" => "EMAIL")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@amazon.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Org Master Account", "Type" => "MASTER_NAME"), Dict{String, Any}("Value" => "FULL", "Type" => "ORGANIZATION_FEATURE_SET")]), Dict{String, Any}("Value" => "juan@example.com", "Type" => "EMAIL")], "ExpirationTimestamp" => "2017-02-16T09:36:05.02Z", "Action" => "INVITE", "State" => "OPEN"))
end

list_accounts_1_patch = @patch function list_accounts(a...; b...)
    return Dict{String, Any}("Accounts" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/111111111111", "Id" => "111111111111", "Email" => "bill@example.com", "JoinedMethod" => "INVITED", "Status" => "ACTIVE", "JoinedTimestamp" => "20161215T193015Z", "Name" => "Master Account"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/222222222222", "Id" => "222222222222", "Email" => "alice@example.com", "JoinedMethod" => "INVITED", "Status" => "ACTIVE", "JoinedTimestamp" => "20161215T210221Z", "Name" => "Developer Account"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/333333333333", "Id" => "333333333333", "Email" => "juan@example.com", "JoinedMethod" => "INVITED", "Status" => "ACTIVE", "JoinedTimestamp" => "20161215T210347Z", "Name" => "Test Account"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/444444444444", "Id" => "444444444444", "Email" => "anika@example.com", "JoinedMethod" => "INVITED", "Status" => "ACTIVE", "JoinedTimestamp" => "20161215T210332Z", "Name" => "Production Account")])
end

list_accounts_for_parent_1_patch = @patch function list_accounts_for_parent(a...; b...)
    return Dict{String, Any}("Accounts" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/333333333333", "Id" => "333333333333", "Email" => "juan@example.com", "JoinedMethod" => "INVITED", "Status" => "ACTIVE", "JoinedTimestamp" => 1.481835795536e9, "Name" => "Development Account"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/444444444444", "Id" => "444444444444", "Email" => "anika@example.com", "JoinedMethod" => "INVITED", "Status" => "ACTIVE", "JoinedTimestamp" => 1.481835812143e9, "Name" => "Test Account")])
end

list_children_1_patch = @patch function list_children(a...; b...)
    return Dict{String, Any}("Children" => Any[Dict{String, Any}("Id" => "ou-examplerootid111-exampleouid111", "Type" => "ORGANIZATIONAL_UNIT"), Dict{String, Any}("Id" => "ou-examplerootid111-exampleouid222", "Type" => "ORGANIZATIONAL_UNIT")])
end

list_create_account_status_1_patch = @patch function list_create_account_status(a...; b...)
    return Dict{String, Any}("CreateAccountStatuses" => Any[Dict{String, Any}("Id" => "car-exampleaccountcreationrequestid1", "RequestedTimestamp" => "2017-01-15T13:45:23.01Z", "AccountId" => "444444444444", "State" => "SUCCEEDED", "CompletedTimestamp" => "2017-01-15T13:45:23.6Z", "AccountName" => "Developer Test Account")])
end

list_create_account_status_2_patch = @patch function list_create_account_status(a...; b...)
    return Dict{String, Any}("CreateAccountStatuses" => Any[Dict{String, Any}("Id" => "car-exampleaccountcreationrequestid2", "RequestedTimestamp" => "2017-01-15T13:45:23.01Z", "State" => "IN_PROGRESS", "AccountName" => "Production Account")])
end

list_handshakes_for_account_1_patch = @patch function list_handshakes_for_account(a...; b...)
    return Dict{String, Any}("Handshakes" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "2017-01-13T14:35:23.3Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION"), Dict{String, Any}("Id" => "juan@example.com", "Type" => "EMAIL")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@amazon.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Org Master Account", "Type" => "MASTER_NAME"), Dict{String, Any}("Value" => "FULL", "Type" => "ORGANIZATION_FEATURE_SET")]), Dict{String, Any}("Value" => "juan@example.com", "Type" => "EMAIL")], "ExpirationTimestamp" => "2017-01-28T14:35:23.3Z", "Action" => "INVITE", "State" => "OPEN")])
end

list_handshakes_for_organization_1_patch = @patch function list_handshakes_for_organization(a...; b...)
    return Dict{String, Any}("Handshakes" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid111", "RequestedTimestamp" => "2017-01-13T14:35:23.3Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION"), Dict{String, Any}("Id" => "juan@example.com", "Type" => "EMAIL")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@amazon.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Org Master Account", "Type" => "MASTER_NAME"), Dict{String, Any}("Value" => "FULL", "Type" => "ORGANIZATION_FEATURE_SET")]), Dict{String, Any}("Value" => "juan@example.com", "Type" => "EMAIL")], "ExpirationTimestamp" => "2017-01-28T14:35:23.3Z", "Action" => "INVITE", "State" => "OPEN"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:handshake/o-exampleorgid/invite/h-examplehandshakeid111", "Id" => "h-examplehandshakeid222", "RequestedTimestamp" => "2017-01-13T14:35:23.3Z", "Parties" => Any[Dict{String, Any}("Id" => "o-exampleorgid", "Type" => "ORGANIZATION"), Dict{String, Any}("Id" => "anika@example.com", "Type" => "EMAIL")], "Resources" => Any[Dict{String, Any}("Value" => "o-exampleorgid", "Type" => "ORGANIZATION", "Resources" => Any[Dict{String, Any}("Value" => "bill@example.com", "Type" => "MASTER_EMAIL"), Dict{String, Any}("Value" => "Master Account", "Type" => "MASTER_NAME")]), Dict{String, Any}("Value" => "anika@example.com", "Type" => "EMAIL"), Dict{String, Any}("Value" => "This is an invitation to Anika's account to join Bill's organization.", "Type" => "NOTES")], "ExpirationTimestamp" => "2017-01-28T14:35:23.3Z", "Action" => "INVITE", "State" => "ACCEPTED")])
end

list_organizational_units_for_parent_1_patch = @patch function list_organizational_units_for_parent(a...; b...)
    return Dict{String, Any}("OrganizationalUnits" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:ou/o-exampleorgid/ou-examlerootid111-exampleouid111", "Id" => "ou-examplerootid111-exampleouid111", "Name" => "Development"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:ou/o-exampleorgid/ou-examlerootid111-exampleouid222", "Id" => "ou-examplerootid111-exampleouid222", "Name" => "Production")])
end

list_parents_1_patch = @patch function list_parents(a...; b...)
    return Dict{String, Any}("Parents" => Any[Dict{String, Any}("Id" => "ou-examplerootid111-exampleouid111", "Type" => "ORGANIZATIONAL_UNIT")])
end

list_policies_1_patch = @patch function list_policies(a...; b...)
    return Dict{String, Any}("Policies" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:policy/o-exampleorgid/service_control_policy/p-examplepolicyid111", "Id" => "p-examplepolicyid111", "Description" => "Enables account admins to delegate permissions for any S3 actions to users and roles in their accounts.", "AwsManaged" => false, "Type" => "SERVICE_CONTROL_POLICY", "Name" => "AllowAllS3Actions"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:policy/o-exampleorgid/service_control_policy/p-examplepolicyid222", "Id" => "p-examplepolicyid222", "Description" => "Enables account admins to delegate permissions for any EC2 actions to users and roles in their accounts.", "AwsManaged" => false, "Type" => "SERVICE_CONTROL_POLICY", "Name" => "AllowAllEC2Actions"), Dict{String, Any}("Arn" => "arn:aws:organizations::aws:policy/service_control_policy/p-FullAWSAccess", "Id" => "p-FullAWSAccess", "Description" => "Allows access to every operation", "AwsManaged" => true, "Type" => "SERVICE_CONTROL_POLICY", "Name" => "FullAWSAccess")])
end

list_policies_for_target_1_patch = @patch function list_policies_for_target(a...; b...)
    return Dict{String, Any}("Policies" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:policy/o-exampleorgid/service_control_policy/p-examplepolicyid222", "Id" => "p-examplepolicyid222", "Description" => "Enables account admins to delegate permissions for any EC2 actions to users and roles in their accounts.", "AwsManaged" => false, "Type" => "SERVICE_CONTROL_POLICY", "Name" => "AllowAllEC2Actions")])
end

list_roots_1_patch = @patch function list_roots(a...; b...)
    return Dict{String, Any}("Roots" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:root/o-exampleorgid/r-examplerootid111", "Id" => "r-examplerootid111", "PolicyTypes" => Any[Dict{String, Any}("Status" => "ENABLED", "Type" => "SERVICE_CONTROL_POLICY")], "Name" => "Root")])
end

list_targets_for_policy_1_patch = @patch function list_targets_for_policy(a...; b...)
    return Dict{String, Any}("Targets" => Any[Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:root/o-exampleorgid/r-examplerootid111", "Type" => "ROOT", "Name" => "Root", "TargetId" => "r-examplerootid111"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:account/o-exampleorgid/333333333333;", "Type" => "ACCOUNT", "Name" => "Developer Test Account", "TargetId" => "333333333333"), Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:ou/o-exampleorgid/ou-examplerootid111-exampleouid111", "Type" => "ORGANIZATIONAL_UNIT", "Name" => "Accounting", "TargetId" => "ou-examplerootid111-exampleouid111")])
end

update_organizational_unit_1_patch = @patch function update_organizational_unit(a...; b...)
    return Dict{String, Any}("OrganizationalUnit" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:ou/o-exampleorgid/ou-examplerootid111-exampleouid111", "Id" => "ou-examplerootid111-exampleouid111", "Name" => "AccountingOU"))
end

update_policy_1_patch = @patch function update_policy(a...; b...)
    return Dict{String, Any}("Policy" => Dict{String, Any}("Content" => "{ \"Version\": \"2012-10-17\", \"Statement\": { \"Effect\": \"Allow\", \"Action\": \"ec2:*\", \"Resource\": \"*\" } }", "PolicySummary" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:policy/o-exampleorgid/service_control_policy/p-examplepolicyid111", "Id" => "p-examplepolicyid111", "Description" => "This description replaces the original.", "AwsManaged" => false, "Type" => "SERVICE_CONTROL_POLICY", "Name" => "Renamed-Policy")))
end

update_policy_2_patch = @patch function update_policy(a...; b...)
    return Dict{String, Any}("Policy" => Dict{String, Any}("Content" => "{ \\\"Version\\\": \\\"2012-10-17\\\", \\\"Statement\\\": { \\\"Effect\\\": \\\"Allow\\\", \\\"Action\\\": \\\"s3:*\\\", \\\"Resource\\\": \\\"*\\\" } }", "PolicySummary" => Dict{String, Any}("Arn" => "arn:aws:organizations::111111111111:policy/o-exampleorgid/service_control_policy/p-examplepolicyid111", "Id" => "p-examplepolicyid111", "Description" => "This description replaces the original.", "AwsManaged" => false, "Type" => "SERVICE_CONTROL_POLICY", "Name" => "Renamed-Policy")))
end


text:
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/build_literate/Analyze_NMR_PreSat_Data.ipynb)

# # Continuous Wave Saturation Experiments
# The following code analyzes data from a steady-state experiment similar to the original work of [Henkelman et al](https://doi.org/10.1002/mrm.1910290607). In this experiment, the magnetization of the coupled spin system is saturated with off-resonant continuous waves of the exponentially spaced frequencies:
Δ = exp.(range(log(0.01e3), log(100e3), length=20)) * 2π # rad/s
# and the amplitudes:
ω1_dB = -60:5:-5 # dB
ω1 = @. 10^(ω1_dB / 20) * π / 2 / 11.4e-6 # rad/s
# The waves were applied for 7 seconds to ensure a steady state. Thereafter, the magnetization was excited with a π/2-pulse and an FID was acquired. The repetition times was 30s to ensure full recovery to thermal equilibrium.

# We fit the data with Henkelman's closed form solution to this steady-state problem while assuming a Lorentzian lineshape for the free spin pool, and different lineshapes for the semi-solid spin pool:
g_Lorentzian(Δ, T2) = T2 / π / (1 + (T2 * Δ)^2)
g_Gaussian(Δ, T2) = T2 / sqrt(2π) * exp(-(T2 * Δ)^2 / 2)
g_superLorentzian(Δ, T2) = T2 * sqrt(2 / π) * quadgk(ct -> exp(- 2 *  (T2 * Δ / (3 * ct^2 - 1))^2) / abs(3 * ct^2 - 1), 0, sqrt(1 / 3), 1)[1];

# For this data analysis we need the following packages:
using MRIgeneralizedBloch
using QuadGK
using LsqFit
using LinearAlgebra
using Statistics
using Printf
using Plots
plotlyjs(bg=RGBA(31 / 255, 36 / 255, 36 / 255, 1.0), ticks=:native); # hide #!nb

# The raw data is stored in a separate [github repository](https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData) and the following functions return the URL to the individual files:
MnCl2_data(ω1_dB) = string("https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData/blob/main/20210419_1mM_MnCl2/ja_PreSat_v2%20(", ω1_dB, ")/1/data.2d?raw=true")
BSA_data(ω1_dB)   = string("https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData/blob/main/20210416_15%25BSA_2ndBatch/ja_PreSat_v2%20(", ω1_dB, ")/1/data.2d?raw=true");

# which can be loaded with functions implemented in this file:
include(string(pathof(MRIgeneralizedBloch), "/../../docs/src/load_NMR_data.jl"));

# We store the off-resonance frequencies and wave amplitudes in a single matrix for convenience:
x = zeros(Float64, length(ω1) * length(Δ), 2)
x[:,1] = repeat(Δ, length(ω1))
x[:,2] = vec(repeat(ω1, 1, length(Δ))');

# ## MnCl``_2`` Sample
# We load the first data point of each FID:
M = zeros(Float64, length(Δ), length(ω1))
for i = 1:length(ω1_dB)
    M[:,i] = load_first_datapoint(MnCl2_data(ω1_dB[i]); set_phase=:abs)
end
M ./= maximum(M);
# In contrast to the inversion-recovery experiment, the phase of the signal was not stable. Therefore, we took the absolute value of the signal by setting the flag `set_phase=:abs`.

# The MnCl``_2``-data can be described with a single compartment model:
function single_compartment_model(x, p)
    (m0, R1, T2) = p

    Δ  = @view x[:,1]
    ω1 = @view x[:,2]

    Rrf = @. π * ω1^2 * g_Lorentzian(Δ, T2)
    m = @. m0 * R1 / (R1 + Rrf)
    return m
end;
# (cf. Eqs. (14) and (15) in the [paper](https://arxiv.org/pdf/2107.11000.pdf)).

# As this model is merely a function of the relaxation times T₁ and T₂, we forgo a fitting routine and use the estimates from the [Inversion Recovery Experiments](@ref) instead:
R1 = 1.479 # 1/s
T2 = 0.075 # s
nothing # hide #md

# Visually, this model describes the data well:
p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(single_compartment_model(x, [1,R1,T2]), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]
display(p) #!md
#md Main.HTMLPlot(p) #hide

# ## Bovine Serum Albumin Sample
# We acquired the same data for the BSA sample, which we load:
for i = 1:length(ω1_dB)
    M[:,i] = load_first_datapoint(BSA_data(ω1_dB[i]); set_phase=:abs)
end
M ./= maximum(M);

# We model the steady-state magnetization as described by [Henkelman et al.](https://doi.org/10.1002/mrm.1910290607):
function Henkelman_model(x, p; lineshape=:superLorentzian)
    (m0, m0s, R1f, R1s, T2f, T2s, Rx) = p

    m0s /= 1 - m0s # switch from m0s + m0f = 1 to m0f = 1 normalization

    Δ  = @view x[:,1]
    ω1 = @view x[:,2]

    Rrf_f = @. π * ω1^2 * g_Lorentzian(Δ, T2f)

    if lineshape == :Lorentzian
        Rrf_s = @. π * ω1^2 * g_Lorentzian(Δ, T2s)
    elseif lineshape == :Gaussian
        Rrf_s = @. π * ω1^2 * g_Gaussian(Δ, T2s)
    elseif lineshape == :superLorentzian
        Rrf_s = @. π * ω1^2 * g_superLorentzian(Δ, T2s)
    end

    m = @. m0 * (R1s * Rx * m0s + Rrf_s * R1f + R1f * R1s + R1f * Rx) / ((R1f + Rrf_f + Rx * m0s) * (R1s + Rrf_s + Rx) - Rx^2 * m0s)
    return m
end;

# Here, we use a fitting routine to demonstrate the best possible fit with each of the three lineshapes. We define an initialization for the fitting routine `p0 = [m0, m0s, R1f, R1s, T2f, T2s, Rx]` and set some reasonable bounds:
p0   = [  1,0.01,   1,   5,0.052,  1e-5, 40]
pmin = [  0,   0,   0,   0,0.052,  1e-6,  1]
pmax = [Inf,   1, Inf, Inf,0.052, 10e-3,100];
# Note that we fixed T₂ᶠ = 52ms to the value estimated with the [Inversion Recovery Experiments](@ref) as T₂ᶠ is poorly defined by this saturation experiment.

# ### Super-Lorentzian Lineshape
# Fitting the model with a super-Lorentzian lineshape to the data achieves good concordance:
fit = curve_fit((x, p) -> Henkelman_model(x, p; lineshape=:superLorentzian), x, vec(M), p0, lower=pmin, upper=pmax)
fit_std = stderror(fit)

p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(Henkelman_model(x, fit.param), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]
display(p) #!md
#md Main.HTMLPlot(p) #hide

# ### Lorentzian Lineshape
# The Lorentzian lineshape, on the other hand, does not fit the data well:
fit = curve_fit((x, p) -> Henkelman_model(x, p; lineshape=:Lorentzian), x, vec(M), p0, lower=pmin, upper=pmax)
fit_std = stderror(fit)

p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(Henkelman_model(x, fit.param), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]
display(p) #!md
#md Main.HTMLPlot(p) #hide

# ### Gaussian Lineshape
# And the Gaussian lineshape does not not fit the data well either:
fit = curve_fit((x, p) -> Henkelman_model(x, p; lineshape=:Lorentzian), x, vec(M), p0, lower=pmin, upper=pmax)
fit_std = stderror(fit)

p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(Henkelman_model(x, fit.param), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]
display(p) #!md
#md Main.HTMLPlot(p) #hide


text:
#### SolveTurbConv

function solve_tridiag_wrapper!(
    grid::Grid,
    sv::StateVec,
    ϕ::Symbol,
    i::Int,
    tri_diag::StateVec,
)
    f = [tri_diag[:f, k] for k in over_elems_real(grid)]
    a = [tri_diag[:a, k] for k in over_elems_real(grid)[2:end]]
    b = [tri_diag[:b, k] for k in over_elems_real(grid)]
    c = [tri_diag[:c, k] for k in over_elems_real(grid)[1:(end - 1)]]
    A = Tridiagonal(a, b, c)
    x = inv(A) * f
    assign_real!(sv, ϕ, grid, x, i)
end

function solve!(
    q_new::StateVec,
    q::StateVec,
    q_tendencies::StateVec,
    ∑tendencies_params,
)
    @unpack params, grid, aux_O2, aux, tri_diag, case = ∑tendencies_params
    gm, en, ud, sd, al = allcombinations(q)

    compute_new_ud_a!(grid, q_new, q, q_tendencies, aux, params)
    apply_bcs!(grid, q_new, aux, params, case)

    compute_new_ud_w!(grid, q_new, q, q_tendencies, aux, params)
    compute_new_ud_scalars!(grid, q_new, q, q_tendencies, aux, params)

    apply_bcs!(grid, q_new, aux, params, case)

    solve_tridiag_wrapper!(grid, q_new, :tke, en, tri_diag[:tke])
    solve_tridiag_wrapper!(grid, q_new, :q_tot, gm, tri_diag[:q_tot])
    solve_tridiag_wrapper!(grid, q_new, :θ_liq, gm, tri_diag[:θ_liq])

    assign_values_to_new!(grid, q, q_new, aux)
    apply_bcs!(grid, q, aux, params, case)
    extrap_0th_order!(q, (:θ_liq, :q_tot), grid, gm)

end


text:
using ZipFile
resultspath="/home/ubuntu/data/dat"
resultspathagg="/home/ubuntu/data/agg"

const   agents = Dict{Int64,String}()
const   nodes = Dict{Int64,Array{Int64,1}}()
const   coords = Dict{Int64,Tuple{Float64,Float64}}()
const   routes1 = Dict{Int64,String}()
const   routes2 = Dict{Int64,String}()

const   routes1coords = Dict{Int64,String}()
const   routes2coords = Dict{Int64,String}()

function split_line(line::AbstractString, delim::AbstractString, firstc::Int64)
    res1=String[]
    ix = 1
    while length(res1) < firstc
        ix2 = findnext(delim,line,ix)[1]
        push!(res1,line[ix:(ix2-1)])
        ix = ix2+1
    end
    last_item = line[(findlast(";",line)[1]+1):end]
    (res1,line[ix:end],last_item)
end

function aggregate2!(resultspath,agents,nodes,coords)
    local lineheader_r
    f_count = 0
    for ff in readdir(resultspath)
        global ffn = joinpath(resultspath,ff)

        if isfile(ffn) && endswith(ff,"_nodes.csv.zip") && f_count < 30
            f_count += 1
            println("$f_count reading $ff")
            local line
            r = ZipFile.Reader(ffn);
            f = r.files[1]

            lineheader_s, lineheader_r = split_line(strip(readline(f)),";",4)
            while (!eof(f))
                line = strip(readline(f))
                line_s,line_r,line_last = split_line(line,";",4)
                id = parse(Int64,line_last)

                agents[id] = line_r
                node_id = parse(Int64,line_s[2])
                if !haskey(nodes,node_id)
                    nodes[node_id] = Int64[]
                    coords[node_id] = (parse(Float64,line_s[3]),parse(Float64,line_s[4]))
                end
                push!(nodes[node_id],id)
            end
            close(f)
            close(r)
        end
    end
    return lineheader_r
end

const node_line_header = aggregate2!(resultspath,agents,nodes,coords)

function collect_routes!(routes1,routes2,routes1coords,routes2coords)
    f_count = 0
    for ff in readdir(resultspath)
        ffn = joinpath(resultspath,ff)
        if isfile(ffn) && endswith(ff,"_routes.csv") && f_count < 30
            f_count += 1
            println("$f_count reading $ff")
            f = open(ffn,"r");
            local line
            while (!eof(f))
                line = strip(readline(f))
                if length(line) > 0
                    els = split(line,";")
                    agentid = parse(Int64,els[1])
                    nodes = parse.(Int64,els[3:end])
                    if els[2] == "towork"
                        routes1[agentid] = join(nodes,"#")
                        routes1coords[agentid] = join([join(coords[route], ",") for route in nodes],"#")
                    else
                        routes2[agentid] = join(nodes,"#")
                        routes2coords[agentid] = join([join(coords[route], ",") for route in nodes],"#")
                    end
                end
            end
            close(f)
        end
    end
end

collect_routes!(routes1,routes2,routes1coords,routes2coords)

function export_data(agents,nodes,coords,routes1,routes2,node_line_header)
    fns = open(joinpath(resultspathagg,"Node_Stats_.csv"),"w")
    println(fns,"\"NODE_ID\";\"longitude\";\"latitude\";\"vehicle_count\"")
    routes_tuple = (routes1,routes2)
    coords_tuple = (routes1coords,routes2coords)
    for node in sort!(collect(keys(nodes)))
        println(fns,"$(node);$(coords[node][1]);$(coords[node][2]);$(length(nodes[node]))")
        ftr = open(joinpath(resultspathagg,"Travel_$(node)_$(coords[node][1])_$(coords[node][2]).csv"),"w")
        fn = open(joinpath(resultspathagg,"Routes_n_$(node)_$(coords[node][1])_$(coords[node][2]).csv"),"w")
        fc = open(joinpath(resultspathagg,"Routes_c_$(node)_$(coords[node][1])_$(coords[node][2]).csv"),"w")
        println(ftr,"\"NODE_ID\";\"longitude\";\"latitude\";"*node_line_header)
        for agentid in nodes[node]
            println(ftr,"$(node);$(coords[node][1]);$(coords[node][2]);"*agents[agentid])
            routetype = "towork"
            for i in 1:2
                println(fn,"$agentid;$routetype;$(routes_tuple[i][agentid])")
                # dolozyc do linii pozynizej odpytywanie z coords
                #println(fc,"$agentid;$routetype;"*join([coords[route] for route in routes[agentid]],";"))
                println(fc, "$agentid;$routetype;$(coords_tuple[i][agentid])")
                routetype = "tohome"
            end
        end
        close(ftr)
        close(fn)
        close(fc)
        flush(fns)
    end
    close(fns)
end

export_data(agents,nodes,coords,routes1,routes2,node_line_header)


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# This is a short script that transforms an MPS file by dropping all
# integrality constraints.

import SCIP

# Note: The SCIP version bundled with Julia can't process .gz files.
function drop_integrality(input_file, output_file)

  scip = Ref{Ptr{SCIP.SCIP_}}(C_NULL)
  SCIP.@SC SCIP.SCIPcreate(scip)

  @assert scip[] != C_NULL
  SCIP.@SC SCIP.SCIPincludeDefaultPlugins(scip[])

  SCIP.@SC SCIP.SCIPreadProb(scip[], input_file, C_NULL)

  num_vars = SCIP.SCIPgetNVars(scip[])

  vars = Vector{Ptr{SCIP.SCIP_VAR}}(undef, num_vars)

  raw_vars = SCIP.SCIPgetVars(scip[])
  # Copy the vars to an array we own.
  # The SCIP documentation warns that SCIPchgVarType() can invalidate the result
  # of SCIPGetVars().
  vars .= unsafe_wrap(Array, raw_vars, num_vars)

  for var in vars
    # This is an output of chgVarType indicating whether the change makes the
    # problem infeasible. We ignore this value.
    infeasible = Ref{SCIP.SCIP_Bool}()
    SCIP.@SC SCIP.SCIPchgVarType(
      scip[],
      var,
      SCIP.SCIP_VARTYPE_CONTINUOUS,
      infeasible,
    )
  end
  @assert SCIP.SCIPgetNBinVars(scip[]) == 0
  @assert SCIP.SCIPgetNIntVars(scip[]) == 0

  SCIP.@SC SCIP.SCIPwriteOrigProblem(scip[], output_file, C_NULL, false)

  SCIP.@SC SCIP.SCIPfree(scip)
end

if length(ARGS) != 2
  @error "Usage: drop_integrality.jl input_file output_file"
end
drop_integrality(ARGS[1], ARGS[2])


text:
using Distributed

pids = addprocs(4)

@everywhere include((@__DIR__)*"/sharedarrays_setup.jl")

function shared_BandedBlockBandedMatrix(::UndefInitializer, bs::BlockSizes, (l, u), (λ, μ); pids=Int[])
    bs = BandedBlockBandedSizes(bs, l, u, λ, μ)
    data = SharedMatrix{Float64}((l+u+1)*(λ+μ+1), size(A,2); pids=pids)
    _BandedBlockBandedMatrix(PseudoBlockArray(data, bs.data_block_sizes), bs)
end

shared_BandedBlockBandedMatrix(::UndefInitializer, (N,M), (l, u), (λ, μ); pids=Int[]) =
    shared_BandedBlockBandedMatrix(undef, BlockSizes(N,M), (l, u), (λ, μ); pids=pids)

function populate!(ret, A::AbstractMatrix)
    pids = procs(ret.data.blocks)
    N,M = nblocks(A)
    bl = N ÷ length(pids)
    @sync begin
        for k = 1:length(pids)
            JR = Block.(((k-1)*bl+1):(k*bl))  #TODO: case where length(pids)  does not divide N
            @async remotecall_wait(blockcolbuild!, pids[k], ret, A, JR)
        end
    end
    ret
end

function shared_BandedBlockBandedMatrix(A::AbstractMatrix; pids=Int[])
    ret = shared_BandedBlockBandedMatrix(undef, blocksizes(A), (1,1), (0,0); pids=pids)
    populate!(ret, A)
end

N = 4 * 1000
    h = 1/N
    D² = FiniteDifference(N)
    A = Kron(D², Eye(N))
@time BandedBlockBandedMatrix(A)           # 5s
@time shared_BandedBlockBandedMatrix(A)    # 0.9s


text:
using Revise
using DistCtrl4DistMan: runExp
using Statistics
using StatsBase
using MAT

DEP_params = Dict{String,Any}("maxDist" => (550e-6, 550e-6));
MAG_params = Dict{String,Any}("maxDist" => (75e-3, 75e-3));
ACU_params = Dict{String,Any}("maxDist" => (50e-3, 50e-3));

MAG - distributed
print("++++++++++++++++++++++++ MAG ++++++++++++++++++++++")
N_agnts_vals = collect(5:5:35)
actuators_per_agent_vals = collect(6:1:12)
nagents2arrsize = (N_agnts, actuators_per_agent) -> (Int64(ceil(sqrt(actuators_per_agent*N_agnts))), Int64(ceil(sqrt(actuators_per_agent*N_agnts))))
N_exps = 200

t_elapsed_admm = zeros(length(N_agnts_vals), length(actuators_per_agent_vals), N_exps)
t_elapsed_admm_mean = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_admm_std = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_centralized = zeros(length(N_agnts_vals), length(actuators_per_agent_vals), N_exps)
t_elapsed_centralized_mean = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_centralized_std = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
for (k, N_agnts) in enumerate(N_agnts_vals)
    for (l, actuators_per_agent) in enumerate(actuators_per_agent_vals)
        N_acts = nagents2arrsize(N_agnts, actuators_per_agent)

        print(k, "/", length(N_agnts_vals), " # Agents ", N_agnts, " | ", l, "/", length(actuators_per_agent_vals)," agent density values\n")

        print("ADMM - ")
        params, exp_data, t_elapsed = runExp(platform=:MAG, algorithm=:admm, N_exps=N_exps, N_iter=1_000_000, N_agnts=N_agnts, N_acts=N_acts, convanalysis=true,
                                plotConvergenceRates=false, saveplots=false, params=MAG_params,
                                plotConvergenceRatesIndividual=true, errstats=false, stopping_value=1e-3, showplots=false, display=:final);
        t_elapsed_admm[k, l, :] = t_elapsed
        t_elapsed_admm_mean[k, l] = mean(trim(t_elapsed, prop=0.25))
        t_elapsed_admm_std[k, l] = std(trim(t_elapsed, prop=0.25))

        iters = [length(exp["conv_data"].expData) for exp in exp_data]
        print("N_iters: mean: ", mean(iters), ", std: ", std(iters), "\n")

        print("Centralized - ")
        params, exp_data, t_elapsed = runExp(platform=:MAG, algorithm=:centralized, N_exps=N_exps, N_iter=1_000_000, N_agnts=N_agnts, N_acts=N_acts, convanalysis=true,
                                plotConvergenceRates=false, saveplots=false, params=MAG_params,
                                plotConvergenceRatesIndividual=true, errstats=false, stopping_value=1e-3, showplots=false, display=:final);
        t_elapsed_centralized[k, l, :] = t_elapsed
        t_elapsed_centralized_mean[k, l] = mean(trim(t_elapsed, prop=0.25))
        t_elapsed_centralized_std[k, l] = std(trim(t_elapsed, prop=0.25))

        iters = [length(exp["conv_data"].expData) for exp in exp_data]
        print("N_iters: mean: ", mean(iters), ", std: ", std(iters), "\n\n")
    end
end
##
file = matopen("/Users/martingurtner/Downloads/MAG_conv_trim_M1.mat", "w")
write(file, "t_elapsed_admm", t_elapsed_admm)
write(file, "t_elapsed_admm_mean", t_elapsed_admm_mean)
write(file, "t_elapsed_admm_std", t_elapsed_admm_std)
write(file, "t_elapsed_centralized", t_elapsed_centralized)
write(file, "t_elapsed_centralized_mean", t_elapsed_centralized_mean)
write(file, "t_elapsed_centralized_std", t_elapsed_centralized_std)
write(file, "N_agnts_vals", N_agnts_vals)
write(file, "actuators_per_agent_vals", actuators_per_agent_vals)
close(file)

## ACU - distributed
print("++++++++++++++++++++++++ ACU ++++++++++++++++++++++")
N_agnts_vals = collect(20:20:120)
actuators_per_agent_vals = collect(40:20:120)
# actuators_per_agent_vals = collect(100:10:100)
nagents2arrsize = (N_agnts, actuators_per_agent) -> (Int64(ceil(sqrt(actuators_per_agent*N_agnts))), Int64(ceil(sqrt(actuators_per_agent*N_agnts))))
N_exps = 100

t_elapsed_admm = zeros(length(N_agnts_vals), length(actuators_per_agent_vals), N_exps)
t_elapsed_admm_mean = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_admm_std = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_centralized = zeros(length(N_agnts_vals), length(actuators_per_agent_vals), N_exps)
t_elapsed_centralized_mean = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_centralized_std = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
for (k, N_agnts) in enumerate(N_agnts_vals)
    for (l, actuators_per_agent) in enumerate(actuators_per_agent_vals)
        N_acts = nagents2arrsize(N_agnts, actuators_per_agent)

        print(k, "/", length(N_agnts_vals), " # Agents ", N_agnts, " | ", l, "/", length(actuators_per_agent_vals)," agent density values\n")

        print("ADMM - ")
        params, exp_data, t_elapsed = runExp(platform=:ACU, algorithm=:admm, N_exps=N_exps, N_iter=1_000_000, N_agnts=N_agnts, N_acts=N_acts, convanalysis=true,
                                plotConvergenceRates=false, saveplots=false, params=ACU_params,
                                plotConvergenceRatesIndividual=true, errstats=false, stopping_value=1e-3, showplots=false, display=:final);
        t_elapsed_admm[k, l, :] = t_elapsed
        t_elapsed_admm_mean[k, l] = mean(trim(t_elapsed, prop=0.25))
        t_elapsed_admm_std[k, l] = std(trim(t_elapsed, prop=0.25))

        iters = [length(exp["conv_data"].expData) for exp in exp_data]
        print("N_iters: mean: ", mean(iters), ", std: ", std(iters), "\n")

        print("Centralized - ")
        params, exp_data, t_elapsed = runExp(platform=:ACU, algorithm=:centralized, N_exps=N_exps, N_iter=1_000_000, N_agnts=N_agnts, N_acts=N_acts, convanalysis=true,
                                plotConvergenceRates=false, saveplots=false, params=ACU_params,
                                plotConvergenceRatesIndividual=true, errstats=false, stopping_value=1e-3, showplots=false, display=:final);
        t_elapsed_centralized[k, l, :] = t_elapsed
        t_elapsed_centralized_mean[k, l] = mean(trim(t_elapsed, prop=0.25))
        t_elapsed_centralized_std[k, l] = std(trim(t_elapsed, prop=0.25))

        iters = [length(exp["conv_data"].expData) for exp in exp_data]
        print("N_iters: mean: ", mean(iters), ", std: ", std(iters), "\n\n")
    end
end
##
file = matopen("/Users/martingurtner/Downloads/ACU_conv_trim_3.mat", "w")
write(file, "t_elapsed_admm", t_elapsed_admm)
write(file, "t_elapsed_admm_mean", t_elapsed_admm_mean)
write(file, "t_elapsed_admm_std", t_elapsed_admm_std)
write(file, "t_elapsed_centralized", t_elapsed_centralized)
write(file, "t_elapsed_centralized_mean", t_elapsed_centralized_mean)
write(file, "t_elapsed_centralized_std", t_elapsed_centralized_std)
write(file, "N_agnts_vals", N_agnts_vals)
write(file, "actuators_per_agent_vals", actuators_per_agent_vals)
close(file)

# print("++++++++++++++++++ DEP +++++++++++++++")
# DEP - distributed
N_agnts_vals = collect(40:20:160)
actuators_per_agent_vals = collect(16:6:40)
nagents2arrsize = (N_agnts, actuators_per_agent) -> (Int64(ceil(sqrt(actuators_per_agent*N_agnts))), Int64(ceil(sqrt(actuators_per_agent*N_agnts))))
N_exps = 10

t_elapsed_admm = zeros(length(N_agnts_vals), length(actuators_per_agent_vals), N_exps)
t_elapsed_admm_mean = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_admm_std = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_centralized = zeros(length(N_agnts_vals), length(actuators_per_agent_vals), N_exps)
t_elapsed_centralized_mean = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
t_elapsed_centralized_std = zeros(length(N_agnts_vals), length(actuators_per_agent_vals))
for (k, N_agnts) in enumerate(N_agnts_vals)
   for (l, actuators_per_agent) in enumerate(actuators_per_agent_vals)
       N_acts = nagents2arrsize(N_agnts, actuators_per_agent)

       print(k, "/", length(N_agnts_vals), " # Agents ", N_agnts, " | ", l, "/", length(actuators_per_agent_vals)," agent density values\n")

       print("ADMM - ")
       params, exp_data, t_elapsed = runExp(platform=:DEP, algorithm=:admm, N_exps=N_exps, N_iter=1_000_000, N_agnts=N_agnts, N_acts=N_acts, convanalysis=true,
                               plotConvergenceRates=false, saveplots=false, params=DEP_params,
                               plotConvergenceRatesIndividual=true, errstats=false, stopping_value=1e-3, showplots=false, display=:final);
       t_elapsed_admm[k, l, :] = t_elapsed
       t_elapsed_admm_mean[k, l] = mean(trim(t_elapsed, prop=0.25))
       t_elapsed_admm_std[k, l] = std(trim(t_elapsed, prop=0.25))

       iters = [length(exp["conv_data"].expData) for exp in exp_data]
       print("N_iters: mean: ", mean(iters), ", std: ", std(iters), "\n")

       print("Centralized - ")
       params, exp_data, t_elapsed = runExp(platform=:DEP, algorithm=:centralized, N_exps=N_exps, N_iter=1_000_000, N_agnts=N_agnts, N_acts=N_acts, convanalysis=true,
                               plotConvergenceRates=false, saveplots=false, params=DEP_params,
                               plotConvergenceRatesIndividual=true, errstats=false, stopping_value=1e-3, showplots=false, display=:final);
       t_elapsed_centralized_mean[k, l] = mean(trim(t_elapsed, prop=0.25))
       t_elapsed_centralized_std[k, l] = std(trim(t_elapsed, prop=0.25))

       iters = [length(exp["conv_data"].expData) for exp in exp_data]
       print("N_iters: mean: ", mean(iters), ", std: ", std(iters), "\n\n")
   end
   file = matopen("DEP_conv_trim.mat", "w")
   write(file, "t_elapsed_admm", t_elapsed_admm)
   write(file, "t_elapsed_admm_mean", t_elapsed_admm_mean)
   write(file, "t_elapsed_admm_std", t_elapsed_admm_std)
   write(file, "t_elapsed_centralized", t_elapsed_centralized)
   write(file, "t_elapsed_centralized_mean", t_elapsed_centralized_mean)
   write(file, "t_elapsed_centralized_std", t_elapsed_centralized_std)
   close(file)
end
#


text:
# docs/make.jl
# Make Documentation for AbstractWallets.jl

push!(LOAD_PATH, "../src/")

using Documenter, AbstractWallets

makedocs(
    modules   = [AbstractWallets],
    doctest   = true,
    clean     = false,
    linkcheck = false,
    format    = Documenter.HTML(prettyurls=!("local" in ARGS)),
    sitename  = "AbstractWallets.jl",
    authors   = "Brandon H Gomes",
    pages     = [
        "Home" => "index.md",
        "Implementing Wallet API" => [
            "Guide" => "impl/guide.md",
            "Example Projects" => "impl/examples.md",
        ],
        "Library" => [
            "Wallets and Portfolios" => "lib/wallets.md",
            "Transactions" => "lib/transactions.md",
            "Custom Wallet Lists" => "lib/lists.md",
        ]
    ],
)

deploydocs(
    repo = "github.com/bhgomes/AbstractWallets.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
using Random
using Distributions
using DSP: filt

export RedGaussian, PinkGaussian

# TODO: doesn't handle rand(::RedGaussian, m) for m samples of n-vectors yet (see _rand!)

"""
$(TYPEDEF)
Red Gaussian noise distribution for random variate generation.

# Example:
```julia-repl
julia> rand(RedGaussian(1000))
1000-element Array{Float64,1}:
   ⋮

julia> rand(RedGaussian(1000, 2.0))
1000-element Array{Float64,1}:
   ⋮
```
"""
Base.@kwdef struct RedGaussian{T} <: Distributions.ContinuousMultivariateDistribution
  n::Int
  σ::T = 1.0
end

RedGaussian(n) = RedGaussian(n=n)

"""
$(TYPEDEF)
Pink Gaussian noise distribution for random variate generation.

# Example:
```julia-repl
julia> rand(PinkGaussian(1000))
1000-element Array{Float64,1}:
   ⋮

julia> rand(PinkGaussian(1000, 2.0))
1000-element Array{Float64,1}:
   ⋮
```
"""
Base.@kwdef struct PinkGaussian{T} <: Distributions.ContinuousMultivariateDistribution
  n::Int
  σ::T = 1.0
end

PinkGaussian(n) = PinkGaussian(n=n)

function Random.rand!(rng::AbstractRNG, d::RedGaussian{T}, x::AbstractVector{T}) where {T<:Real}
  length(x) >= d.n || throw(ArgumentError("length of x must be at least n"))
  extra = 100
  v = rand(rng, Normal{T}(zero(T), d.σ), length(x)+extra)
  for j = 2:length(v)
    v[j] += v[j-1]
  end
  removedc!(v)
  x .= v[extra+1:end]/3.20377
  return x
end

function Random.rand!(rng::AbstractRNG, d::PinkGaussian{T}, x::AbstractVector{T}) where {T<:Real}
  length(x) >= d.n || throw(ArgumentError("length of x must be at least n"))
  # based on https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html
  hb = [0.049922035, -0.095993537, 0.050612699, -0.004408786]
  ha = [1.0, -2.494956002, 2.017265875, -0.522189400]
  extra = 1430
  v = rand(rng, Normal{T}(zero(T), d.σ), length(x)+extra)
  v = filt(hb, ha, v)/0.08680587859687908
  x .= v[extra+1:end]
  return x
end

Base.length(d::RedGaussian) = d.n
Base.length(d::PinkGaussian) = d.n

Base.rand(rng::AbstractRNG, d::RedGaussian) = Random.rand!(rng, d, zeros(typeof(d.σ), d.n))
Base.rand(rng::AbstractRNG, d::PinkGaussian) = Random.rand!(rng, d, zeros(typeof(d.σ), d.n))


text:
"""Shift integer to push the mantissa in the right position. Used to determine
round up or down in the tie case. `keepbits` is the number of mantissa bits to
be kept (i.e. not zero-ed) after rounding."""
function get_shift(::Type{T},keepbits::Integer) where {T<:Base.IEEEFloat}
    return Base.significand_bits(T) - keepbits
end

"""Returns for a Float-type `T` and `keepbits`, the number of mantissa bits to be
kept/non-zeroed after rounding, half of the unit in the last place as unsigned integer.
Used in round (nearest) to add ulp/2 just before round down to achieve round nearest.
Technically ulp/2 here is just smaller than ulp/2 which rounds down the ties. For
a tie round up +1 is added in `round(T,keepbits)`."""
function get_ulp_half(  ::Type{T},
                        keepbits::Integer
                        ) where {T<:Base.IEEEFloat}
    # convert to signed for arithmetic bitshift
    # trick to push in 0s from the left and 1s from the right
    return ~unsigned(signed(~Base.significand_mask(T)) >> (keepbits+1))
end

"""Returns a mask that's 1 for all bits that are kept after rounding and 0 for the
discarded trailing bits. E.g.
```
julia> get_keep_mask(Float16,5)
0xffe0
```."""
function get_keep_mask( ::Type{T},
                        keepbits::Integer
                        ) where {T<:Base.IEEEFloat}
    # convert to signed for arithmetic bitshift
    # trick to push in 1s from the left and 0s from the right
    return unsigned(signed(~Base.significand_mask(T)) >> keepbits)
end

"""Returns a mask that's `1` for a given `mantissabit` and `0` else. Mantissa bits
are positive for the mantissa (`mantissabit = 1` is the first mantissa bit), `mantissa = 0`
is the last exponent bit, and negative for the other exponent bits."""
function get_bit_mask(  ::Type{T},
                        mantissabit::Integer
                        ) where {T<:Base.IEEEFloat}
    # push the sign mask (1000....) in the right position
    return Base.sign_mask(T) >> (Base.exponent_bits(T) + mantissabit)
end

"""IEEE's round to nearest tie to even for Float16/32/64."""
function Base.round(x::T,               # Float to be rounded
                    ulp_half::UIntT,    # obtain from get_ulp_half,
                    shift::Integer,     # get_shift, and
                    keepmask::UIntT     # get_keep_mask
                    ) where {T<:Base.IEEEFloat,UIntT<:Unsigned}
    ui = reinterpret(UIntT,x)                       # bitpattern as uint
    ui += ulp_half + ((ui >> shift) & one(UIntT))   # add ulp/2 with tie to even
    return reinterpret(T,ui & keepmask)             # set trailing bits to zero
end

"""Scalar version of `round(::Float,keepbits)` that first obtains
`shift, ulp_half, keep_mask` and then rounds."""
function Base.round(x::T,
                    keepbits::Integer
                    ) where {T<:Base.IEEEFloat}
    return round(x,get_ulp_half(T,keepbits),get_shift(T,keepbits),get_keep_mask(T,keepbits))
end

"""IEEE's round to nearest tie to even for a float array `X` in-place. Calculates from `keepbits`
only once the variables `ulp_half`, `shift` and `keep_mask` and loops over every element of the
array."""
function round!(X::AbstractArray{T},            # any array with element type T
                keepbits::Integer               # how many mantissa bits to keep
                ) where {T<:Base.IEEEFloat}     # constrain element type to Float16/32/64

    ulp_half = get_ulp_half(T,keepbits)         # half of unit in last place (ulp)
    shift = get_shift(T,keepbits)               # integer used for bitshift
    keep_mask = get_keep_mask(T,keepbits)       # mask to zero trailing mantissa bits

    @inbounds for i in eachindex(X)             # apply rounding to each element
        X[i] = round(X[i],ulp_half,shift,keep_mask)
    end

    return X
end

"""IEEE's round to nearest tie to even for a float array `X` which returns a rounded copy of `X`."""
function Base.round(X::AbstractArray{T},        # any array with element type T
                    keepbits::Integer           # mantissa bits to keep
                    ) where {T<:Base.IEEEFloat} # constrain element type to Float32/64

    Xcopy = copy(X)                             # copy array to avoid in-place changes
    round!(Xcopy,keepbits)                      # in-place round the copied array
    return Xcopy
end

"""Checks a given `mantissabit` of `x` for eveness. 1=odd, 0=even. Mantissa bits
are positive for the mantissa (`mantissabit = 1` is the first mantissa bit), `mantissa = 0`
is the last exponent bit, and negative for the other exponent bits."""
function Base.iseven(x::T,
                    mantissabit::Integer
                    ) where {T<:Base.IEEEFloat}

    bitmask = get_bit_mask(T,mantissabit)
    return 0x0 == reinterpret(typeof(bitmask),x) & bitmask
end

"""Checks a given `mantissabit` of `x` for oddness. 1=odd, 0=even. Mantissa bits
are positive for the mantissa (`mantissabit = 1` is the first mantissa bit), `mantissa = 0`
is the last exponent bit, and negative for the other exponent bits."""
Base.isodd(x::T,mantissabit::Integer) where {T<:Base.IEEEFloat} = ~iseven(x,mantissabit)


text:
"""
Creates RF Pulses based on JSON files.
see testing/ for example files
"""
module RFPulse

# requirements
using JSON;
using BlochSim;
using PyPlot; pygui(true);

import Base: getindex, display;

export  AbstractRFPulse,
	HyperbolicSecant,
	ExternalWaveform,
	display,
	read_pulse,
	write_pulse,
	getindex,
	simulate!,
	plot_pulse,
	plot_simulation;

#
# TODOs:
#	Other Pulses
#	VERSE
#
#

"""
 Abstract RF Pulse:
 Dictionary of RF Properties
 Interface:
  "RFShape" => Complex RF pulse samples
  "GradShape" => Floating Point Gradient Samples
  "Time" => Vector of time points
"""
abstract type AbstractRFPulse end;

type HyperbolicSecant <: AbstractRFPulse
	properties::Dict{String,Any};

"""
  Create the HS Pulse
"""
  function HyperbolicSecant(config::Dict{String,Any})

    sysconfig   = config["SystemParameters"];
    pulseconfig = config["PulseParameters"];
    rfshape     = config["RFShape"];

    # system pars
    Gamma    = sysconfig["Gamma"];
    sys_max_grad = sysconfig["MaxGrad"];

    #RF pars
    ns    = pulseconfig["Samples"];
    dur   = pulseconfig["Duration"];
    slice = pulseconfig["SliceThickness"];
    Beta  = pulseconfig["HyperbolicSecantParameters"]["Beta"];
    Mu    = pulseconfig["HyperbolicSecantParameters"]["Mu"];
    #ref_b1= pulseconfig["RefB1"];
    b1_factor = pulseconfig["HyperbolicSecantParameters"]["B1ScaleFactor"];

    #adiabacity
    b1_hs_limit = sqrt(Mu) * Beta / (Gamma*2*pi);  #B1 must be above this limit for adiabacity
    ref_b1= b1_hs_limit * b1_factor;

    # create the time vector
    t = linspace(-dur/2.,dur/2,ns);

    b1_hs = ref_b1.*sech.(Beta.*t);         # Amplitude
    df_hs = -1.0.*Mu.*Beta.*tanh.(Beta.*t); # Frequency

    # bandwidth
    bw = Mu*Beta/pi;

    # Gradient Strength
    g_hs = bw/(Gamma*slice);  #mT/m

    # gradient vector
    refgrad  = g_hs;
    grad     = ones(ns); # const. gradient for HS pulses

    # rf vector
    rf_amp   = b1_hs./ref_b1;
    rf_phase = (df_hs .* t);
    rf       = cat(2,rf_amp,rf_phase);

    # set some properties
    config["PulseParameters"]["Bandwidth"] = bw;
    config["PulseParameters"]["RefGrad"] = refgrad;
    config["PulseParameters"]["MinSlice"] = 1.0e-3; # default
    config["PulseParameters"]["MaxSlice"] = bw/(Gamma*sys_max_grad);
    config["PulseParameters"]["RefB1"] = ref_b1;

    # add to dictionary
    props = config;
    props["RFShape"] = rf;
    props["GradShape"] = grad;
    props["Time"] = Array(linspace(0.0,1.0,ns));

    #finally calculate flip angle
    calculate_flip_angle!(props);

    new(props);
  end

end # HyperbolicSecant

type ExternalWaveform <: AbstractRFPulse
  properties::Dict{String,Any};

  """
  create it
  """
  function ExternalWaveform(config::Dict{String,Any});
    new(config);
  end
end

#getindex for all abstract RF pulses
getindex(p::AbstractRFPulse,idx) = p.properties[idx];

"""
 Reads a JSON config file and returns the requested RF Pulse Type
"""
function read_pulse(filename::String)
  fid = open(filename,"r");
  rawdata = read(fid);
  strdata = String(rawdata);
  config = JSON.parse(strdata);

  # create the RF pulse
  if config["PulseType"] == "HyperbolicSecant"
    pulse = HyperbolicSecant(config);
  else
    error("UnKnown Pulse type");
  end
  return pulse;
end

"""
  Writes the RFPulse type to a JSON file
"""
function write_pulse(pulse::AbstractRFPulse,filename::String)

	fid = open(filename,"w");

	JSON.print(fid,pulse.properties);

	close(fid);

end

"""
RF Pulse Interface to BlochSim.jl
performs an inplace modification to the RF pulse
"""
function simulate!(pulse::AbstractRFPulse)
  # parse the config dict for the simulation parameters

  ref_b1   = pulse["PulseParameters"]["RefB1"];
  ref_grad = pulse["PulseParameters"]["RefGrad"];
  dur      = pulse["PulseParameters"]["Duration"];

  Gamma    = pulse["SystemParameters"]["Gamma"];
  sim_pars = pulse["SimulationParameters"];
  ns       = sim_pars["PositionSamples"];
  nf       = sim_pars["FrequencySamples"];
  t1       = sim_pars["T1"];
  t2       = sim_pars["T2"];
  frange   = sim_pars["FrequencyRange"];
  prange   = sim_pars["PositionRange"];

  rf_magn  = pulse["RFShape"][:,1] .* ref_b1;
  rf_phase = pulse["RFShape"][:,2];
  rf       = rf_magn .* exp.(im .* rf_phase);
  grad     = pulse["GradShape"] .* ref_grad;
  t        = pulse["Time"] .* dur;

  pos_vec  = linspace(-prange/2.,prange/2.,ns);

  #TODO, fix frequency range and think about how to represent a 2D sim in the JSON file
  freq_vec = 0.0; #linspace(-frange/2.,frange/2.,nf);

  (m,m2) = BlochSim.sliceprofile(rf,grad,t,t1,t2,pos_vec,freq_vec,Gamma);

  output = Dict{String,Any}("Mx" => m[1,:],
			    "My" => m[2,:],
			    "Mz" => m[3,:],
			    "X"  => pos_vec,
			    "F"  => freq_vec);
  pulse["SimulationParameters"]["Output"] = output;
end

"""
Plots the RF Pulse
"""
function plot_pulse(pulse::AbstractRFPulse)
  ref_b1   = pulse["PulseParameters"]["RefB1"] * 1000.0; # convert to uT
  ref_grad = pulse["PulseParameters"]["RefGrad"];
  dur      = pulse["PulseParameters"]["Duration"] *1000.0; # convert to ms for viewing

  # get pulse shape from the dictionary
  t    = pulse["Time"];
  rf_magn    = pulse["RFShape"][:,1] .* ref_b1;
  rf_phase   = pulse["RFShape"][:,2];
  rf = rf_magn .* exp.(im .* rf_phase);
  grad = pulse["GradShape"] .* ref_grad;

  subplot(1,2,1);
  plot(t,abs.(rf));
  plot(t,real.(rf));
  plot(t,imag.(rf));
  title("RF Pulse");
  ylabel("B1 [uT]")
  xlabel("Time [ms]");
  subplot(1,2,2);
  plot(t,grad);
  title("gradient");
  ylabel("Grad [mT]")
  xlabel("Time [ms]");
end

"""
Plots the simulation results
"""
function plot_simulation(pulse::AbstractRFPulse)
  # get results from the dictionary
  sim = pulse["SimulationParameters"]["Output"];
  st  = pulse["PulseParameters"]["SliceThickness"];
  x   = sim["X"];
  mz = sim["Mz"];
  my = sim["My"];
  mx = sim["Mx"];
  mt = abs.(mx + im.*my);

  # ideal slice profile
  ideal_sp = ones(size(x));
  ideal_sp[ find((f)->abs.(f)<st/2.0,x)] = -1.0;

  subplot(1,2,1);
  plot(x,mz);
  plot(x,ideal_sp);
  title("Slice Profile");
  subplot(1,2,2);
  plot(x,mt);
  title("Signal");
end

"""
output pulse parameters to something more human readable
"""
function display(pulse::AbstractRFPulse)
  @printf("RF Pulse: %s \n",pulse["PulseType"]);
  @printf("\nSystem Parameters: \n");
  syspars = pulse["SystemParameters"];
  for (k,v) in syspars
    println(k,": ", v);
  end
  @printf("\nPulseParameters:\n");
  pulsepars = pulse["PulseParameters"];
  for (k,v) in pulsepars
    println(k,": ", v);
  end
end

"""
 Calculates the Ffip angle by integrating the RF Pulse envelope
 args can be RFpulse of config dictionary
"""
function calculate_flip_angle!(config::Union{AbstractRFPulse,Dict{String,Any}})
  Gamma  = config["SystemParameters"]["Gamma"];
  rf_amp = config["RFShape"][:,1] .* config["PulseParameters"]["RefB1"];
  t      = config["Time"] .* config["PulseParameters"]["Duration"];

  fa = 2pi* Gamma * sum(rf_amp .* t);
  config["PulseParameters"]["FlipAngleDeg"] = rad2deg(fa);
end

end #module


text:
@inline function User!(kwargs_dict::Dict)
    if haskey(kwargs_dict, :user)
        user = User(kwargs_dict[:user])
        delete!(kwargs_dict, :user)
        return user
    end
    throw(UserNotProvidedError("Username of registering user not provided"))
end


text:
using Test, NonconvexCore
import JuMP

@testset "Symbol variable names" begin
    model = JuMP.Model()
    JuMP.@variable(model, x >= 0)
    JuMP.@variable(model, 0 <= y <= 3)
    JuMP.@objective(model, Min, 20y)
    JuMP.@constraint(model, c1, 6x + 8y >= 100)
    JuMP.@constraint(model, c2, 7x + 12y >= 120)
    JuMP.@constraint(model, c3, x + y == 25)
    JuMP.@constraint(model, c4, x + y <= 50)

    dict_model = DictModel(model)
    obj = NonconvexCore.getobjective(dict_model)
    ineq = NonconvexCore.getineqconstraints(dict_model)
    eq = NonconvexCore.geteqconstraints(dict_model)
    x = NonconvexCore.getinit(dict_model)

    @test getmin(dict_model) == OrderedDict(:x => 0.0, :y => 0.0)
    @test getmax(dict_model) == OrderedDict(:x => Inf, :y => 3.0)
    @test obj(x) == 20 * x[:y]
    @test ineq(x) == [
        100 - 6 * x[:x] - 8 * x[:y],
        120 - 7 * x[:x] - 12 * x[:y],
        x[:x] + x[:y] - 50,
    ]
    @test eq(x) == [x[:x] + x[:y] - 25]

    vec_model, _, _ = NonconvexCore.tovecmodel(dict_model)
    vec_x = NonconvexCore.getinit(vec_model)
    vec_obj = NonconvexCore.getobjective(vec_model)
    vec_ineq = NonconvexCore.getineqconstraints(vec_model)
    vec_eq = NonconvexCore.geteqconstraints(vec_model)

    val0, grad0 = NonconvexCore.value_gradient(vec_obj, vec_x)
    val1, jac1 = NonconvexCore.value_jacobian(vec_ineq, vec_x)
    val2, jac2 = NonconvexCore.value_jacobian(vec_eq, vec_x)

    @test val0 == obj(x)
    @test val1 == ineq(x)
    @test val2 == eq(x)

    @test grad0 == [0, 20]
    @test jac1 == [-6 -8; -7 -12; 1 1;]
    @test jac2 == [1 1;]
end

@testset "x[ind] variable names 1" begin
    model = JuMP.Model()
    lbs = [0.0, 0.0]
    ubs = [Inf, 3.0]
    JuMP.@variable(model, x[i=1:2], lower_bound = lbs[i], upper_bound = ubs[i])
    JuMP.@objective(model, Min, 12x[1])
    JuMP.@constraint(model, c1, 6x[1] + 8x[2] >= 100)
    JuMP.@constraint(model, c2, 7x[1] + 12x[2] >= 120)
    JuMP.@constraint(model, c3, x[1] + x[2] == 25)
    JuMP.@constraint(model, c4, x[1] + x[2] <= 50)

    dict_model = DictModel(model)
    obj = NonconvexCore.getobjective(dict_model)
    ineq = NonconvexCore.getineqconstraints(dict_model)
    eq = NonconvexCore.geteqconstraints(dict_model)
    x = NonconvexCore.getinit(dict_model)

    @test getmin(dict_model) == OrderedDict(:x => [0.0, 0.0])
    @test getmax(dict_model) == OrderedDict(:x => [Inf, 3.0])
    @test obj(x) == 12 * x[:x][1]
    @test ineq(x) == [
        100 - 6 * x[:x][1] - 8 * x[:x][2],
        120 - 7 * x[:x][1] - 12 * x[:x][2],
        x[:x][1] + x[:x][2] - 50,
    ]
    @test eq(x) == [x[:x][1] + x[:x][2] - 25]

    vec_model, _, _ = NonconvexCore.tovecmodel(dict_model)
    vec_x = NonconvexCore.getinit(vec_model)
    vec_obj = NonconvexCore.getobjective(vec_model)
    vec_ineq = NonconvexCore.getineqconstraints(vec_model)
    vec_eq = NonconvexCore.geteqconstraints(vec_model)

    val0, grad0 = NonconvexCore.value_gradient(vec_obj, vec_x)
    val1, jac1 = NonconvexCore.value_jacobian(vec_ineq, vec_x)
    val2, jac2 = NonconvexCore.value_jacobian(vec_eq, vec_x)

    @test val0 == obj(x)
    @test val1 == ineq(x)
    @test val2 == eq(x)

    @test grad0 == [12, 0]
    @test jac1 == [-6 -8; -7 -12; 1 1;]
    @test jac2 == [1 1;]
end

@testset "x[ind] variable names 2" begin
    model = JuMP.Model()
    lbs = [0.0, 0.0]
    ubs = [Inf, 3.0]
    JuMP.@variable(model, x[i=2:-1:1], lower_bound = lbs[i], upper_bound = ubs[i])
    JuMP.@objective(model, Min, 12x[1])
    JuMP.@constraint(model, c1, 6x[1] + 8x[2] >= 100)
    JuMP.@constraint(model, c2, 7x[1] + 12x[2] >= 120)
    JuMP.@constraint(model, c3, x[1] + x[2] == 25)
    JuMP.@constraint(model, c4, x[1] + x[2] <= 50)

    dict_model = DictModel(model)
    obj = NonconvexCore.getobjective(dict_model)
    ineq = NonconvexCore.getineqconstraints(dict_model)
    eq = NonconvexCore.geteqconstraints(dict_model)
    x = NonconvexCore.getinit(dict_model)

    @test getmin(dict_model) == OrderedDict(:x => JuMP.Containers.DenseAxisArray(lbs[end:-1:1], 2:-1:1))
    @test getmax(dict_model) == OrderedDict(:x => JuMP.Containers.DenseAxisArray(ubs[end:-1:1], 2:-1:1))
    @test obj(x) == 12 * x[:x][1]
    @test ineq(x) == [
        100 - 6 * x[:x][1] - 8 * x[:x][2],
        120 - 7 * x[:x][1] - 12 * x[:x][2],
        x[:x][1] + x[:x][2] - 50,
    ]
    @test eq(x) == [x[:x][1] + x[:x][2] - 25]

    vec_model, _, _ = NonconvexCore.tovecmodel(dict_model)
    vec_x = NonconvexCore.getinit(vec_model)
    vec_obj = NonconvexCore.getobjective(vec_model)
    vec_ineq = NonconvexCore.getineqconstraints(vec_model)
    vec_eq = NonconvexCore.geteqconstraints(vec_model)

    val0, grad0 = NonconvexCore.value_gradient(vec_obj, vec_x)
    val1, jac1 = NonconvexCore.value_jacobian(vec_ineq, vec_x)
    val2, jac2 = NonconvexCore.value_jacobian(vec_eq, vec_x)

    @test val0 == obj(x)
    @test val1 == ineq(x)
    @test val2 == eq(x)

    @test grad0 == [0, 12]
    @test jac1 == [-8 -6; -12 -7; 1 1;]
    @test jac2 == [1 1;]
end

@testset "x[(ind1, ind2)] variable names" begin
    model = JuMP.Model()
    lbs = [0.0, 0.0]
    ubs = [Inf, 3.0]
    inds = [(1,1), (2,1)]
    JuMP.@variable(model, x[ind = inds], lower_bound = lbs[ind[1]], upper_bound = ubs[ind[1]])
    JuMP.@objective(model, Min, 12x[(1,1)] + 20x[(2,1)])
    JuMP.@constraint(model, c1, 6x[(1,1)] + 8x[(2,1)] >= 100)
    JuMP.@constraint(model, c2, 7x[(1,1)] + 12x[(2,1)] >= 120)
    JuMP.@constraint(model, c3, x[(1,1)] + x[(2,1)] == 25)
    JuMP.@constraint(model, c4, x[(1,1)] + x[(2,1)] <= 50)

    dict_model = DictModel(model)
    obj = NonconvexCore.getobjective(dict_model)
    ineq = NonconvexCore.getineqconstraints(dict_model)
    eq = NonconvexCore.geteqconstraints(dict_model)
    x = NonconvexCore.getinit(dict_model)

    @test getmin(dict_model) == OrderedDict(:x => JuMP.Containers.DenseAxisArray(lbs, [(1, 1), (2, 1)]))
    @test getmax(dict_model) == OrderedDict(:x => JuMP.Containers.DenseAxisArray(ubs, [(1, 1), (2, 1)]))
    @test obj(x) == 12 * x[:x][(1, 1)] + 20 * x[:x][(2, 1)]
    @test ineq(x) == [
        100 - 6 * x[:x][(1, 1)] - 8 * x[:x][(2, 1)],
        120 - 7 * x[:x][(1, 1)] - 12 * x[:x][(2, 1)],
        x[:x][(1, 1)] + x[:x][(2, 1)] - 50,
    ]
    @test eq(x) == [x[:x][(1, 1)] + x[:x][(2, 1)] - 25]

    vec_model, _, _ = NonconvexCore.tovecmodel(dict_model)
    vec_x = NonconvexCore.getinit(vec_model)
    vec_obj = NonconvexCore.getobjective(vec_model)
    vec_ineq = NonconvexCore.getineqconstraints(vec_model)
    vec_eq = NonconvexCore.geteqconstraints(vec_model)

    val0, grad0 = NonconvexCore.value_gradient(vec_obj, vec_x)
    val1, jac1 = NonconvexCore.value_jacobian(vec_ineq, vec_x)
    val2, jac2 = NonconvexCore.value_jacobian(vec_eq, vec_x)

    @test val0 == obj(x)
    @test val1 == ineq(x)
    @test val2 == eq(x)

    @test grad0 == [12, 20]
    @test jac1 == [-6 -8; -7 -12; 1 1;]
    @test jac2 == [1 1;]
end


text:
using TerminalLoggers
using REPL
using REPL:Terminals
using Logging: global_logger

global_logger(TerminalLogger())

const log_term = Terminals.TTYTerminal("", stdin, stdout, stderr)

function replprint(output::String;
    newline::Int=0, clearline::Int=1, color::Symbol=:white, bold::Bool=false,
    prefix::String="", prefix_color::Symbol=:green, prefix_bold::Bool=true, indent_level::Int=0)

    # clear line
    if clearline > 0
        for i in 1:(clearline + 1)
            Terminals.clear_line(log_term)
        end
    end

    # indent
    print(" "^indent_level)

    # print info
    if !isempty(prefix)
        printstyled(prefix, color=prefix_color, bold=prefix_bold)
    end
    printstyled(output, color=color, bold=bold)

    # print new lines
    [println() for _ ∈ 1:newline]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 41d2dbc6-5d31-11eb-12d5-d9c5db4e22c0
begin
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots
end

# ╔═╡ 11ec86dc-5d31-11eb-2d63-ddafcd73800d
md"""
# Chapter12 Random Walk
"""

# ╔═╡ 53376968-5d31-11eb-1c73-d9d50607cab9
N = 21

# ╔═╡ 539f2eea-5d31-11eb-1969-25432dee8823
true_values = -1:0.1:1

# ╔═╡ 56dfb4a8-5d31-11eb-217e-d75aa440a47d
Base.@kwdef struct RecordRMS <: AbstractHook
    rms::Vector{Float64}=[]
end

# ╔═╡ 5c1854ca-5d31-11eb-013b-a1aea66007e1
(h::RecordRMS)(::PostEpisodeStage, agent, env) = push!(h.rms, sqrt(mean((agent.policy.learner.approximator.table[2:end-1] - true_values[2:end-1]).^2)))

# ╔═╡ 82f0ec74-5d31-11eb-0bbc-979653a27bf4
function create_agent_env(α, λ)
    env = RandomWalk1D(N=21)
    ns, na =  length(state_space(env)),  length(action_space(env))
    agent = Agent(
        policy=VBasedPolicy(
            learner=TDλReturnLearner(
                approximator=TabularVApproximator(;n_state=ns, opt=Descent(α)),
                γ=1.0,
                λ=λ
            ),
            mapping = (env, V) -> rand(1:na)
        ),
        trajectory=VectorSARTTrajectory()
    )
    agent, env
end

# ╔═╡ 863aba54-5d31-11eb-21d1-ab6ec717c0c4
function records(α, λ, nruns=10)
    rms = []
    for _ in 1:nruns
        hook = RecordRMS()
        run(create_agent_env(α, λ)..., StopAfterEpisode(10, is_show_progress=false),hook)
        push!(rms, mean(hook.rms))
    end
    mean(rms)
end

# ╔═╡ 8a1b50b4-5d31-11eb-26bd-03f9dbaf549d
begin
	As = [0:0.1:1, 0:0.1:1, 0:0.1:1, 0:0.1:1, 0:0.1:1, 0:0.05:0.5, 0:0.02:0.2, 0:0.01:0.1]
	Λ = [0., 0.4, .8, 0.9, 0.95, 0.975, 0.99, 1.]
	p = plot(legend=:topright)
	for (A, λ) in zip(As, Λ)
		plot!(p, A, [records(α, λ) for α in A], label="lambda = $λ")
	end
	ylims!(p, (0.25, 0.55))
	p
end

# ╔═╡ Cell order:
# ╟─11ec86dc-5d31-11eb-2d63-ddafcd73800d
# ╠═41d2dbc6-5d31-11eb-12d5-d9c5db4e22c0
# ╠═53376968-5d31-11eb-1c73-d9d50607cab9
# ╠═539f2eea-5d31-11eb-1969-25432dee8823
# ╠═56dfb4a8-5d31-11eb-217e-d75aa440a47d
# ╠═5c1854ca-5d31-11eb-013b-a1aea66007e1
# ╠═82f0ec74-5d31-11eb-0bbc-979653a27bf4
# ╠═863aba54-5d31-11eb-21d1-ab6ec717c0c4
# ╠═8a1b50b4-5d31-11eb-26bd-03f9dbaf549d


text:
using LinearAlgebra
using Plots
using Printf

include("helpers.jl")
include("simple.jl")

# Create an abstract type for first-order descent methods
abstract type FirstOrderMethod end

#*******************************************************************************
# Adam Implementation
mutable struct Adam <: FirstOrderMethod
    α; γv; γs; ϵ; k; v; s
end

adam_init!(A::Adam) = A

function adam_step!(A::Adam, ∇f, x)
    α, γv, γs, ϵ, k, s, v, g = A.α , A.γv, A.γs, A.ϵ, A.k, A.s, A.v, ∇f(x)
    v[:] = γv*v + (1-γv)*g
    s[:] = γs*s + (1-γs)*g.*g
    A.k = k += 1
    v̂ = v ./ (1 - γv^k)
    ŝ = s ./ (1 - γs^k)
    return x - α*v̂ ./ (sqrt.(ŝ) .+ ϵ)
end

function optimize_and_store_path(M, f, g, x_best, n)
    x1_arr = Float64[]; x2_arr = Float64[];
    push!(x1_arr, x_best[1]); push!(x2_arr, x_best[2]);
    k = 1
    while k < n
        x_best = adam_step!(M,g,x_best)
        push!(x1_arr, x_best[1]);push!(x2_arr, x_best[2])
        k = k + 1
    end
    return x_best, x1_arr, x2_arr
end

function optimize_and_store_y_vals(M, f, g, x_best, n)
    y_arr = Float64[];
    push!(y_arr, f(x_best))
    k = 1
    while k < n
        x_best = adam_step!(M,g,x_best)
        push!(y_arr, f(x_best))
        k = k + 1
    end
    return y_arr
end

function plot_rosenbrock_path(f, g, n)
    x01 = [3.0, -3.0]; x02 = [0.0, 3.0]; x03 = [-1.7, 2.8]
    M = Adam(0.5, 0.4, 0.999, 1e-8, 0.0, [0.0,0.0], [0.0,0.0])
    x1_best, x11_arr, x21_arr = optimize_and_store_path(M, f, g, x01, n);
    M = Adam(0.5, 0.4, 0.999, 1e-8, 0.0, [0.0,0.0], [0.0,0.0])
    x2_best, x12_arr, x22_arr = optimize_and_store_path(M, f, g, x02, n);
    M = Adam(0.5, 0.4, 0.999, 1e-8, 0.0, [0.0,0.0], [0.0,0.0])
    x3_best, x13_arr, x23_arr = optimize_and_store_path(M, f, g, x03, n);

    # Plot Rosenbrock Plot with select contours
    x1 = -4.0:0.01:4.0
    x2 = -4.0:0.01:4.0
    plot_f(x1, x2) = begin
        (1.0 - x1)^2 + 100.0 * (x2 - x1^2)^2
    end
    X = repeat(reshape(x1, 1, :), length(x2), 1)
    Y = repeat(x2, 1, length(x1))
    Z = map(plot_f, X, Y)
    #h = [minimum(Z), 2.5, 5, 10, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000, 25000, maximum(Z)]
    h = [minimum(Z), 2.5, 10, 25, 100, 250, 500, 1000, 2500, 5000, 10000, 15000, 25000, 35000,maximum(Z)]
    #h = [minimum(Z), 2.5, 25, 250, 500, 1000, 2500, 5000,10000,25000,50000,maximum(Z)]
    contour(x1, x2, Z, levels=h, fill=false, c=cgrad(:viridis, rev = true), linewidth=2.0)

    # Plot first path
    plot!(x11_arr,x21_arr,color=:magenta, label = "x0: (3.0, -3.0)",legend = :bottomleft,
          title = "Adam Method to Minimize Rosenbrock's Function", titlefontsize = 12,
          xlabel = "x1", ylabel = "x2", dpi=300, colorbar = false)
    scatter!([x11_arr], [x21_arr], color=:magenta, label = "", markersize = 3, outline = :magenta)
    scatter!([x01[1]], [x01[2]], color=:magenta,  marker = :rect,label = "", markersize = 5)
    scatter!([x1_best[1]], [x1_best[2]], color=:magenta, marker = :star4, label = "", markersize = 7)

    # Plot second path
    plot!(x12_arr,x22_arr,color=:red, label = "x0: (0.0, 3.0)")
    scatter!([x12_arr], [x22_arr], color=:red, label = "", markersize = 3, outline = false)
    scatter!([x02[1]], [x02[2]], color=:red, marker = :rect, label = "", markersize = 5)
    scatter!([x2_best[1]], [x2_best[2]], color=:red, marker = :star4, label = "", markersize = 7)

    # Plot third path
    plot!(x13_arr,x23_arr,color=:darkorange, label = "x0: (-1.7, 2.8)")
    scatter!([x13_arr], [x23_arr], color=:darkorange, label = "", markersize = 3, outline = false)
    scatter!([x03[1]], [x03[2]], color=:darkorange,  marker = :rect, label = "", markersize = 5)
    scatter!([x3_best[1]], [x3_best[2]], color=:darkorange, marker = :star4, label = "", markersize = 7)

    scatter!([1], [1], color=:dodgerblue, marker = :star5, label = "", markersize = 10)

    png("adam.png")

end

function return_y_vals(f, g, x0, n, prob)
    x_best = x0; x_arr = Float64[]; y_arr = Float64[]
    push!(x_arr, x0[1]); push!(y_arr, x0[2])
    if prob == "simple1"
        M = Adam(0.6, 0.4, 0.999, 1e-8, 0.0, [0.0,0.0], [0.0,0.0])
        y_arr = optimize_and_store_y_vals(M, f, g, x0, n)
    elseif prob == "simple2"
        M = Adam(0.5, 0.4, 0.999, 1e-8, 0.0, [0.0,0.0], [0.0,0.0])
        y_arr = optimize_and_store_y_vals(M, f, g, x0, n)
    elseif prob == "simple3"
        M = Adam(0.4, 0.6, 0.999, 1e-8, 0.0, [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0])
        y_arr = optimize_and_store_y_vals(M, f, g, x0, n)
    else
        print("Input a simple problem")
        return
    end

    return y_arr
end

n = 100
my_x = plot_rosenbrock_path(rosenbrock, rosenbrock_gradient, n)

x0_r = rosenbrock_init()
x0_r_x1 = round(x0_r[1], digits=3); x0_r_x2 = round(x0_r[2], digits=3)
rosenbrock_convergence =
    return_y_vals(rosenbrock, rosenbrock_gradient, x0_r, n, "simple1");

x0_h = himmelblau_init()
x0_h_x1 = round(x0_h[1], digits=3); x0_h_x2 = round(x0_h[2], digits=3)
himmelblau_convergence =
    return_y_vals(himmelblau, himmelblau_gradient, x0_h, n, "simple2");

x0_p = powell_init()
x0_p_x1 = round(x0_p[1], digits=3); x0_p_x2 = round(x0_p[2], digits=3)
x0_p_x3 = round(x0_p[3], digits=3); x0_p_x4 = round(x0_p[4], digits=3)
powell_convergence =
    return_y_vals(powell, powell_gradient, x0_p, n, "simple3");

rosenbrock_convergence1 =
    return_y_vals(rosenbrock, rosenbrock_gradient, [3.0, -3.0], n, "simple1");
rosenbrock_convergence2 =
    return_y_vals(rosenbrock, rosenbrock_gradient, [0.0, 3.0], n, "simple1");
rosenbrock_convergence3 =
    return_y_vals(rosenbrock, rosenbrock_gradient, [-1.7, 2.8], n, "simple1");

himmelblau_convergence1 =
    return_y_vals(himmelblau, himmelblau_gradient, [3.0, -3.0], n, "simple2");
himmelblau_convergence2 =
    return_y_vals(himmelblau, himmelblau_gradient, [0.0, 3.0], n, "simple2");
himmelblau_convergence3 =
    return_y_vals(himmelblau, himmelblau_gradient, [-1.7, 2.8], n, "simple2");

powell_convergence1 =
    return_y_vals(powell, powell_gradient, [3.0, -3.0, 3.0, -3.0], n, "simple3");
powell_convergence2 =
    return_y_vals(powell, powell_gradient, [0.0, 3.0, 0.0, 3.0], n, "simple3");
powell_convergence3 =
    return_y_vals(powell, powell_gradient, [-1.7, 2.8, -1.7, 2.8], n, "simple3");

x_vals = 1:n

plot(x_vals, rosenbrock_convergence1, lw = 3, label = "x0: (3.0, -3.0)",legend = :topright,
    title="Convergence for Rosenbrock Function",xlabel="Iteration",ylabel="f(x)",dpi=300)
plot!(x_vals, rosenbrock_convergence2, lw = 3, label = "x0: (0.0, 3.0)")
plot!(x_vals, rosenbrock_convergence3, lw = 3, label = "x0: (-1.7, 2.8)")
png("rosenbrock_convergence.png")

plot(x_vals, himmelblau_convergence1, lw = 3, label = "x0: (3.0, -3.0)",legend = :topright,
    title="Convergence for Himmelblau Function",xlabel="Iteration",ylabel="f(x)",dpi=300)
plot!(x_vals, himmelblau_convergence2, lw = 3, label = "x0: (0.0, 3.0)")
plot!(x_vals, himmelblau_convergence3, lw = 3, label = "x0: (-1.7, 2.8)")
png("himmelblau_convergence.png")

plot(x_vals, powell_convergence1, lw = 3, label = "x0: (3.0, -3.0, 3.0, -3.0)",legend = :topright,
    title="Convergence for Powell Function",xlabel="Iteration",ylabel="f(x)",dpi=300)
plot!(x_vals, powell_convergence2, lw = 3, label = "x0: (0.0, 3.0, 0.0, 3.0)")
plot!(x_vals, powell_convergence3, lw = 3, label = "x0: (-1.7, 2.8, -1.7, 2.8)")
png("powell_convergence.png")

##

print("HELLO")


text:
# Everything about intersecting
export meet_check, meet, ∧, stab_count, is_simple

"""
`meet_check(L::HLine,LL::HLine)` determines if two lines intersect.
Also for any combination of lines, segments, or rays.
"""
function meet_check(L::HLine, LL::HLine)::Bool
    if L == LL
        return false
    end
    s = L.s
    t = L.t
    ss = LL.s
    tt = LL.t

    if s < ss < t < tt
        return true
    end
    if ss < s < tt < t
        return true
    end
    return false
end

function meet_check(L::HLine, S::HSegment)::Bool
    LL = HLine(S)  # extend S
    if !meet_check(L,LL)
        return false
    end

    P = meet(L,LL)
    return in(P,S)
end

meet_check(S::HSegment,L::HLine)::Bool = meet_check(L,S)

function meet_check(S::HSegment, SS::HSegment)::Bool
    L = HLine(S)
    LL = HLine(SS)

    if L==LL  # see if they share exactly one endpoint
        A,B = endpoints(S)
        AA,BB = endpoints(SS)
        if A==AA && between(B,A,BB)
            return true
        end
        if A==BB && between(A,B,BB)
            return true
        end
        if B==AA && between(A,B,BB)
            return true
        end
        if B==BB && between(A,B,AA)
            return true
        end
        return false
    end

    if !meet_check(L,LL)
        return false
    end

    P = meet(L,LL)
    if in(P,S) && in(P,SS)
        return true
    end
    return false
end

function meet_check(L::HLine, R::HRay)::Bool
    LL = HLine(R)
    if !meet_check(L,LL)
        return false
    end
    P = meet(L,LL)
    return in(P,R)
end
meet_check(R::HRay, L::HLine)::Bool = meet_check(L,R)

function meet_check(R::HRay, S::HSegment)::Bool
    LR = HLine(R)
    LS = HLine(S)

    if LR == LS   # check if endpoints match up
        A,B = endpoints(S)
        X = get_vertex(R)
        if A==X
            return !in(B,R)
        end
        if B==X
            return !in(A,R)
        end
        return false
    end
    # R and S noncollinear
    if !meet_check(LR,LS)
        return false
    end
    P = meet(LR,LS)
    return in(P,R) && in(P,S)
end
meet_check(S::HSegment, R::HRay)::Bool = meet_check(R,S)

function meet_check(R::HRay, RR::HRay)::Bool
    L = HLine(R)
    LL = HLine(RR)

    if L==LL   # make sure same vertex, opposite directions
        if get_vertex(R) != get_vertex(RR)
            return false
        end
        if abs(R.t - RR.t) <= THRESHOLD * eps(1.0)
            return false
        end
        return true
    end

    if !meet_check(L,LL)
        return false
    end
    P = meet(L,LL)
    return in(P,R) && in(P,RR)
end

################################################################

"""
`meet(L,LL)` finds a point on lines `L` and `LL` or throws an
error if they don't intersect.  Also for a line and a segment, or
two segments.

See `meet_check`.
"""
function meet(L::HLine, LL::HLine)::HPoint
    @assert meet_check(L,LL) "The lines do not intersect"

    s = L.s   # artificially rotate to 0
    t = L.t
    ss = LL.s
    tt = LL.t

    a = exp(im*(t-s))
    b = exp(im*(ss-s))
    c = exp(im*(tt-s))

    A = real(in_up(a))
    B = real(in_up(b))
    C = real(in_up(c))

    R = abs(B-C)/2
    Z = (B+C)/2

    y = sqrt(R^2 - (A-Z)^2)

    p = up_in(A + im*y)*exp(im*s)

    return HPoint(p)
end

function meet(L::HLine, S::HSegment)::HPoint
    @assert meet_check(L,S) "The line and segment do not intersect at a unique point"
    p = meet(L, HLine(S))
    return p
end

meet(S::HSegment,L::HLine)::HPoint = meet(L,S)

function meet(S::HSegment,SS::HSegment)
    @assert meet_check(S,SS) "The segments do not intersect at a unique point"

    L = HLine(S)
    LL = HLine(SS)

    if L != LL
        p = meet(HLine(S),HLine(SS))
        return p
    end

    # special case: segments overlap in an end point
    A,B = endpoints(S)
    AA,BB = endpoints(SS)

    if A==AA || A==B
        return A
    end
    if B==AA || B==BB
        return B
    end
    # Should never get here
    @error "Programming error in meet(HSegment, HSegment)"
end

function meet(L::HLine, R::HRay)::HPoint
    @assert meet_check(L,R) "The line and ray do not intersect at a unique point"
    return meet(L, HLine(R))
end
meet(R::HRay, L::HLine)::HPoint = meet(L,R)

function meet(R::HRay, S::HSegment)::HPoint
    @assert meet_check(R,S) "The ray and segment do not intersect at a unique point"
    LR = HLine(R)
    LS = HLine(S)
    if LR != LS
        return meet(HLine(R, HLine(S)))
    end
    # must have an end of S as the common point
    A,B = endpoints(S)
    if in(A,R)
        return A
    end
    return B
end
meet(S::HSegment, R::HRay)::HPoint = meet(R,S)

function meet(R::HRay, RR::HRay)::HPoint
    @assert meet_check(R,RR) "The two rays do not instersect at a unique point"
    if HLine(R) == HLine(RR)
        return get_vertex(R)
    end
    return meet(HLine(R,HLine(RR)))
end

(∧)(L::HLinear, LL::HLinear)::HPoint = meet(L,LL)

"""
`stab_count(R,targets)` is used to count how many items in `targets` are
intersected by `R` where `targets` is a list of `HLinear` objects and `R`
is an `HLinear`.
"""
function stab_count(R::S, targets::Array{T,1}) where {S<:HLinear, T <: HLinear}
    count = 0
    nt = length(targets)
    for j = 1:nt
        if meet_check(R,targets[j])
            count += 1
        end
    end
    return count
end

"""
`in(p::HPoint, X::HPolygon)` determines if the point lies on the boundary or,
or is interior to, the polygon. Be sure that `polygon_check(X)` and
`is_simple(X)` both return `true`.
"""
function in(p::HPoint, X::HPolygon)::Bool
    for a in X.plist
        if p==a
            return true
        end
    end

    ss = sides(X)
    for s in ss
        if in(p,s)
            return true
        end
    end

    # Get a ray that doesn't go through a vertex (and therefore not through a side)
    R = RandomHRay(p)
    while true
        for a in X.plist
            if in(a,R)
                R = RandomHRay()
                continue
            end
        end
        break
    end

    return stab_count(R,ss)%2 == 1
end