text:
function mpc = c118swf
%CASE118    Power flow data for IEEE 118 bus test case.
%   Please see CASEFORMAT for details on the case file format.
%   This data was converted from IEEE Common Data Format
%   (ieee118cdf.txt) on 20-Sep-2004 by cdf2matp, rev. 1.11
%   See end of file for warnings generated during conversion.
%
%   Converted from IEEE CDF file from:
%       https://labs.ece.uw.edu/pstca/
%   With baseKV data take from the PSAP format file from the same site,
%   added manually on 10-Mar-2006.
%
%  08/25/93 UW ARCHIVE           100.0  1961 W IEEE 118 Bus Test Case
%
% Updates: Sep 15'2010 by CM and CZ. Buses grouped on four areas;
%        Branch thermal limits incorporated based on: KV level, same
%        cable, and safe security factor (assuming bundled config on 345KV
%        and double circuit on 161 and 138KV lines).
%        Transformer capacity for 9 units also set.
%        Split data for parallel lines that initially had aggregated params
%        RAMP_AGC,RAMP_10 & RAMP_30 set to 20% of PMAX for each generator
%        function change of name to c118
% Aug 13'2012 by CM and DM.
%         The 118-bus test system was modified in order to allow the
%         loading factor to be considerably increased without compromising
%         the feasibility of the system. The procedure followed is
%         described here.
%
%         First, all the synchronous capacitors were removed from the
%         system in order to decrease the size of the model. Those were
%         initially 22 generators. Two of them were added again (generators
%         at buses 110 and 104) because there are currently hydro
%         generators close to the geographic locations where the
%         synchronous capacitors were located originally. This overall
%         retirement yielded a base case consisting of 34 generators, that
%         is the original 54 generators minus 20 synchronous generators
%         retired.
%
%         Second, the set of relevant contingencies was chosen from all
%         possible contingencies. The generator contingencies selected were
%         those of generators located at buses 10, 80, and 89. The branch
%         contingencies selected were those of branches located between
%         buses 17-30, 30-38, and 5-11. Those contingencies showed to have
%         important impact over the system when using a dispatch deviation
%         criterion. The criterion was applied using both DC and AC
%         versions of the OPF.
%
%         Third, the loading factor range  was obtained for the set of
%         contingencies for both the DC and AC versions of the model. For
%         the DC case, all the loading factors of the base case (no
%         contingencies) and all the contingencies were above 1.6. For the
%         AC case, however, the loading factor of the base case was around
%         1.3, and even worse, the loading factor of the worst contingency
%         was around 1.1. Since similar loading factor conditions were
%         desired for both DC and AC cases, some changes or upgrades were
%         required.
%
%         In the AC base case there were voltage problems in buses 70
%         and 76, and that was causing the system to become infeasible when
%         exceeding a loading factor of 1.3. Apparently, there was not
%         enough reactive power to be dispatched in that zone since all the
%         synchronous capacitors were removed. Hence, two synchronous
%         capacitors were added in that particular zone, namely generators
%         at buses 76 and 77. Additionally, the reactive power limits of the
%         generator at bus 77 were increased from [-20, 70] MVar to [-50,
%         330] MVar in order to obtain a loading factor above 1.6 for the
%         base case.
%
%         Furthermore, in the case under the contingency of the generator at
%         bus 10, again with the AC case, the loading factor was the one
%         close to 1.1. To improve that factor, the synchronous capacitor
%         at bus 6 was added and its reactive power limits increased from
%         [-13, 50] MVar to [-50, 250] MVar. With these changes, the
%         voltage level deficiencies at bus 6 were solved, at least for
%         loading factors up to 1.5. Beyond this loading factor it gets
%         more difficult to improve it without resorting to the addition of
%         more synchronous generators. Therefore, the maximum loading
%         factor for the AC case was defined to be 1.5 whereas for the DC
%         case it was defined to be 1.6.

%   MOST

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
% bus_i	type Pd	Qd Gs  Bs area	Vm	    Va 	baseKV zone	Vmax	Vmin
mpc.bus = [
	1	2	51	27	0	0	1	0.955	10.67	138	1	1.06	0.94;
	2	1	20	9	0	0	1	0.971	11.22	138	1	1.06	0.94;
	3	1	39	10	0	0	1	0.968	11.56	138	1	1.06	0.94;
	4	2	39	12	0	0	1	0.998	15.28	138	1	1.06	0.94;
	5	1	0	0	0	-40	1	1.002	15.73	138	1	1.06	0.94;
	6	2	52	22	0	0	1	0.99	13	    138	1	1.06	0.94;
	7	1	19	2	0	0	1	0.989	12.56	138	1	1.06	0.94;
	8	2	28	0	0	0	1	1.015	20.77	345	1	1.06	0.94;
	9	1	0	0	0	0	1	1.043	28.02	345	1	1.06	0.94;
	10	2	0	0	0	0	1	1.05	35.61	345	1	1.06	0.94;
	11	1	70	23	0	0	1	0.985	12.72	138	1	1.06	0.94;
	12	2	47	10	0	0	1	0.99	12.2	138	1	1.06	0.94;
	13	1	34	16	0	0	1	0.968	11.35	138	1	1.06	0.94;
	14	1	14	1	0	0	1	0.984	11.5	138	1	1.06	0.94;
	15	2	90	30	0	0	1	0.97	11.23	138	1	1.06	0.94;
	16	1	25	10	0	0	1	0.984	11.91	138	1	1.06	0.94;
	17	1	11	3	0	0	1	0.995	13.74	138	1	1.06	0.94;
	18	2	60	34	0	0	1	0.973	11.53	138	1	1.06	0.94;
	19	2	45	25	0	0	1	0.963	11.05	138	1	1.06	0.94;
	20	1	18	3	0	0	1	0.958	11.93	138	1	1.06	0.94;
	21	1	14	8	0	0	2	0.959	13.52	138	1	1.06	0.94;
	22	1	10	5	0	0	2	0.97	16.08	138	1	1.06	0.94;
	23	1	7	3	0	0	2	1	    21	    138	1	1.06	0.94;
	24	2	13	0	0	0	4	0.992	20.89	138	1	1.06	0.94;
	25	2	0	0	0	0	2	1.05	27.93	138	1	1.06	0.94;
	26	2	0	0	0	0	2	1.015	29.71	345	1	1.06	0.94;
	27	2	71	13	0	0	2	0.968	15.35	138	1	1.06	0.94;
	28	1	17	7	0	0	2	0.962	13.62	138	1	1.06	0.94;
	29	1	24	4	0	0	2	0.963	12.63	138	1	1.06	0.94;
	30	1	0	0	0	0	1	0.968	18.79	345	1	1.06	0.94;
	31	2	43	27	0	0	2	0.967	12.75	138	1	1.06	0.94;
	32	2	59	23	0	0	2	0.964	14.8	138	1	1.06	0.94;
	33	1	23	9	0	0	1	0.972	10.63	138	1	1.06	0.94;
	34	2	59	26	0	14	3   0.986	11.3	138	1	1.06	0.94;
	35	1	33	9	0	0	3	0.981	10.87	138	1	1.06	0.94;
	36	2	31	17	0	0	3	0.98	10.87	138	1	1.06	0.94;
	37	1	0	0	0	-25	3	0.992	11.77	138	1	1.06	0.94;
	38	1	0	0	0	0	3	0.962	16.91	345	1	1.06	0.94;
	39	1	27	11	0	0	3	0.97	8.41	138	1	1.06	0.94;
	40	2	66	23	0	0	3	0.97	7.35	138	1	1.06	0.94;
	41	1	37	10	0	0	3	0.967	6.92	138	1	1.06	0.94;
	42	2	96	23	0	0	3	0.985	8.53	138	1	1.06	0.94;
	43	1	18	7	0	0	3	0.978	11.28	138	1	1.06	0.94;
	44	1	16	8	0	10	3	0.985	13.82	138	1	1.06	0.94;
	45	1	53	22	0	10	3	0.987	15.67	138	1	1.06	0.94;
	46	2	28	10	0	10	3	1.005	18.49	138	1	1.06	0.94;
	47	1	34	0	0	0	3	1.017	20.73	138	1	1.06	0.94;
	48	1	20	11	0	15	3	1.021	19.93	138	1	1.06	0.94;
	49	2	87	30	0	0	3	1.025	20.94	138	1	1.06	0.94;
	50	1	17	4	0	0	3	1.001	18.9	138	1	1.06	0.94;
	51	1	17	8	0	0	3	0.967	16.28	138	1	1.06	0.94;
	52	1	18	5	0	0	3	0.957	15.32	138	1	1.06	0.94;
	53	1	23	11	0	0	3	0.946	14.35	138	1	1.06	0.94;
	54	2	113	32	0	0	3	0.955	15.26	138	1	1.06	0.94;
	55	2	63	22	0	0	3	0.952	14.97	138	1	1.06	0.94;
	56	2	84	18	0	0	3	0.954	15.16	138	1	1.06	0.94;
	57	1	12	3	0	0	3	0.971	16.36	138	1	1.06	0.94;
	58	1	12	3	0	0	3	0.959	15.51	138	1	1.06	0.94;
	59	2	277	113	0	0	3	0.985	19.37	138	1	1.06	0.94;
	60	1	78	3	0	0	3	0.993	23.15	138	1	1.06	0.94;
	61	2	0	0	0	0	3	0.995	24.04	138	1	1.06	0.94;
	62	2	77	14	0	0	3	0.998	23.43	138	1	1.06	0.94;
	63	1	0	0	0	0	3	0.969	22.75	345	1	1.06	0.94;
	64	1	0	0	0	0	3	0.984	24.52	345	1	1.06	0.94;
	65	2	0	0	0	0	3	1.005	27.65	345	1	1.06	0.94;
	66	2	39	18	0	0	3	1.05	27.48	138	1	1.06	0.94;
	67	1	28	7	0	0	3	1.02	24.84	138	1	1.06	0.94;
	68	1	0	0	0	0	4	1.003	27.55	345	1	1.06	0.94;
	69	3	0	0	0	0	4	1.035	30	    138	1	1.06	0.94;
	70	2	66	20	0	0	4	0.984	22.58	138	1	1.06	0.94;
	71	1	0	0	0	0	4	0.987	22.15	138	1	1.06	0.94;
	72	2	12	0	0	0	4	0.98	20.98	138	1	1.06	0.94;
	73	2	6	0	0	0	4	0.991	21.94	138	1	1.06	0.94;
	74	2	68	27	0	12	4	0.958	21.64	138	1	1.06	0.94;
	75	1	47	11	0	0	4	0.967	22.91	138	1	1.06	0.94;
	76	2	68	36	0	0	4	0.943	21.77	138	1	1.06	0.94;
	77	2	61	28	0	0	4	1.006	26.72	138	1	1.06	0.94;
	78	1	71	26	0	0	4	1.003	26.42	138	1	1.06	0.94;
	79	1	39	32	0	20	4	1.009	26.72	138	1	1.06	0.94;
	80	2	130	26	0	0	4	1.04	28.96	138	1	1.06	0.94;
	81	1	0	0	0	0	4	0.997	28.1	345	1	1.06	0.94;
	82	1	54	27	0	20	4	0.989	27.24	138	1	1.06	0.94;
	83	1	20	10	0	10	4	0.985	28.42	138	1	1.06	0.94;
	84	1	11	7	0	0	4	0.98	30.95	138	1	1.06	0.94;
	85	2	24	15	0	0	4	0.985	32.51	138	1	1.06	0.94;
	86	1	21	10	0	0	4	0.987	31.14	138	1	1.06	0.94;
	87	2	0	0	0	0	4	1.015	31.4	161	1	1.06	0.94;
	88	1	48	10	0	0	4	0.987	35.64	138	1	1.06	0.94;
	89	2	0	0	0	0	4	1.005	39.69	138	1	1.06	0.94;
	90	2	163	42	0	0	4	0.985	33.29	138	1	1.06	0.94;
	91	2	10	0	0	0	4	0.98	33.31	138	1	1.06	0.94;
	92	2	65	10	0	0	4	0.993	33.8	138	1	1.06	0.94;
	93	1	12	7	0	0	4	0.987	30.79	138	1	1.06	0.94;
	94	1	30	16	0	0	4	0.991	28.64	138	1	1.06	0.94;
	95	1	42	31	0	0	4	0.981	27.67	138	1	1.06	0.94;
	96	1	38	15	0	0	4	0.993	27.51	138	1	1.06	0.94;
	97	1	15	9	0	0	4	1.011	27.88	138	1	1.06	0.94;
	98	1	34	8	0	0	4	1.024	27.4	138	1	1.06	0.94;
	99	2	42	0	0	0	4	1.01	27.04	138	1	1.06	0.94;
	100	2	37	18	0	0	4	1.017	28.03	138	1	1.06	0.94;
	101	1	22	15	0	0	4	0.993	29.61	138	1	1.06	0.94;
	102	1	5	3	0	0	4	0.991	32.3	138	1	1.06	0.94;
	103	2	23	16	0	0	4	1.001	24.44	138	1	1.06	0.94;
	104	2	38	25	0	0	4	0.971	21.69	138	1	1.06	0.94;
	105	2	31	26	0	20	4	0.965	20.57	138	1	1.06	0.94;
	106	1	43	16	0	0	4	0.962	20.32	138	1	1.06	0.94;
	107	2	50	12	0	6	4	0.952	17.53	138	1	1.06	0.94;
	108	1	2	1	0	0	4	0.967	19.38	138	1	1.06	0.94;
	109	1	8	3	0	0	4	0.967	18.93	138	1	1.06	0.94;
	110	2	39	30	0	6	4	0.973	18.09	138	1	1.06	0.94;
	111	2	0	0	0	0	4	0.98	19.74	138	1	1.06	0.94;
	112	2	68	13	0	0	4	0.975	14.99	138	1	1.06	0.94;
	113	2	6	0	0	0	2	0.993	13.74	138	1	1.06	0.94;
	114	1	8	3	0	0	2	0.96	14.46	138	1	1.06	0.94;
	115	1	22	7	0	0	2	0.96	14.46	138	1	1.06	0.94;
	116	2	184	0	0	0	4	1.005	27.12	138	1	1.06	0.94;
	117	1	20	8	0	0	1	0.974	10.67	138	1	1.06	0.94;
	118	1	33	15	0	0	4	0.949	21.92	138	1	1.06	0.94;
];

%% text names of buses (buses with repeated names are usually transformers)

mpc.bus_name = { ...
    'RIVERSIDE','POKAGON','HICKORY CREEK','NEW CARLISLE','OLIVE','KANKAKEE','JACKSON RD','OLIVE','BEQUINE','BREEDER', ...
    'BEND','TWIN BRANCH','CONCORD','GOSHEN JT','FORT WAYNE','NORTH END', 'SORENSON SUB FORT WAYNE','MCKINLEY AVE FORT WAYNE','LINCOLN','ADAMS', ...
    'JAY','RANDOLPH','COLLEGE CORNER','TRENTON','TANNERS CREEK','TANNERS CREEK','MADISON','MULLIN','GRANT','SORENSON SUB FORT WAYNE', ...
    'DEER CREEK','DELAWARE','HAVILAND','ROCK HILL','WEST LIMA','STERLING','EAST LIMA','EAST LIMA','NEW LIBERTY','WEST END', ...
    'SOUTH TIFFIN','HOWARD','SOUTH KENTON','WEST MOUNT VERNON','NORTH NEWARK','WEST LANCASTER','CROOKSVILLE','ZANESVILLE','PHILO','WEST CAMBRIDGE', ...
    'NEW COMERSTOWN','SOUTH COSHOCTON','WOOSTER','TORREY','WAGENHALS','SUNNYSIDE','WEST NEW PHILADELPHIA','WEST NEW PHILADELPHIA 2','TIDD','SOUTH WEST KAMMER', ...
    'KAMMER','NATRIUM','TIDD','KAMMER','MUSKINGUM','MUSKINGUM','SUMMERFIELD','PHILIP SPORN','PHILIP SPORN','PORTSMOUTH', ...
    'NORTH PORTSMOUTH','HILLSBORO','SARGENTS','BELLEFONTE','SOUTH POINT','DARRAH','TURNER','CHEMICAL','CAPITOL HILL','CABIN CREEK', ...
    'KANAWHA','LOGAN','SPRIGG','BETSY LAYNE','BEAVER CREEK','HAZZARD','PINEVILLE','FREMONT','CLINCH RIVER','HOLSTON', ...
    'HOLSTON T','SALTVILLE','TAZEWELL','SWITCHBACK','CALDWELL','BAILEYSVILLE','SUNDIAL','BRADLEY','HINTON','GLEN LYN', ...
    'WYTHE','SMYTHE','CLAYTOR','HANCOCK','ROANOKE','CLOVERDALE','REUSENS','BLAINE','FRANKLIN','FIELDALE', ...
    'DAN RIVER','DANVILLE','DEER CREEK','WEST MEDFORD','MEDFORD','KYGER CREEK','COREY','WEST HUNTINGTON' ...
    }';

%% generator data
% bus  Pg  Qg  Qmax	Qmin Vg	 mBase sts Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
mpc.gen = [
% 	1	0	0	15	-5	0.955	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
% 	4	0	0	300	-300 0.998	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	6	0	0	250	-50	0.99	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;  % Generator critical for voltage support on bus 6
	8	0	0	300	-300 1.015	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;  % in higher load conditions
	10	450	0	200	-147 1.05	100	1 550	55	0	0	0	0	0	0       110	110	110	0	0;
	12	85	0	120	-35	 0.99	100	1 185	20	0	0	0	0	0	0       37	37	37	0	0;
% 	15	0	0	30	-10	 0.97	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
% 	18	0	0	50	-16	 0.973	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
% 	19	0	0	24	-8	 0.962	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	24	0	0	300	-300 0.992	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
	25	220	0	140	-47	 1.05	100	1 320	30	0	0	0	0	0	0       64	64  64	0	0;
	26	314	0  1000 -1000 1.015 100	1 414	40	0	0	0	0	0	0       83	83	83	0	0;
	27	0	0	300	-300 0.968	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
% 	31	7	0	300	-300 0.967	100	1 107	0	0	0	0	0	0	0       21	21	21	0	0;
% 	32	0	0	42	-14	 0.963	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
% 	34	0	0	24	-8	 0.984	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
% 	36	0	0	24	-8	 0.98	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	40	0	0	300	-300 0.97	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
	42	0	0	300	-300 0.985	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
	46	19	0	100	-100 1.005	100	1 119	12	0	0	0	0	0	0       24	24	24	0	0;
	49	204	0	210	-85	 1.025	100	1 304	30	0	0	0	0	0	0       61	61	61	0	0;
	54	48	0	300	-300 0.955	100	1 148	15	0	0	0	0	0	0       30	30	30	0	0;
% 	55	0	0	23	-8	 0.952	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
% 	56	0	0	15	-8	 0.954	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	59	155	0	180	-60	 0.985	100	1 255	25	0	0	0	0	0	0       51	51	51	0	0;
	61	160	0	300	-100 0.995	100	1 260	25	0	0	0	0	0	0       52	52	52	0	0;
	62	0	0	20	-20	 0.998	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
	65	391	0	200	-67	 1.005	100	1 491	50	0	0	0	0	0	0       98	98	98	0	0;
	66	392	0	200	-67	 1.05	100	1 492	50	0	0	0	0	0	0       98	98	98	0	0;
	69	516.4 0	300 -300 1.035	100	1 805.2 80	0	0	0	0	0	0       161	161	161	0	0;
% 	70	0	0	32	-10	 0.984	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	72	0	0	100	-100 0.98	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
	73	0	0	100	-100 0.991	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
% 	74	0	0	9	-6	 0.958	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	76	0	0	23	-8	 0.943	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;  % These two critical for voltage support of bus 76
	77	0	0	330	-50	 1.006	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;  % in higher load conditions
	80	477	0	280	-165 1.04	100	1 577	60	0	0	0	0	0	0       115	115	115	0	0;
% 	85	0	0	23	-8	 0.985	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	87	4	0	1000 -100 1.015	100	1 104	10	0	0	0	0	0	0       20	20	20	0	0;
	89	607	0	300	-210 1.005	100	1 707	70	0	0	0	0	0	0       141	141	141	0	0;
	90	0	0	300	-300 0.985	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
	91	0	0	100	-100 0.98	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
% 	92	0	0	9	-3	 0.99	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	99	0	0	100	-100 1.01	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	100	252	0	155	-50	 1.017	100	1 352	35	0	0	0	0	0	0       70	70	70	0	0;
	103	40	0	40	-15	 1.01	100	1 140	0	0	0	0	0	0	0       28	28	28	0	0;
	104	0	0	23	-8	 0.971	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
% 	105	0	0	23	-8	 0.965	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	107	0	0	200	-200 0.952	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	110	0	0	23	-8	 0.973	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	111	36	0	1000 -100 0.98	100	1 136	14	0	0	0	0	0	0       27	27	27	0	0;
% 	112	0	0	1000 -100 0.975	100	1 100	0	0	0	0	0	0	0       20	20	20	0	0;
	113	0	0	200	-100 0.993	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
	116	0	0  1000 -1000 1.005	100	1 100	10	0	0	0	0	0	0       20	20	20	0	0;
    % from here on, wind units
    1   0   0   200  -100 1.00  100 1 65   0   0   0   0   0   0   0      65 65 65  0   0; % TIFFIN
    6   0   0  2000 -1000 1.035 100 1 626  0   0   0   0   0   0   0     626 626 626 0   0; %
    9   0   0  2000 -1000 1.035 100 1 742  0   0   0   0   0   0   0     742 742 742 0   0; %
   18   0   0   200  -100 1.035 100 1 212  0   0   0   0   0   0   0     212 212 212 0   0; %
   19   0   0   200  -100 1.035 100 1  63  0   0   0   0   0   0   0      63  63  63 0   0; %
   41   0   0  2000 -1000 1.035 100 1 512  0   0   0   0   0   0   0     512 512 512 0   0; %
   43   0   0   200  -100 1.035 100 1  66  0   0   0   0   0   0   0      66  66  66 0   0; %
   62   0   0   500  -300 1.035 100 1 153  0   0   0   0   0   0   0     153 153 153 0   0; %
   63   0   0   600  -300 1.035 100 1  87  0   0   0   0   0   0   0      87  87  87 0   0; %
   72   0   0   200  -100 1.035 100 1  45  0   0   0   0   0   0   0      45  45  45 0   0; %
   80   0   0   200  -100 1.035 100 1  40  0   0   0   0   0   0   0      40  40  40 0   0; %
    % From here on, storage units
    59  0   0    10   -10 1.00  100 1 400 -400  0   0   0   0   0   0       20  50 100  0   0; % TIDD
   116  0   0    10   -10 1.035 100 1 400 -400  0   0   0   0   0   0       20  50 100  0   0; % KYGER CREEK
    90  0   0    10   -10 1.035 100 1 400 -400  0   0   0   0   0   0       20  50 100  0   0; % HOLSTON
    80  0   0    10   -10 1.035 100 1 400 -400  0   0   0   0   0   0       20  50 100  0   0; % CABIN CREEK
];

%% branch data
% fbus tbus	  r	      x	     b	  rateA	rateB rateC ratio angle status	angmin angmax
mpc.branch = [
	1	2	0.0303	0.0999	0.0254	175	0	0	0	0	1	-360	360;
	1	3	0.0129	0.0424	0.01082	175 0	0	0	0	1	-360	360;
	4	5	0.00176	0.00798	0.0021	175	0	0	0	0	1	-360	360;
	3	5	0.0241	0.108	0.0284	175	0	0	0	0	1	-360	360;
	5	6	0.0119	0.054	0.01426	175	0	0	0	0	1	-360	360;
	6	7	0.00459	0.0208	0.0055	175	0	0	0	0	1	-360	360;
	8	9	0.0049  0.0610  0.5810  440	0	0	0	0	1	-360	360;%HD
	8	9	0.0049  0.0610  0.5810	440	0	0	0	0	1	-360	360;%HD
	8	5	0	    0.0267	0	    550 0	0 0.985	0	1	-360	360;%trafo
	9	10	0.0052  0.0644  0.6150	440	0	0	0	0	1	-360	360;%HD
	9	10	0.0052  0.0644  0.6150	440	0	0	0	0	1	-360	360;%HD
	4	11	0.0209	0.0688	0.01748	175	0	0	0	0	1	-360	360;
	5	11	0.0203	0.0682	0.01738	175	0	0	0	0	1	-360	360;
	11	12	0.00595	0.0196	0.00502	175	0	0	0	0	1	-360	360;
	2	12	0.0187	0.0616	0.01572	175	0	0	0	0	1	-360	360;
	3	12	0.0484	0.16	0.0406	175	0	0	0	0	1	-360	360;
	7	12	0.00862	0.034	0.00874	175	0	0	0	0	1	-360	360;
	11	13	0.02225	0.0731	0.01876	175	0	0	0	0	1	-360	360;
	12	14	0.0215	0.0707	0.01816	175	0	0	0	0	1	-360	360;
	13	15	0.0744	0.2444	0.06268	175	0	0	0	0	1	-360	360;
	14	15	0.0595	0.195	0.0502	175	0	0	0	0	1	-360	360;
	12	16	0.0212	0.0834	0.0214	175	0	0	0	0	1	-360	360;
	15	17	0.0264  0.0874  0.0222	175	0	0	0	0	1	-360	360;%Doub C
	15	17	0.0264  0.0874  0.0222	175	0	0	0	0	1	-360	360;%Doub C
	16	17	0.0454	0.1801	0.0466	175	0	0	0	0	1	-360	360;
	17	18	0.0123	0.0505	0.01298	175	0	0	0	0	1	-360	360;
	18	19	0.01119	0.0493	0.01142	175	0	0	0	0	1	-360	360;
	19	20	0.0252	0.117	0.0298	175	0	0	0	0	1	-360	360;
	15	19	0.012	0.0394	0.0101	175	0	0	0	0	1	-360	360;
	20	21	0.0183	0.0849	0.0216	175	0	0	0	0	1	-360	360;
	21	22	0.0209	0.097	0.0246	175	0	0	0	0	1	-360	360;
	22	23	0.0342	0.159	0.0404	175	0	0	0	0	1	-360	360;
	23	24	0.0270  0.0984  0.0249	175	0	0	0	0	1	-360	360;%D circ
	23	24	0.0270  0.0984  0.0249	175	0	0	0	0	1	-360	360;%D circ
	23	25	0.0312  0.1600  0.0432	175	0	0	0	0	1	-360	360;%D circ
	23	25	0.0312  0.1600  0.0432	175	0	0	0	0	1	-360	360;%D circ
	26	25	0       0.0382	0	    200	0	0  0.96 0	1	-360	360;%trafo
	25	27	0.0636  0.3260  0.0882	175	0	0	0	0	1	-360	360;%D circ
	25	27	0.0636  0.3260  0.0882	175	0	0	0	0	1	-360	360;%D circ
	27	28	0.01913	0.0855	0.0216	175	0	0	0	0	1	-360	360;
	28	29	0.0237	0.0943	0.0238	175	0	0	0	0	1	-360	360;
	30	17	0	    0.0388	0   	300	0	0 0.96	0	1	-360	360;%trafo
	8	30	0.00431	0.0504	0.514	440	0	0	0	0	1	-360	360;
	26	30	0.00799	0.086	0.908	440	0	0	0	0	1	-360	360;
	17	31	0.0474	0.1563	0.0399	175	0	0	0	0	1	-360	360;
	29	31	0.0108	0.0331	0.0083	175	0	0	0	0	1	-360	360;
	23	32	0.0634  0.2306  0.0587	175	0	0	0	0	1	-360	360;%D circ
	23	32	0.0634  0.2306  0.0587	175	0	0	0	0	1	-360	360;%D circ
	31	32	0.0298	0.0985	0.0251	175	0	0	0	0	1	-360	360;
	27	32	0.0229	0.0755	0.01926	175	0	0	0	0	1	-360	360;
	15	33	0.038	0.1244	0.03194	175	0	0	0	0	1	-360	360;
	19	34	0.0752	0.247	0.0632	175	0	0	0	0	1	-360	360;
	35	36	0.00224	0.0102	0.00268	175	0	0	0	0	1	-360	360;
	35	37	0.011	0.0497	0.01318	175	0	0	0	0	1	-360	360;
	33	37	0.0415	0.142	0.0366	175	0	0	0	0	1	-360	360;
	34	36	0.00871	0.0268	0.00568	175	0	0	0	0	1	-360	360;
	34	37	0.0051  0.0188  0.0049	175	0	0	0	0	1	-360	360;%D CIRC
	34	37	0.0051  0.0188  0.0049	175	0	0	0	0	1	-360	360;%D CIRC
	38	37	0	    0.0375	0	    350	0	0 0.935	0	1	-360	360;%trafo
	37	39	0.0321	0.106	0.027	175	0	0	0	0	1	-360	360;
	37	40	0.0593	0.168	0.042	175	0	0	0	0	1	-360	360;
	30	38	0.00464	0.054	0.422	440	0	0	0	0	1	-360	360;
	39	40	0.0184	0.0605	0.01552	175	0	0	0	0	1	-360	360;
	40	41	0.0145	0.0487	0.01222	175	0	0	0	0	1	-360	360;
	40	42	0.0555	0.183	0.0466	175	0	0	0	0	1	-360	360;
	41	42	0.041	0.135	0.0344	175	0	0	0	0	1	-360	360;
	43	44	0.0608	0.2454	0.06068	175	0	0	0	0	1	-360	360;
	34	43	0.0413	0.1681	0.04226	175	0	0	0	0	1	-360	360;
	44	45	0.0224	0.0901	0.0224	175	0	0	0	0	1	-360	360;
	45	46	0.04	0.1356	0.0332	175	0	0	0	0	1	-360	360;
	46	47	0.0760  0.2540  0.0158	175	0	0	0	0	1	-360	360;%D CIRC
	46	47	0.0760  0.2540  0.0158	175	0	0	0	0	1	-360	360;%D CIRC
	46	48	0.0601	0.189	0.0472	175	0	0	0	0	1	-360	360;
	47	49	0.0191	0.0625	0.01604	175	0	0	0	0	1	-360	360;
	42	49	0.0715	0.323	0.086	175	0	0	0	0	1	-360	360;
	42	49	0.0715	0.323	0.086	175	0	0	0	0	1	-360	360;
	45	49	0.0684	0.186	0.0444	175	0	0	0	0	1	-360	360;
	48	49	0.0179	0.0505	0.01258	175	0	0	0	0	1	-360	360;
	49	50	0.0267	0.0752	0.01874	175	0	0	0	0	1	-360	360;
	49	51	0.0486	0.137	0.0342	175	0	0	0	0	1	-360	360;
	51	52	0.0203	0.0588	0.01396	175	0	0	0	0	1	-360	360;
	52	53	0.0405	0.1635	0.04058	175	0	0	0	0	1	-360	360;
	53	54	0.0263	0.122	0.031	175	0	0	0	0	1	-360	360;
	49	54	0.073	0.289	0.0738	175	0	0	0	0	1	-360	360;
	49	54	0.0869	0.291	0.073	175	0	0	0	0	1	-360	360;
	54	55	0.0169	0.0707	0.0202	175	0	0	0	0	1	-360	360;
	54	56	0.0055  0.0191  0.0037	175	0	0	0	0	1	-360	360;%D circ
	54	56	0.0055  0.0191  0.0037	175	0	0	0	0	1	-360	360;%D circ
	55	56	0.00488	0.0151	0.00374	175	0	0	0	0	1	-360	360;
	56	57	0.0343	0.0966	0.0242	175	0	0	0	0	1	-360	360;
	50	57	0.0474	0.134	0.0332	175	0	0	0	0	1	-360	360;
	56	58	0.0343	0.0966	0.0242	175	0	0	0	0	1	-360	360;
	51	58	0.0255	0.0719	0.01788	175	0	0	0	0	1	-360	360;
	54	59	0.0503	0.2293	0.0598	175	0	0	0	0	1	-360	360;
	56	59	0.0825	0.251	0.0569	175	0	0	0	0	1	-360	360;
	56	59	0.0803	0.239	0.0536	175	0	0	0	0	1	-360	360;
	55	59	0.04739	0.2158	0.05646	175	0	0	0	0	1	-360	360;
	59	60	0.0317	0.145	0.0376	175	0	0	0	0	1	-360	360;
	59	61	0.0328	0.15	0.0388	175	0	0	0	0	1	-360	360;
	60	61	0.0053  0.0270  0.0073	175	0	0	0	0	1	-360	360;%D circ
	60	61	0.0053  0.0270  0.0073	175	0	0	0	0	1	-360	360;%D circ
	60	62	0.0123	0.0561	0.01468	175	0	0	0	0	1	-360	360;
	61	62	0.00824	0.0376	0.0098	175	0	0	0	0	1	-360	360;
	63	59	0	    0.0386	0	    400	0	0 0.96	0	1	-360	360;%trafo
	63	64	0.00172	0.02	0.216	440	0	0	0	0	1	-360	360;%H
	64	61	0   	0.0268	0	    200	0	0 0.985	0	1	-360	360;%trafo
	38	65	0.00901	0.0986	1.046	440	0	0	0	0	1	-360	360;%H
	64	65	0.00269	0.0302	0.38	440	0	0	0	0	1	-360	360;%H
	49	66	0.018	0.0919	0.0248	175	0	0	0	0	1	-360	360;
	49	66	0.018	0.0919	0.0248	175	0	0	0	0	1	-360	360;
	62	66	0.0482	0.218	0.0578	175	0	0	0	0	1	-360	360;
	62	67	0.0258	0.117	0.031	175	0	0	0	0	1	-360	360;
	65	66	0	    0.037	0	    200	0	0 0.935	0	1	-360	360;%trafo
	66	67	0.0224	0.1015	0.02682	175	0	0	0	0	1	-360	360;
	65	68	0.0028  0.0320  0.3190	440	0	0	0	0	1	-360	360;%HD
	65	68	0.0028  0.0320  0.3190	440	0	0	0	0	1	-360	360;%HD
	47	69	0.0844	0.2778	0.07092	175	0	0	0	0	1	-360	360;
	49	69	0.0985	0.324	0.0828	175	0	0	0	0	1	-360	360;
	68	69	0	    0.037	0	    350	0	0 0.935	0	1	-360	360;%trafo
	69	70	0.0600  0.2540  0.0610	175	0	0	0	0	1	-360	360;%D circ
	69	70	0.0600  0.2540  0.0610	175	0	0	0	0	1	-360	360;%D circ
	24	70	0.00221	0.4115	0.10198	175	0	0	0	0	1	-360	360;
	70	71	0.00882	0.0355	0.00878	175	0	0	0	0	1	-360	360;
	24	72	0.0488	0.196	0.0488	175	0	0	0	0	1	-360	360;
	71	72	0.0446	0.18	0.04444	175	0	0	0	0	1	-360	360;
	71	73	0.00866	0.0454	0.01178	175	0	0	0	0	1	-360	360;
	70	74	0.0401	0.1323	0.03368	175	0	0	0	0	1	-360	360;
	70	75	0.0428	0.141	0.036	175	0	0	0	0	1	-360	360;
	69	75	0.0810  0.2440  0.0620	175	0	0	0	0	1	-360	360;%D circ
	69	75	0.0810  0.2440  0.0620	175	0	0	0	0	1	-360	360;%D circ
	74	75	0.0123	0.0406	0.01034	175	0	0	0	0	1	-360	360;
	76	77	0.0444	0.148	0.0368	175	0	0	0	0	1	-360	360;
	69	77	0.0618  0.2020  0.0519	175	0   0	0	0	1	-360	360;%D circ
	69	77	0.0618  0.2020  0.0519	175	0   0	0	0	1	-360	360;%D circ
	75	77	0.0601	0.1999	0.04978	175	0	0	0	0	1	-360	360;
	77	78	0.0075  0.0248  0.0063	175	0	0	0	0	1	-360	360;%D circ
	77	78	0.0075  0.0248  0.0063	175	0	0	0	0	1	-360	360;%D circ
	78	79	0.00546	0.0244	0.00648	175	0	0	0	0	1	-360	360;
	77	80	0.017	0.0485	0.0472	175	0	0	0	0	1	-360	360;
	77	80	0.0294	0.105	0.0228	175	0	0	0	0	1	-360	360;
	79	80	0.0156	0.0704	0.0187	175	0	0	0	0	1	-360	360;
	68	81	0.0035  0.0404  0.4040	440	0	0	0	0	1	-360	360;%HD
	68	81	0.0035  0.0404  0.4040	440	0	0	0	0	1	-360	360;%HD
	81	80	0	    0.037	0	    200	0   0 0.935	0	1	-360	360;%Trafo
	77	82	0.0596  0.1706  0.0409	175	0	0	0	0	1	-360	360;%D circ
	77	82	0.0596  0.1706  0.0409	175	0	0	0	0	1	-360	360;%D circ
	82	83	0.0224  0.0733  0.0190	175	0	0	0	0	1	-360	360;%D circ
	82	83	0.0224  0.0733  0.0190	175	0	0	0	0	1	-360	360;%D circ
	83	84	0.0625	0.132	0.0258	175	0	0	0	0	1	-360	360;
	83	85	0.043	0.148	0.0348	175	0	0	0	0	1	-360	360;
	84	85	0.0302	0.0641	0.01234	175	0	0	0	0	1	-360	360;
	85	86	0.035	0.123	0.0276	175	0	0	0	0	1	-360	360;
	86	87	0.02828	0.2074	0.0445	200	0	0	0	0	1	-360	360;%161Kv
	85	88	0.02	0.102	0.0276	175	0	0	0	0	1	-360	360;
	85	89	0.0239	0.173	0.047	175	0	0	0	0	1	-360	360;
	88	89	0.0139	0.0712	0.01934	175	0	0	0	0	1	-360	360;
	89	90	0.0518	0.188	0.0528	175	0	0	0	0	1	-360	360;
	89	90	0.0238	0.0997	0.106	175	0	0	0	0	1	-360	360;
	90	91	0.0254	0.0836	0.0214	175	0	0	0	0	1	-360	360;
	89	92	0.0099	0.0505	0.0548	175	0	0	0	0	1	-360	360;
	89	92	0.0393	0.1581	0.0414	175	0	0	0	0	1	-360	360;
	91	92	0.0387	0.1272	0.03268	175	0	0	0	0	1	-360	360;
	92	93	0.0258	0.0848	0.0218	175	0	0	0	0	1	-360	360;
	92	94	0.0481	0.158	0.0406	175	0	0	0	0	1	-360	360;
	93	94	0.0223	0.0732	0.01876	175	0	0	0	0	1	-360	360;
	94	95	0.0132	0.0434	0.0111	175	0	0	0	0	1	-360	360;
	80	96	0.0356	0.182	0.0494	175	0	0	0	0	1	-360	360;
	82	96	0.0324  0.1060  0.0272	175	0	0	0	0	1	-360	360;%D circ
	82	96	0.0324  0.1060  0.0272	175	0	0	0	0	1	-360	360;%D circ
	94	96	0.0269	0.0869	0.023	175	0	0	0	0	1	-360	360;
	80	97	0.0183	0.0934	0.0254	175	0	0	0	0	1	-360	360;
	80	98	0.0238	0.108	0.0286	175	0	0	0	0	1	-360	360;
	80	99	0.0454	0.206	0.0546	175	0	0	0	0	1	-360	360;
	92	100	0.0648	0.295	0.0472	175	0	0	0	0	1	-360	360;
	94	100	0.0356  0.1160  0.0302	175	0	0	0	0	1	-360	360;%D circ
	94	100	0.0356  0.1160  0.0302	175	0	0	0	0	1	-360	360;%D circ
	95	96	0.0171	0.0547	0.01474	175	0	0	0	0	1	-360	360;
	96	97	0.0173	0.0885	0.024	175	0	0	0	0	1	-360	360;
	98	100	0.0397	0.179	0.0476	175	0	0	0	0	1	-360	360;
	99	100	0.018	0.0813	0.0216	175	0	0	0	0	1	-360	360;
	100	101	0.0277	0.1262	0.0328	175	0	0	0	0	1	-360	360;
	92	102	0.0123	0.0559	0.01464	175	0	0	0	0	1	-360	360;
	101	102	0.0246	0.112	0.0294	175	0	0	0	0	1	-360	360;
	100	103	0.0320  0.1050  0.0268	175	0	0	0	0	1	-360	360;%D circ
	100	103	0.0320  0.1050  0.0268	175	0	0	0	0	1	-360	360;%D circ
	100	104	0.0451	0.204	0.0541	175	0	0	0	0	1	-360	360;
	103	104	0.0466	0.1584	0.0407	175	0	0	0	0	1	-360	360;
	103	105	0.0535	0.1625	0.0408	175	0	0	0	0	1	-360	360;
	100	106	0.0605	0.229	0.062	175	0	0	0	0	1	-360	360;
	104	105	0.00994	0.0378	0.00986	175	0	0	0	0	1	-360	360;
	105	106	0.014	0.0547	0.01434	175	0	0	0	0	1	-360	360;
	105	107	0.053	0.183	0.0472	175	0	0	0	0	1	-360	360;
	105	108	0.0261	0.0703	0.01844	175	0	0	0	0	1	-360	360;
	106	107	0.053	0.183	0.0472	175	0	0	0	0	1	-360	360;
	108	109	0.0105	0.0288	0.0076	175	0	0	0	0	1	-360	360;
	103	110	0.03906	0.1813	0.0461	175	0	0	0	0	1	-360	360;
	109	110	0.0278	0.0762	0.0202	175	0	0	0	0	1	-360	360;
	110	111	0.022	0.0755	0.02	175	0	0	0	0	1	-360	360;
	110	112	0.0494  0.1280  0.0310	175	0	0	0	0	1	-360	360;%D circ
	110	112	0.0494  0.1280  0.0310	175	0	0	0	0	1	-360	360;%D circ
	17	113	0.00913	0.0301	0.00768	175	0	0	0	0	1	-360	360;
	32	113	0.0615	0.203	0.0518	175	0	0	0	0	1	-360	360;
	32	114	0.0135	0.0612	0.01628	175	0	0	0	0	1	-360	360;
	27	115	0.0164	0.0741	0.01972	175	0	0	0	0	1	-360	360;
	114	115	0.0023	0.0104	0.00276	175	0	0	0	0	1	-360	360;
	68	116	0.0007  0.0081  0.0820	440	0	0	0	0	1	-360	360;%HD
    68	116	0.0007  0.0081  0.0820	440	0	0	0	0	1	-360	360;%HD
	12	117	0.0329	0.14	0.0358	175	0	0	0	0	1	-360	360;
	75	118	0.0145	0.0481	0.01198	175	0	0	0	0	1	-360	360;
	76	118	0.0164	0.0544	0.01356	175	0	0	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown n	x1	y1	...	xn yn % For piecewise linear case
%	2	startup	shutdown	n	c(n-1)	...	c0 %for polynomial case
mpc.gencost = [
% 	2	0	0	2	 35.39 	0	0	0	0	0;	%
% 	2	0	0	2	 31.45 	0	0	0	0	0;	%
	2	2000	0	2	 20.22 	0	0	0	0	0;	%
	2	3000	0	2	 28.40 	300	0	0	0	0;	%
	2   16000	0	2	 29.36 	1600	0	0	0	0;	%
	2	7000	0	2	 36.85 	70	0	0	0	0;	%
% 	2	0	0	2	 45.58 	0	0	0	0	0;	%
% 	2	0	0	2	 213.32 	0	0	0	0	0;	%
% 	2	0	0	2	 62.15 	0	0	0	0	0;	%
	2	2700	0	2	 26.64 	300	0	0	0	0;	%
	2	8000	0	2	 26.64 	800	0	0	0	0;	%
	2	11000	0	2	 26.64 	1100	0	0	0	0;	%
	2	5000	0	2	 45.58 	500	0	0	0	0;	%
% 	2	0	0	2	 31.77 	0	0	0	0	0;	%
% 	2	0	0	2	 37.34 	0	0	0	0	0;	%
% 	2	0	0	2	 212.03 	0	0	0	0	0;	%
% 	2	0	0	2	 329.73 	0	0	0	0	0;	%
	2	4000	0	2	 37.34 	400	0	0	0	0;	%
	2	3000	0	2	 33.39 	300	0	0	0	0;	%
	2	4000	0	2	 33.39 	400	0	0	0	0;	%
	2	8000	0	2	 25.87 	800	0	0	0	0;	%
	2	4000	0	2	 53.21 	200	0	0	0	0;	%
% 	2	0	0	2	 2.00 	0	0	0	0	0;	%
% 	2	0	0	2	 612.62 	0	0	0	0	0;	%
	2	7000	0	2	 26.57 	700	0	0	0	0;	%
	2	8000	0	2	 31.77 	800	0	0	0	0;	%
	2	2000	0	2	 34.60 	100	0	0	0	0;	%
	2	12000	0	2	 24.44 	1200	0	0	0	0;	%
	2	15000	0	2	 30.81 	1500	0	0	0	0;	%
	2	40000	0	2	 48.13 	4000	0	0	0	0;	%
% 	2	0	0	2	 1.29 	0	0	0	0	0;	%
	2	2500	0	2	 24.44 	250	0	0	0	0;	%
	2	2500	0	2	 24.01 	250	0	0	0	0;	%
% 	2	0	0	2	 25.87 	0	0	0	0	0;	%
	2	20  	0	2	  1.29 	0	0	0	0	0;	%
	2	3200	0	2	 31.45 	0	0	0	0	0;	%
	2	30000	0	2	 45.66 	300	0	0	0	0;	%
% 	2	0	0	2	 66.19 	0	0	0	0	0;	%
	2	3200	0	2	 30.21 	320	0	0	0	0;	%
	2	22000	0	2	 30.09 	2200	0	0	0	0;	%
	2	3400	0	2	 33.39 	340	0	0	0	0;	%
	2	3000	0	2	 30.09 	300	0	0	0	0;	%
% 	2	0	0	2	 215.89 	0	0	0	0	0;	%
	2	0	0	2	 1.29 	0	0	0	0	0;	%
	2	10000	0	2	 28.99 	1000	0	0	0	0;	%
	2	0	0	2	 2.00 	0	0	0	0	0;	%
	2	0	0	2	 2.00 	0	0	0	0	0;	%
% 	2	0	0	2	 1.29 	0	0	0	0	0;	%
	2	0	0	2	 1.29 	0	0	0	0	0;	%
	2	0	0	2	 1.29 	0	0	0	0	0;	%
	2	5000	0	2	 33.39 	500	0	0	0	0;	%
% 	2	0	0	2	 59.82 	0	0	0	0	0;	%
	2	3100	0	2	 30.87 	300	0	0	0	0;	%
	2	2700	0	2	 26.64 	300	0	0	0	0;	%

    % wind units
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;

    % storage units
    %1   0   0   3   -200   2    0    0    200    2;
    %1   0   0   3   -200   2    0    0    200    2;
	%2	0	0	2	0.0	 0.0	0    0      0    0;
	%2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2   0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
	2	0	0	2   0.0	 0.0	0    0      0    0;
	2	0	0	2	0.0	 0.0	0    0      0    0;
];

%% generator fuel type
mpc.genfuel = { ...
% 	    'syncgen'	;
% 	    'syncgen'	;
	    'syncgen'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
% 	    'syncgen'	;
% 	    'oil'	;
% 	    'syncgen'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
% 	    'syncgen'	;
% 	    'syncgen'	;
% 	    'syncgen'	;
% 	    'syncgen'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
	    'ng'	;
% 	    'syncgen'	;
% 	    'syncgen'	;
	    'coal'	;
	    'coal'	;
	    'ng'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
% 	    'syncgen'	;
	    'coal'	;
	    'coal'	;
% 	    'syncgen'	;
	    'syncgen'	;
	    'syncgen'	;
	    'coal'	;
% 	    'syncgen'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
	    'coal'	;
% 	    'oil'	;
	    'hydro'	;
	    'coal'	;
	    'hydro'	;
	    'hydro'	;
% 	    'syncgen'	;
	    'hydro'	;
	    'hydro'	;
	    'coal'	;
% 	    'syncgen'	;
	    'coal'	;
	    'coal'	;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'wind'  ;
	    'ess'   ;
	    'ess'   ;
	    'ess'   ;
	    'ess'   ;
};

icoal = [];
ioil = [];
ing = [];
ihydro = [];
inuke = [];
iwind = [];
iess = [];
isyncgen = [];
for i = 1:length(mpc.genfuel)
  if strcmp('coal', mpc.genfuel{i})
    icoal = [ icoal; i];
  elseif strcmp('oil', mpc.genfuel{i})
    ioil = [ ioil; i];
  elseif strcmp('hydro', mpc.genfuel{i})
    ihydro = [ ihydro; i];
  elseif strcmp('ng', mpc.genfuel{i})
    ing = [ ing; i];
  elseif strcmp('nuke', mpc.genfuel{i})
    inuke = [ inuke; i];
  elseif strcmp('wind', mpc.genfuel{i})
    iwind = [ iwind; i];
  elseif strcmp('ess', mpc.genfuel{i})
    iess = [ iess; i];
  elseif strcmp('syncgen', mpc.genfuel{i})
    isyncgen = [ isyncgen; i];
  else
    error('unrecognized fuel type')
  end
end
mpc.icoal = icoal;
mpc.ioil = ioil;
mpc.ing = ing;
mpc.ihydro = ihydro;
mpc.inuke = inuke;
mpc.iwind = iwind;
mpc.iess = iess;
mpc.isyncgen = isyncgen;

define_constants;
mpc.NamePlateWind = mpc.gen(iwind, PMAX);

% Warnings from cdf2matp conversion:
%
% ***** area data conversion not yet implemented (creating dummy area data)
% ***** negative Pg at bus 4 treated as Pd
% ***** negative Pg at bus 8 treated as Pd
% ***** negative Pg at bus 24 treated as Pd
% ***** negative Pg at bus 27 treated as Pd
% ***** negative Pg at bus 40 treated as Pd
% ***** negative Pg at bus 42 treated as Pd
% ***** negative Pg at bus 72 treated as Pd
% ***** negative Pg at bus 73 treated as Pd
% ***** negative Pg at bus 90 treated as Pd
% ***** negative Pg at bus 91 treated as Pd
% ***** negative Pg at bus 99 treated as Pd
% ***** negative Pg at bus 107 treated as Pd
% ***** negative Pg at bus 112 treated as Pd
% ***** negative Pg at bus 113 treated as Pd
% ***** negative Pg at bus 116 treated as Pd
% ***** Insufficient generation, setting Pmax at slack bus (bus 69) to 805.2
% ***** MVA limit of branch 1 - 2 not given, set to 9900
% ***** MVA limit of branch 1 - 3 not given, set to 9900
% ***** MVA limit of branch 4 - 5 not given, set to 9900
% ***** MVA limit of branch 3 - 5 not given, set to 9900
% ***** MVA limit of branch 5 - 6 not given, set to 9900
% ***** MVA limit of branch 6 - 7 not given, set to 9900
% ***** MVA limit of branch 8 - 9 not given, set to 9900
% ***** MVA limit of branch 8 - 5 not given, set to 9900
% ***** MVA limit of branch 9 - 10 not given, set to 9900
% ***** MVA limit of branch 4 - 11 not given, set to 9900
% ***** MVA limit of branch 5 - 11 not given, set to 9900
% ***** MVA limit of branch 11 - 12 not given, set to 9900
% ***** MVA limit of branch 2 - 12 not given, set to 9900
% ***** MVA limit of branch 3 - 12 not given, set to 9900
% ***** MVA limit of branch 7 - 12 not given, set to 9900
% ***** MVA limit of branch 11 - 13 not given, set to 9900
% ***** MVA limit of branch 12 - 14 not given, set to 9900
% ***** MVA limit of branch 13 - 15 not given, set to 9900
% ***** MVA limit of branch 14 - 15 not given, set to 9900
% ***** MVA limit of branch 12 - 16 not given, set to 9900
% ***** MVA limit of branch 15 - 17 not given, set to 9900
% ***** MVA limit of branch 16 - 17 not given, set to 9900
% ***** MVA limit of branch 17 - 18 not given, set to 9900
% ***** MVA limit of branch 18 - 19 not given, set to 9900
% ***** MVA limit of branch 19 - 20 not given, set to 9900
% ***** MVA limit of branch 15 - 19 not given, set to 9900
% ***** MVA limit of branch 20 - 21 not given, set to 9900
% ***** MVA limit of branch 21 - 22 not given, set to 9900
% ***** MVA limit of branch 22 - 23 not given, set to 9900
% ***** MVA limit of branch 23 - 24 not given, set to 9900
% ***** MVA limit of branch 23 - 25 not given, set to 9900
% ***** MVA limit of branch 26 - 25 not given, set to 9900
% ***** MVA limit of branch 25 - 27 not given, set to 9900
% ***** MVA limit of branch 27 - 28 not given, set to 9900
% ***** MVA limit of branch 28 - 29 not given, set to 9900
% ***** MVA limit of branch 30 - 17 not given, set to 9900
% ***** MVA limit of branch 8 - 30 not given, set to 9900
% ***** MVA limit of branch 26 - 30 not given, set to 9900
% ***** MVA limit of branch 17 - 31 not given, set to 9900
% ***** MVA limit of branch 29 - 31 not given, set to 9900
% ***** MVA limit of branch 23 - 32 not given, set to 9900
% ***** MVA limit of branch 31 - 32 not given, set to 9900
% ***** MVA limit of branch 27 - 32 not given, set to 9900
% ***** MVA limit of branch 15 - 33 not given, set to 9900
% ***** MVA limit of branch 19 - 34 not given, set to 9900
% ***** MVA limit of branch 35 - 36 not given, set to 9900
% ***** MVA limit of branch 35 - 37 not given, set to 9900
% ***** MVA limit of branch 33 - 37 not given, set to 9900
% ***** MVA limit of branch 34 - 36 not given, set to 9900
% ***** MVA limit of branch 34 - 37 not given, set to 9900
% ***** MVA limit of branch 38 - 37 not given, set to 9900
% ***** MVA limit of branch 37 - 39 not given, set to 9900
% ***** MVA limit of branch 37 - 40 not given, set to 9900
% ***** MVA limit of branch 30 - 38 not given, set to 9900
% ***** MVA limit of branch 39 - 40 not given, set to 9900
% ***** MVA limit of branch 40 - 41 not given, set to 9900
% ***** MVA limit of branch 40 - 42 not given, set to 9900
% ***** MVA limit of branch 41 - 42 not given, set to 9900
% ***** MVA limit of branch 43 - 44 not given, set to 9900
% ***** MVA limit of branch 34 - 43 not given, set to 9900
% ***** MVA limit of branch 44 - 45 not given, set to 9900
% ***** MVA limit of branch 45 - 46 not given, set to 9900
% ***** MVA limit of branch 46 - 47 not given, set to 9900
% ***** MVA limit of branch 46 - 48 not given, set to 9900
% ***** MVA limit of branch 47 - 49 not given, set to 9900
% ***** MVA limit of branch 42 - 49 not given, set to 9900
% ***** MVA limit of branch 42 - 49 not given, set to 9900
% ***** MVA limit of branch 45 - 49 not given, set to 9900
% ***** MVA limit of branch 48 - 49 not given, set to 9900
% ***** MVA limit of branch 49 - 50 not given, set to 9900
% ***** MVA limit of branch 49 - 51 not given, set to 9900
% ***** MVA limit of branch 51 - 52 not given, set to 9900
% ***** MVA limit of branch 52 - 53 not given, set to 9900
% ***** MVA limit of branch 53 - 54 not given, set to 9900
% ***** MVA limit of branch 49 - 54 not given, set to 9900
% ***** MVA limit of branch 49 - 54 not given, set to 9900
% ***** MVA limit of branch 54 - 55 not given, set to 9900
% ***** MVA limit of branch 54 - 56 not given, set to 9900
% ***** MVA limit of branch 55 - 56 not given, set to 9900
% ***** MVA limit of branch 56 - 57 not given, set to 9900
% ***** MVA limit of branch 50 - 57 not given, set to 9900
% ***** MVA limit of branch 56 - 58 not given, set to 9900
% ***** MVA limit of branch 51 - 58 not given, set to 9900
% ***** MVA limit of branch 54 - 59 not given, set to 9900
% ***** MVA limit of branch 56 - 59 not given, set to 9900
% ***** MVA limit of branch 56 - 59 not given, set to 9900
% ***** MVA limit of branch 55 - 59 not given, set to 9900
% ***** MVA limit of branch 59 - 60 not given, set to 9900
% ***** MVA limit of branch 59 - 61 not given, set to 9900
% ***** MVA limit of branch 60 - 61 not given, set to 9900
% ***** MVA limit of branch 60 - 62 not given, set to 9900
% ***** MVA limit of branch 61 - 62 not given, set to 9900
% ***** MVA limit of branch 63 - 59 not given, set to 9900
% ***** MVA limit of branch 63 - 64 not given, set to 9900
% ***** MVA limit of branch 64 - 61 not given, set to 9900
% ***** MVA limit of branch 38 - 65 not given, set to 9900
% ***** MVA limit of branch 64 - 65 not given, set to 9900
% ***** MVA limit of branch 49 - 66 not given, set to 9900
% ***** MVA limit of branch 49 - 66 not given, set to 9900
% ***** MVA limit of branch 62 - 66 not given, set to 9900
% ***** MVA limit of branch 62 - 67 not given, set to 9900
% ***** MVA limit of branch 65 - 66 not given, set to 9900
% ***** MVA limit of branch 66 - 67 not given, set to 9900
% ***** MVA limit of branch 65 - 68 not given, set to 9900
% ***** MVA limit of branch 47 - 69 not given, set to 9900
% ***** MVA limit of branch 49 - 69 not given, set to 9900
% ***** MVA limit of branch 68 - 69 not given, set to 9900
% ***** MVA limit of branch 69 - 70 not given, set to 9900
% ***** MVA limit of branch 24 - 70 not given, set to 9900
% ***** MVA limit of branch 70 - 71 not given, set to 9900
% ***** MVA limit of branch 24 - 72 not given, set to 9900
% ***** MVA limit of branch 71 - 72 not given, set to 9900
% ***** MVA limit of branch 71 - 73 not given, set to 9900
% ***** MVA limit of branch 70 - 74 not given, set to 9900
% ***** MVA limit of branch 70 - 75 not given, set to 9900
% ***** MVA limit of branch 69 - 75 not given, set to 9900
% ***** MVA limit of branch 74 - 75 not given, set to 9900
% ***** MVA limit of branch 76 - 77 not given, set to 9900
% ***** MVA limit of branch 69 - 77 not given, set to 9900
% ***** MVA limit of branch 75 - 77 not given, set to 9900
% ***** MVA limit of branch 77 - 78 not given, set to 9900
% ***** MVA limit of branch 78 - 79 not given, set to 9900
% ***** MVA limit of branch 77 - 80 not given, set to 9900
% ***** MVA limit of branch 77 - 80 not given, set to 9900
% ***** MVA limit of branch 79 - 80 not given, set to 9900
% ***** MVA limit of branch 68 - 81 not given, set to 9900
% ***** MVA limit of branch 81 - 80 not given, set to 9900
% ***** MVA limit of branch 77 - 82 not given, set to 9900
% ***** MVA limit of branch 82 - 83 not given, set to 9900
% ***** MVA limit of branch 83 - 84 not given, set to 9900
% ***** MVA limit of branch 83 - 85 not given, set to 9900
% ***** MVA limit of branch 84 - 85 not given, set to 9900
% ***** MVA limit of branch 85 - 86 not given, set to 9900
% ***** MVA limit of branch 86 - 87 not given, set to 9900
% ***** MVA limit of branch 85 - 88 not given, set to 9900
% ***** MVA limit of branch 85 - 89 not given, set to 9900
% ***** MVA limit of branch 88 - 89 not given, set to 9900
% ***** MVA limit of branch 89 - 90 not given, set to 9900
% ***** MVA limit of branch 89 - 90 not given, set to 9900
% ***** MVA limit of branch 90 - 91 not given, set to 9900
% ***** MVA limit of branch 89 - 92 not given, set to 9900
% ***** MVA limit of branch 89 - 92 not given, set to 9900
% ***** MVA limit of branch 91 - 92 not given, set to 9900
% ***** MVA limit of branch 92 - 93 not given, set to 9900
% ***** MVA limit of branch 92 - 94 not given, set to 9900
% ***** MVA limit of branch 93 - 94 not given, set to 9900
% ***** MVA limit of branch 94 - 95 not given, set to 9900
% ***** MVA limit of branch 80 - 96 not given, set to 9900
% ***** MVA limit of branch 82 - 96 not given, set to 9900
% ***** MVA limit of branch 94 - 96 not given, set to 9900
% ***** MVA limit of branch 80 - 97 not given, set to 9900
% ***** MVA limit of branch 80 - 98 not given, set to 9900
% ***** MVA limit of branch 80 - 99 not given, set to 9900
% ***** MVA limit of branch 92 - 100 not given, set to 9900
% ***** MVA limit of branch 94 - 100 not given, set to 9900
% ***** MVA limit of branch 95 - 96 not given, set to 9900
% ***** MVA limit of branch 96 - 97 not given, set to 9900
% ***** MVA limit of branch 98 - 100 not given, set to 9900
% ***** MVA limit of branch 99 - 100 not given, set to 9900
% ***** MVA limit of branch 100 - 101 not given, set to 9900
% ***** MVA limit of branch 92 - 102 not given, set to 9900
% ***** MVA limit of branch 101 - 102 not given, set to 9900
% ***** MVA limit of branch 100 - 103 not given, set to 9900
% ***** MVA limit of branch 100 - 104 not given, set to 9900
% ***** MVA limit of branch 103 - 104 not given, set to 9900
% ***** MVA limit of branch 103 - 105 not given, set to 9900
% ***** MVA limit of branch 100 - 106 not given, set to 9900
% ***** MVA limit of branch 104 - 105 not given, set to 9900
% ***** MVA limit of branch 105 - 106 not given, set to 9900
% ***** MVA limit of branch 105 - 107 not given, set to 9900
% ***** MVA limit of branch 105 - 108 not given, set to 9900
% ***** MVA limit of branch 106 - 107 not given, set to 9900
% ***** MVA limit of branch 108 - 109 not given, set to 9900
% ***** MVA limit of branch 103 - 110 not given, set to 9900
% ***** MVA limit of branch 109 - 110 not given, set to 9900
% ***** MVA limit of branch 110 - 111 not given, set to 9900
% ***** MVA limit of branch 110 - 112 not given, set to 9900
% ***** MVA limit of branch 17 - 113 not given, set to 9900
% ***** MVA limit of branch 32 - 113 not given, set to 9900
% ***** MVA limit of branch 32 - 114 not given, set to 9900
% ***** MVA limit of branch 27 - 115 not given, set to 9900
% ***** MVA limit of branch 114 - 115 not given, set to 9900
% ***** MVA limit of branch 68 - 116 not given, set to 9900
% ***** MVA limit of branch 12 - 117 not given, set to 9900
% ***** MVA limit of branch 75 - 118 not given, set to 9900
% ***** MVA limit of branch 76 - 118 not given, set to 9900

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef UF5 < PROBLEM
% <multi> <real> <large/none>
% Unconstrained benchmark MOP

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, S. Zhao, P. N. Suganthan, W. Liu, and S. Tiwari,
% Multiobjective optimization test instances for the CEC 2009 special
% session and competition, School of CS & EE, University of Essex, Working
% Report CES-487, 2009.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = [0,zeros(1,obj.D-1)-1];
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            J1 = 3 : 2 : obj.D;
            J2 = 2 : 2 : obj.D;
            Y  = X - sin(6*pi*repmat(X(:,1),1,obj.D)+repmat(1:obj.D,size(X,1),1)*pi/obj.D);
            hY = 2*Y.^2 - cos(4*pi*Y) + 1;
            PopObj(:,1) = X(:,1)   + (1/20+0.1)*abs(sin(20*pi*X(:,1)))+2*mean(hY(:,J1),2);
            PopObj(:,2) = 1-X(:,1) + (1/20+0.1)*abs(sin(20*pi*X(:,1)))+2*mean(hY(:,J2),2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = (0:1:20)'/20;
            R(:,2) = 1 - R(:,1);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [L,E] = facet_laplacian(V,F)
  % FACET_LAPLACIAN Builds an "edge-based" Laplacian L, which is an #E by #V
  % rectangular matrix which maps scalar functions living at vertices to
  % Laplacian values living at edges. For tets, edges are now facets.
  %
  % [L,E] = facet_laplacian(V,F)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by element-size list of triangle indices
  % Outputs:
  %   L  #E by #V edge-based Laplacian
  %   E  #E by 2 list of edges
  %
  % Examples:
  %   % load some non-convex shape (V,F)
  %   [L,E] = facet_laplacian(V,F);
  %   % find boundary edges
  %   B = is_boundary_edge(E,F);
  %   % Construct mass matrix
  %   [M,mE] = crouzeix_raviart_massmatrix(V,F);
  %   % Be sure same edges are being used.
  %   assert(all(E(:)==mE(:)));
  %   % "Kill" boundary edges
  %   L(B,:) = 0;
  %   % Linear functions are now in the spectrum
  %   [~,~,U] = svd(full(L));
  %   tsurf(F,[V(:,1:2) U(:,end-1)])
  %
  %   % mesh in (V,F)
  %   [Le] = facet_laplacian(V,F);
  %   [Lcr,E,EMAP] = crouzeix_raviart_cotmatrix(V,F);
  %   V2E = sparse(E(:),repmat(1:size(E,1),1,size(E,2))',1,size(V,1),size(E,1))';
  %   V2E = bsxfun(@rdivide,V2E,sum(V2E,2));
  %   LcrV2E = Lcr*V2E;
  %   max(max(abs(LcrV2E--2*Le)))
  %
  %
  % See also: crouzeix_raviart_massmatrix, is_boundary_edge
  %

  %% check for non-manifold edges
  %S = statistics(V,F,'Fast',true);
  %if S.num_nonmanifold_edges > 0
  %  error(sprintf('There are %d non-manifold edges',S.num_nonmanifold_edges));
  %end

  switch size(F,2)
  case 3
    % number of vertices
    n = size(V,1);
    % number of faces
    m = size(F,1);
    % Compute cotangents
    C = cotangent(V,F);
    % Map each face-edge to a unique edge
    F2E = reshape(1:3*m,m,3);
    % Assemble entries
    %             R S S R S S R S S
    LI = [ F2E(:,[1 1 1 2 2 2 3 3 3]) ];
    LJ = [   F(:,[1 2 3 2 3 1 3 1 2]) ];
    LV = [ ...
      C(:,2)+C(:,3), -C(:,3), -C(:,2) ...
      C(:,3)+C(:,1), -C(:,1), -C(:,3) ...
      C(:,1)+C(:,2), -C(:,2), -C(:,1)];

    %warning('only spokes');
    %LI = LI(:,[2 3 5 6 8 9]);
    %LJ = LJ(:,[2 3 5 6 8 9]);
    %LV = LV(:,[2 3 5 6 8 9]);
    %warning('only rims');
    %LI = LI(:,[1 4 7]);
    %LJ = LJ(:,[1 4 7]);
    %LV = LV(:,[1 4 7]);

    assert(all(size(LI)==size(LJ)));
    assert(all(size(LI)==size(LV)));
    % Throw contribution at each edge
    L = sparse(LI,LJ,LV,3*m,n);

    allE = [F(:,[2 3]);F(:,[3 1]);F(:,[1 2])];
    % Map duplicate edges to first instance
    [E,~,EMAP] = unique(sort(allE,2),'rows');

    L = sparse(EMAP,F2E(:),1,size(E,1),3*m) * L;

    % Q: What's going on for boundary edges?
    % A: Interior edges are integrating around butterly. Boundary edges are only
    % integrating around half what would be their butterfly. They "should" also
    % integrate along themselves to enclose an area. Since they don't they
    % amount to computing minus the normal derivative.
  case 4
    [Lcr,E] = crouzeix_raviart_cotmatrix(V,F);
    A = sparse(E(:),repmat(1:size(E,1),1,3)',1,size(V,1),size(E,1))';
    Df = diag(sparse(sum(A,2)));
    % Legacy factor of 2 to match triangle version
    L = 0.5*3*Lcr*(Df\A);
    % Lv == 0.5*A'*Lf;
  end

end

author: alecjacobson

repo: gptoolbox


text:
function M = CGDA_LPP(TrainSample, d, W)

      options.PCARatio =0.9;
      options.ReducedDim=d;
      M= LPP(W, options, TrainSample');

end

author: BehnoodRasti


text:
function [q,Qq1,Qq2] = qProd(q1,q2)

% QPROD Quaternion product.
%   QPROD(Q1,Q2) is the quaternion product Q1 * Q2
%
%   [Q,Qq1,Qq2] = QPROD(...) gives the Jacobians wrt Q1 and Q2.
%
%   See also QUATERNION, QROT, R2Q, Q2E, Q2V.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

q = [...
    q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) - q1(4)*q2(4)
    q1(1)*q2(2) + q1(2)*q2(1) + q1(3)*q2(4) - q1(4)*q2(3)
    q1(1)*q2(3) - q1(2)*q2(4) + q1(3)*q2(1) + q1(4)*q2(2)
    q1(1)*q2(4) + q1(2)*q2(3) - q1(3)*q2(2) + q1(4)*q2(1)];

if nargout > 1
    Qq1 = [...
        [  q2(1), -q2(2), -q2(3), -q2(4)]
        [  q2(2),  q2(1),  q2(4), -q2(3)]
        [  q2(3), -q2(4),  q2(1),  q2(2)]
        [  q2(4),  q2(3), -q2(2),  q2(1)]];

    Qq2 = [...
        [  q1(1), -q1(2), -q1(3), -q1(4)]
        [  q1(2),  q1(1), -q1(4),  q1(3)]
        [  q1(3),  q1(4),  q1(1), -q1(2)]
        [  q1(4), -q1(3),  q1(2),  q1(1)]];
end

return

%%
syms a b c d w x y z real
q1 = [a b c d]';
q2 = [w x y z]';

[q,Qq1,Qq2] = qProd(q1,q2);

simplify(Qq1 - jacobian(q,q1))
simplify(Qq2 - jacobian(q,q2))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [compartmentMetabolites] = findMetFromCompartment(model, compartment)
% Finds all the metabolites and their identifiers in
% a compartment of interest.
%
% USAGE:
%
%    [compartmentMetabolites] = findMetFromCompartment(model,Compartment)
%
% INPUTS:
%    model:                     COBRA model strcture
%    compartment:               compartment of interest (e.g.: '[m]', '[n]', '[e]', etc.)
%
% OUTPUT:
%    compartmentMetabolites:    List of metabolites in the compartment of interest
%
% .. Authors:
%       - written by Diana El Assal 27/10/15

[metabolites]=[model.mets]; % Form a matrix with the metabolites and its identifiers
%Find the metabolites in the compartment of interest (e.g. '[m], '[n]')
compartmentMets=strfind(model.mets, compartment);
index=find(~cellfun(@isempty, compartmentMets));
compartmentMets=model.mets(index);

%Find all metabolite identifiers
for i=1:size(compartmentMets,1);
    num=find(ismember(metabolites(:,1),compartmentMets{i,1}));
    if ~isempty(num);
        compartmentMetabolites(i,:)=metabolites(num,:);
    end
end

author: opencobra

repo: cobratoolbox


text:
% Demo for distributed (locally multivariate) HRF

close all
clear variables

% Choose basic settings for simulations
TR          = 3;                    % sampling time interval (in sec)
n_t         = 60/TR;                % number of time samples (over 40 sec)
deltat      = 5e-2;                 % micro-time resolution
decim       = max([1,round(TR./deltat)]);
u           = zeros(1,n_t*decim);         % input
u(1,1)      = 1e0;
f_fname     = @f_HRF2;               % Ballon model evolution function
g_fname     = @g_HRF_distributed;               % Balloon model observation function
disp('Extracting HRF params...')
[theta,phi] = get_HRFparams(TR,deltat);
disp('Done.')
alpha       = Inf;                  % simulated state noise precision
sigma       = 1e6;                  % simulated data noise precision

% spatial observation parameters
inG.ind_hrf = 1:length(phi);
inG.n_reg = 1;
inG.n_phi = 4;
inG.B = randn(15,inG.n_phi);
phi = [ phi
        1
        2
        -1
        -2  ];
inG.ind_profile{1} = inG.ind_hrf(end)+1:length(phi);

% Build priors for model inversion
priors.muX0         = [0;0;0;0];
priors.SigmaX0      = 0e0*eye(4);
priors.muTheta      = 0.1*ones(length(theta),1);
priors.SigmaTheta   = 1e0*eye(length(theta));
priors.muPhi        = 0*ones(length(phi),1);
priors.SigmaPhi     = 1e0*eye(length(phi));
priors.SigmaPhi(1)  = 0;
priors.a_alpha      = Inf;%1e6;
priors.b_alpha      = 0;%1e2;
priors.a_sigma      = 1e0;
priors.b_sigma      = 1e0;

% Build options and dim structures for model inversion
options.priors      = priors;
options.inF.deltat  = deltat;
options.inF.fullDCM = 0;
options.inF.linearized = 0;
options.inG.TE      = 0.04;
options.decim       = decim;
options.microU      = 1;
options.inG         = inG;
dim.n_theta         = length(theta);
dim.n_phi           = length(phi);
dim.n               = 4;
% options.checkGrads = 1;

% Simulate time series of hidden states and observations
[y,x,x0,eta,e]   = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options);

% Display simulated time series
displaySimulations(y,x,eta,e);
% disp('--paused--')
% pause

% Call inversion routine
% [posterior,out] = VBA_onlineWrapper(y,u,f_fname,g_fname,dim,options);
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

% Display inference results
displayResults(posterior,out,y,x,x0,theta,phi,alpha,sigma)

% Make predictions
try
    options = out.options;
    [xs,ys,xhat,vx,yhat,vy] = VBA_comparePredictions(n_t,theta,phi,u,alpha,sigma,options,posterior,dim);
catch
    disp('------!!Unable to form predictions!!------')
end

author: MBB-team

repo: VBA-toolbox


text:
function [output] = viewTableOnFigure(inputTable,varargin)
	% Adds a table to a figure.
	% Biafra Ahanonu
	% started: 2014.01.03 [19:13:01]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	options.newFig = 1;
	% handle to previous table
	options.tableHandle = [];
	% row, column, value
	options.updateTable = {};
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		output = 1;
		if ~isempty(options.tableHandle)
			% set(t,'Data',data);
			t = options.tableHandle;

			% t.Data{[options.updateTable{1}],[options.updateTable{2}]} = {options.updateTable{3}};
			t.Data([options.updateTable{1}],[options.updateTable{2}]) = options.updateTable{3};

			% jscroll = findjobj(t);
			% jtable = jscroll.getViewport.getComponent(0);
			% jtable.setValueAt(java.lang.String(options.updateTable{3}),options.updateTable{1},options.updateTable{2}); % to insert this value in cell (1,1)

			output = t;
		else
			% uitable('Data',T{:,:},'ColumnName',T.Properties.VariableNames,...
			    % 'RowName',T.Properties.RowNames,'Units', 'Normalized', 'Position',[0, 0, 1, 1]);

			if options.newFig==1
		    	f = figure;
		    	% pause(0.05)
		    	pause(2)
		    else
		    	clf
		    	f = gcf;
		    end
	        fp = get(f, 'Position');
		    t = uitable(f);
		    output = t;
		    % for i = 1:size(T,1)
		    % 	for j = 1:size(T,2)
		    % 			T(i,j) = {num2str(T{i,j})};
		    % 		try
			   %  	catch
			   %  	end
		    % 	end
		    % end
		    % class(T{1,1})
		    T2 = table2cell(inputTable);
		    t.Data = T2;
		    t.RowName = inputTable.Properties.RowNames;
		    t.ColumnName = inputTable.Properties.VariableNames;
		    % Figure-Size
		    t.Position = [0 0 fp(3:4)];
		    % t.Data = char(T{:,:});
		    % set(f, 'ResizeFcn', {@resizeCallback, gcf, t, T2});

			% % Get the table in string form.
			% TString = evalc('disp(T)');
			% % Use TeX Markup for bold formatting and underscores.
			% TString = strrep(TString,'<strong>','\bf');
			% TString = strrep(TString,'</strong>','\rm');
			% TString = strrep(TString,'_','\_');
			% % Get a fixed-width font.
			% FixedWidth = get(0,'FixedWidthFontName');
			% % Output the table using the annotation command.
			% annotation(gcf,'Textbox','String',TString,'Interpreter','Tex',...
			%     'FontName',FixedWidth,'FontSize',7,'Units','Normalized','Position',[0 0 1 1]);
		end
	catch err
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
end
% function resizeCallback(hFig, ~, hAx, t, T2)
% 	pause(0.05)
%     fp = get(hFig, 'Position');
%     t.Position = [0 0 fp(3:4)];
% 	    t.Data = T2;
% end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function Retinex = retinex_mccann99(L, nIterations)

global OPE RRE Maximum
[nrows ncols] = size(L) ;                              % get size of the input image
nLayers = ComputeLayers(nrows, ncols) ;                % compute the number of pyramid layers
nrows = nrows/( 2 ^nLayers) ;                            % size of image to process for layer 0
ncols = ncols/( 2 ^nLayers) ;
if (nrows*ncols > 25 )                                % not processing images of area > 25
  error( ' invalid image size. ' )                       % at first layer
end
Maximum = max(L(:)) ;                                  % maximum color value in the image
OP = Maximum*ones([nrows ncols]) ;                     % initialize Old Product
for layer = 0 :nLayers
   RR = ImageDownResolution(L, 2 ^(nLayers-layer)) ;    % reduce input to required layer size

   OPE = [zeros(nrows, 1 ) OP zeros(nrows, 1 )] ;          % pad OP with additional columns
   OPE = [zeros( 1 ,ncols+ 2 ) ; OPE; zeros(1,ncols+2)];  % and rows
   RRE = [RR(:, 1 ) RR RR(:,end)] ;                      % pad RR with additional columns
   RRE = [RRE( 1 ,:) ; RRE; RRE(end,:)];                % and rows

   for iter = 1 :nIterations
     CompareWithNeighbor(- 1 , 0 ) ;                      % North
     CompareWithNeighbor(- 1 , 1 ) ;                      % North-East
     CompareWithNeighbor( 0 , 1 ) ;                       % East
     CompareWithNeighbor( 1 , 1 ) ;                       % South-East
     CompareWithNeighbor( 1 , 0 ) ;                       % South
     CompareWithNeighbor( 1 , - 1 ) ;                      % South-West
     CompareWithNeighbor( 0 , - 1 ) ;                      % West
     CompareWithNeighbor(- 1 , - 1 ) ;                     % North-West
   end
   NP = OPE( 2 :(end- 1 ), 2 :(end- 1 )) ;
   OP = NP(:, [fix( 1 : 0 . 5 :ncols) ncols]) ;              %%% these two lines are equivalent with
   OP = OP([fix( 1 : 0 . 5 :nrows) nrows], :) ;              %%% OP = imresize(NP, 2) if using Image
   nrows = 2 *nrows ; ncols = 2*ncols;                 % Processing Toolbox in MATLAB
end
Retinex = NP ;

function CompareWithNeighbor(dif_row, dif_col)
global OPE RRE Maximum
% Ratio-Product operation
IP = OPE( 2 + dif_row: (end- 1 +dif_row), 2 + dif_col: (end- 1 +dif_col)) + ...
     RRE( 2 :(end- 1 ), 2 :(end- 1 )) - RRE( 2 + dif_row: (end- 1 +dif_row), 2 + dif_col: (end- 1 +dif_col)) ;

IP(IP > Maximum) = Maximum ;                           % The Reset step

% ignore the results obtained in the rows or columns for which the neighbors are undefined
if (dif_col == - 1 ) IP(:, 1 ) = OPE( 2 :(end- 1 ), 2 ) ; end
if (dif_col == + 1 ) IP(:,end) = OPE( 2 :(end- 1 ),end- 1 ) ; end
if (dif_row == - 1 ) IP( 1 ,:) = OPE( 2 , 2 :(end- 1 )) ; end
if (dif_row == + 1 ) IP(end,:) = OPE(end- 1 , 2 :(end- 1 )) ; end
NP = (OPE( 2 :(end- 1 ), 2 :(end- 1 )) + IP)/ 2 ;               % The Averaging operation
OPE( 2 :(end- 1 ), 2 :(end- 1 )) = NP ;

function Layers = ComputeLayers(nrows, ncols)
power = 2 ^fix(log2(gcd(nrows, ncols))) ;               % start from the Greatest Common Divisor
while(power > 1 & ((rem(nrows, power) ~= 0 ) | (rem(ncols, power) ~= 0 )))
   power = power/ 2 ;                                   % and find the greatest common divisor
end                                                  % that is a power of 2
Layers = log2(power) ;

function Result = ImageDownResolution(A, blocksize)
[rows, cols] = size(A) ;                               % the input matrix A is viewed as
result_rows = rows/blocksize ;                         % a series of square blocks
result_cols = cols/blocksize ;                         % of size = blocksize
Result = zeros([result_rows result_cols]) ;
for crt_row = 1 :result_rows                          % then each pixel is computed as
   for crt_col = 1 :result_cols                       % the average of each such block
      Result(crt_row, crt_col) = mean2(A( 1 +(crt_row- 1 )* blocksize: crt_row*blocksize, ...
                                       1 +(crt_col- 1 )* blocksize: crt_col*blocksize)) ;
   end
en

author: AomanHao

repo: Matlab-Image-Dehaze-Enhance


text:
function [D, Isubj]=pons_cut_dir_adf(dirname,th_pval)
%         For all the subjects in the directory "dirname":
%           Computes the Dice coefficients D=2Nab/Na+Nb
%                           where:
%   Na is the volume of the Cerebellum obtrained trough the volume-based labeling
%     Nb is the volume "filled" obtained from the surface-based stream
%              Nab is the volume of the overlap
%

%
% pons_cut_dir_afd.m
%
% Original Author: Laurence Wastiaux
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

if (nargin<1 | nargin>2)
    msg=sprintf('USAGE: [D]=cc_cut_dir_adf(dirname, th_pval)');
    disp(msg)
end

%%% Get the table's directory %%%
if(getenv('FREESURFER_HOME'))
    fsh=getenv('FREESURFER_HOME');
    fsafdDir=strcat(fsh, '/fsafd');
else
    error(sprintf('Impossible to find FREESURFER_HOME\n'));
end

%%% Load stats from the Buckner data set %%%
%stat_file='/space/okapi/3/data/laurence/ADF/cutting_planes/PonsCutDice.adf';
stat_file=strcat(fsafdDir, '/PonsCutDice.adf');
fid=fopen(stat_file);
if(fid==-1)
    mess=sprintf('Could not find %s', stat_file);
    error(mess)
end
while(strfind(fgetl(fid), '#'))
    pos=ftell(fid);
end
fseek(fid, pos, 'bof');
DDpons=fscanf(fid, '%g');

files=dir(dirname);

D=[];
Isubj=[];
for i=1:length(files)
    SubjectDir=strcat(dirname,'/',files(i).name);
    %disp(SubjectDir)
    CorDir1=strcat(SubjectDir,'/mri/aseg/');
    CorDir2=strcat(SubjectDir,'/mri/filled/');
    d1=dir(CorDir1);
    d2=dir(CorDir2);
    if (length(d1)<3 | length(d2)<3 |(strcmp(files(i).name, '010128_2105')==1) |(strcmp(files(i).name, '010621_vc7110')==1) | ( length(strfind(files(i).name,'0'))==0 &&(length(strfind(files(i).name,'1'))==0 )))
    %if (length(d1)<3 | length(d2)<3 )
        i=i+1; % go to the next subject in the directory
    else
        %%% load the volumes aseg and wm %%%
    [vol1 M1 t1]=load_cor2(SubjectDir,'aseg');
    [vol2 M2 t2]=load_cor2(SubjectDir,'filled');

    if ((t1~=1) | (t2~=1))
        i=i+1;
    else
        %%% Compute the Dice coefficient %%%
        d=compute_dice(vol1,vol2);
        D=[D d];
        Isubj=[Isubj i];
        pval=compute_pval(d, DDpons);
        if (pval<th_pval)
            mess=sprintf('Cerebellum & Brain-stem not succesfully removed for subject: %s (dice=%.3g, pval=%.3g)', files(i).name, d, pval);
            disp(mess)
        else
            mess=sprintf('Cerebellum & Brain-stem succesfully removed for subject: %s (dice=%.3g, pval=%.3g)', files(i).name, d, pval);
            disp(mess)
        end
    end
    end
    vol1=[];
    vol2=[];
end

%%% sub-function [n]=compute_dice(vol1,vol2) %%%

function [n]=compute_dice(vol1,vol2)

Na=0;
Nb=0;
Nab=0;
sz=size(vol1);
sw=size(vol2);
if (sz==sw)
    for i=1:sz(1)
        for j=1:sz(2)
            for k=1:sz(3)
                if (vol1(i,j,k)==7 | vol1(i,j,k)==46 | vol1(i,j,k)==16)
                    Na=Na+1;
                end
                if (vol2(i,j,k)~=0)
                    Nb=Nb+1;
                end
                if ((vol1(i,j,k)==7 | vol1(i,j,k)==46 | vol1(i,j,k)==16) && (vol2(i,j,k)~=0))
                    Nab=Nab+1;
                end
            end
        end
    end
    n=2*Nab/(Na+Nb);
else
    mess=sprintf('The wm and aseg volumes do not have the same size');
    n=0;
    disp(mess)
end

% subfunction compute_pval() %
function [p_sup]=compute_pval(val, Dpons)
%load('/space/okapi/3/data/laurence/ADF/cutting_planes/Dice_pons.mat'); %loads Dpons
pas=0.05;
x=0:pas:1;
[h] = hist(Dpons,x);
p = h/sum(h);
dsup=find(x>val);
xsup=x(dsup);
psup=p(length(x)-length(xsup)+1:end);
if (val>=0 & length(xsup) >1)
    p_sup=trapz(xsup,psup)/pas;
elseif (val>=0 & (length(xsup)<2))
    pas2=pas/10;
    x2=0:pas2:1;
    [h2] = hist(Dpons,x2);
    p2 = h2/sum(h2);
    dsup2=find(x2>val);
    xsup2=x2(dsup2);
    psup2=p2(length(x2)-length(xsup2)+1:end);
    if(length(xsup2)>1)
        p_sup=trapz(xsup2,psup2)/pas2;
    else
        p_sup=0;
    end
else
    p_sup=0;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
clear;
DIRNAME='C:\Benisil';

%Define the System and Generate Stationary Sensor Outputs
[SenErrDef, M, pva, imu_th, dt]=sys_def(DIRNAME);

%load Imu outputs
imu_data=readbin_v000([DIRNAME '\imu.bin'],size(M,1)+2);

%Imu model
[Aimu Bimu Cimu Dimu sPimu]=imu_modTI_v000(SenErrDef);

%compute u and self calibrate
%[u, ximu, Pimu, Pu, Pux]=selfcalib1_v000(imu_data(3:end,:), Aimu, Bimu*Bimu',Cimu, Dimu*Dimu', sPimu*sPimu', M, []);
sensor_ref=zeros(size(M,2),4);
sensor_ref(4:6,:)=ones(3,1)*[1 0 0.001 0.0001];
sensor_ref=[];
[u, ximu, Pimu, Pu, Pux]=selfcalib1_v001(imu_data(3:end,:), Aimu, Bimu*Bimu',Cimu, Dimu*Dimu', inv(sPimu*sPimu'), M, sensor_ref);

%[u1, ximu1, Pimu1, Pu1, Pux1]=selfcalib2_v000(imu_data(3:end,:), Aimu, Bimu*Bimu',Cimu, Dimu*Dimu', sPimu*sPimu', M, sensor_ref);

repo: ekfmonoslam


text:
function h = dtiLoadDt6Gui(h, fname)
% Load a dt6 (tensor) data file into the mrDiffusion GUI handle structure.
%
%    handles = dtiLoadDt6Gui(handles, dt6fname);
%
% This routine adds a number of background images and data values from a
% dt6 file and adds them to the window handle struct.
%
% The name of the main figure is adjusted to be mrdMain <subName>
%
% If you want to load a dt6 struct without attaching the data to the GUI
% (like in a batch script), use dtiLoadDt6.
%
% Example:
%
% To load a file and attach it to the window, but without bringing up the
% load GUI, use
%
%  dtiFiberUI;
%  dt6Name = fullfile(mrvDataRootPath,'diffusion','sampleData','dti40','dt6.mat')
%  f = dtiGet; h = guidata(f);
%  h = dtiLoadDt6Gui(h, dt6Name);
%  guidata(f, h);
%  h = dtiRefreshFigure(h);
%
% See also:  dtiLoadDt6, dtiFiberUI, t_mrd
%
% (c) Stanford VISTA Team

% TODO
% There is another function, handles = updateStandardSpaceValue(handles);
% that Bob calls.  I am not sure why it is there as well. I didn't put it
% in the example.
% Perhaps this routine should include the guidata() call to attach the
% updated handles, rather than have that be done on the return.
%

% Always apply the brain mask
applyBrainMask = true;

% Load the data from the file and apply the mask
if(~exist('fname','var')), fname = []; end
[dt, t1] = dtiLoadDt6(fname, applyBrainMask);
if(isempty(dt)), return; end

% Clear old data
if(isfield(h, 'bg')), h = rmfield(h, {'bg'}); end
set(h.popupBackground, 'String', {'loading...'});
set(h.popupOverlay, 'String', {'loading...'});

% Go through the dt field names and replace the data in handles. This
% attaches the dt6 structure and many other parameters to the figure
% handle.
fn = fieldnames(dt);
for ii=1:length(fn)
    h.(fn{ii}) = dt.(fn{ii});
end
clear dt;

h.dataDir = fileparts(h.dataFile);
h.defaultPath = h.dataDir;

% Add the b=0 image data
h = dtiAddBackgroundImage(h, double(h.b0), 'b0', h.mmPerVoxel, h.xformToAcpc, [0.4 0.99], 1, 'arb');
h = rmfield(h,'b0');

% The dt6 was attached to h above.
[eigVec, eigVal] = dtiSplitTensor(h.dt6);
eigVal(isnan(eigVal)|eigVal<0) = 0;
md = mean(eigVal,4);
h = dtiAddBackgroundImage(h, md, 'Mean Diffusivity', h.mmPerVoxel, h.xformToAcpc, [], 0, h.adcUnits, [0 5 1]);

fa = dtiComputeFA(eigVal);
fa(isnan(fa)) = 0; fa(fa>1) = 1; fa(fa<0) = 0;
[h, faNum] = dtiAddBackgroundImage(h, fa, 'FA', h.mmPerVoxel, h.xformToAcpc, [0 1], 0, 'ratio');
set(h.popupOverlay, 'Value', faNum);

img = squeeze(eigVec(:,:,:,[1 2 3],1));
img(isnan(img)) = 0;
img = abs(img);
for ii=1:3, img(:,:,:,ii) = img(:,:,:,ii).*fa; end
h = dtiAddBackgroundImage(h, img, 'vectorRGB', h.mmPerVoxel, h.xformToAcpc, [0 1], 0, '');

% For some reason, the t1 data are handled separately.
if(~isempty(t1) && isfield(t1,'img'))
	h = dtiAddBackgroundImage(h, double(t1.img), 't1', t1.mmPerVoxel, t1.xformToAcpc, [0.4 0.99], 1, 'arb');
    if(isfield(t1,'brainMask'))
        h.brainMask = t1.brainMask;
        h.brainMaskXform = t1.brainMaskXform;
    end
end

if(isfield(t1,'talairachScale'))
    h.talairachScale = t1.talairachScale;
    h = dtiSet(h, 'addStandardSpace', 'Talairach');
end

if(isfield(t1,'t1NormParams'))
    h.t1NormParams = t1.t1NormParams;
    if(isfield(h.t1NormParams,'name'))
        for ii=1:length(h.t1NormParams)
            h = dtiSet(h, 'addStandardSpace', h.t1NormParams(ii).name);
        end
        if(~isempty(strmatch('MNI',{h.t1NormParams.name})))
            labels = dtiGetBrainLabel;
            for ii=1:length(labels)
                h = dtiSet(h, 'addStandardSpace', labels{ii});
            end
        end
    end
end

% Add the subject name to the figure name.
h.title =  sprintf('mrdMain%d %s',h.fig,h.subName);
set(h.fig,'Name',h.title);

return;


text:
function intervals = ToIntervals(x,in)

%ToIntervals - Convert logical vector to a list of intervals.
%
%  USAGE
%
%    intervals = ToIntervals(x,in)
%
%    x              values, e.g. timestamps
%    in             logical vector (1 = x value is inside, 0 = x value is outside)
%
%  NOTE
%
%    Values can also be omitted, in which case the intervals are defined in terms
%    of indices in the logical vector (see Example below).
%
%  EXAMPLES
%
%
%    x = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1];
%    in = [0 0 1 1 1 0 0 0 1 1 0];
%    ToIntervals(x,in)
%
%     ans =
%               0.3      0.5
%               0.9        1
%
%    ToIntervals(in)
%
%     ans =
%               3      5
%               9     10
%
%  SEE
%
%    See also ConsolidateIntervals, ExcludeIntervals, InIntervals, Restrict,
%    FindInInterval, CountInIntervals, PlotIntervals.
%

% Copyright (C) 2010-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

if nargin < 1,
	error('Incorrect number of parameters (type ''help ToIntervals'' for details).');
end

if nargin == 1,
	in = x;
end

if ~islvector(in),
	error('Incorrect logical vector (type ''help ToIntervals'' for details).');
end

if in(end) == 1, in(end+1) = 0; end
in = in(:);
din = diff([0;in]);

start = din == 1;
stop = din == -1;

intervals = [find(start) find(stop)-1];

if nargin >=2,
	if ~isdvector(x),
		error('Incorrect x values (type ''help ToIntervals'' for details).');
	end
	intervals = x(intervals);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef Rx < adi.AD9361.Base & adi.AD9361.TuneAGC & ...
        adi.common.Rx
    % adi.AD9361.Rx Receive data from the AD9361 transceiver
    %   The adi.AD9361.Rx System object is a signal source that can receive
    %   complex data from the AD9361.
    %
    %   rx = adi.AD9361.Rx;
    %   rx = adi.AD9361.Rx('uri','192.168.2.1');
    %
    %   AD9361 Datasheet
    %
    %   See also adi.FMComms2.Rx, adi.FMComms3.Rx, adi.FMComms5.Rx

    properties
        %CenterFrequency Center Frequency
        %   RF center frequency, specified in Hz as a scalar. The
        %   default is 2.4e9.  This property is tunable.
        CenterFrequency = 2.4e9;
        %SamplingRate Sampling Rate
        %   Baseband sampling rate in Hz, specified as a scalar
        %   from 65105 to 61.44e6 samples per second.
        SamplingRate = 3e6;
        %RFBandwidth RF Bandwidth
        %   RF Bandwidth of front-end analog filter in Hz, specified as a
        %   scalar from 200 kHz to 56 MHz.
        RFBandwidth = 3e6;
    end

    properties
        %GainControlModeChannel0 Gain Control Mode Channel 0
        %   specified as one of the following:
        %   'slow_attack' — For signals with slowly changing power levels
        %   'fast_attack' — For signals with rapidly changing power levels
        %   'manual' — For setting the gain manually with the Gain property
        %   'hybrid' — For configuring hybrid AGC mode
        GainControlModeChannel0 = 'slow_attack';
        %GainChannel0 Gain Channel 0
        %   Channel 0 gain, specified as a scalar from -4 dB to 71 dB. The acceptable
        %   minimum and maximum gain setting depends on the center
        %   frequency.
        GainChannel0 = 10;
        %GainControlModeChannel1 Gain Control Mode Channel 1
        %   specified as one of the following:
        %   'slow_attack' — For signals with slowly changing power levels
        %   'fast_attack' — For signals with rapidly changing power levels
        %   'manual' — For setting the gain manually with the Gain property
        %   'hybrid' — For configuring hybrid AGC mode
        GainControlModeChannel1 = 'slow_attack';
        %GainChannel1 Gain Channel 1
        %   Channel 1 gain, specified as a scalar from -4 dB to 71 dB. The acceptable
        %   minimum and maximum gain setting depends on the center
        %   frequency.
        GainChannel1 = 10;
    end

    properties (Nontunable)
        %DigitalLoopbackMode Digital Loopback Mode
        %   Option to set digital loopback mode, specified as 0,
        %   1 or 2. Allows either to digitally loopback TX data
        %   into the RX path or vice versa.
        %    Value   |    Mode
        %   ---------------------------
        %      0     |   Disable
        %      1     |   Digital TX -> Digital RX
        %      2     |   RF RX -> RF TX
        LoopbackMode = 0;
    end

    properties (Nontunable, Logical) % MUST BE NONTUNABLE OR SIMULINK WARNS
        %EnableQuadratureTracking Enable Quadrature Tracking
        %   Option to enable quadrature tracking, specified as true or
        %   false. When this property is true, IQ imbalance compensation is
        %   applied to the input signal.
        EnableQuadratureTracking = true;
        %EnableRFDCTracking Enable RFDC Tracking
        %   Option to enable RF DC tracking, specified as true or false.
        %   When this property is true, an RF DC blocking filter is applied
        %   to the input signal.
        EnableRFDCTracking = true;
        %EnableBasebandDCTracking Enable Baseband DC Tracking
        %   Option to enable baseband DC tracking, specified as true or
        %   false. When this property is true, a baseband DC blocking
        %   filter is applied to the input signal.
        EnableBasebandDCTracking = true;
    end

    properties
        %RFPortSelect RF Port Select
        %    'A_BALANCED'
        %    'B_BALANCED'
        %    'C_BALANCED'
        %    'A_N'
        %    'A_P'
        %    'B_N'
        %    'B_P'
        %    'C_N'
        %    'C_P'
        %    'TX_MONITOR1'
        %    'TX_MONITOR2'
        %    'TX_MONITOR1_2'
        RFPortSelect = 'A_BALANCED';
    end

    properties(Constant, Hidden)
        GainControlModeChannel0Set = matlab.system.StringSet({ ...
            'manual','fast_attack','slow_attack','hybrid'});
        GainControlModeChannel1Set = matlab.system.StringSet({ ...
            'manual','fast_attack','slow_attack','hybrid'});
        RFPortSelectSet = matlab.system.StringSet({ ...
            'A_BALANCED', 'B_BALANCED', 'C_BALANCED',...
            'A_N', 'A_P', 'B_N', 'B_P', 'C_N', 'C_P',...
            'TX_MONITOR1', 'TX_MONITOR2', 'TX_MONITOR1_2'});
    end

    properties (Hidden, Nontunable, Access = protected)
        isOutput = false;
    end

    properties(Nontunable, Hidden, Constant)
        Type = 'Rx';
        channel_names = {'voltage0','voltage1','voltage2','voltage3'};
    end

    properties (Nontunable, Hidden)
        devName = 'cf-ad9361-lpc';
    end

    methods
        %% Constructor
        function obj = Rx(varargin)
            coder.allowpcode('plain');
            obj = obj@adi.AD9361.Base(varargin{:});
        end
        % Check RFPortSelect
        function set.RFPortSelect(obj, value)
            obj.RFPortSelect = value;
            if obj.ConnectedToDevice
                obj.setAttributeRAW('voltage0','rf_port_select',value,false);
            end
        end
        % Check GainControlModeChannel0
        function set.GainControlModeChannel0(obj, value)
            obj.GainControlModeChannel0 = value;
            if obj.ConnectedToDevice
                id = 'voltage0';
                obj.setAttributeRAW(id,'gain_control_mode',value,false);
            end
        end
        % Check GainControlModeChannel1
        function set.GainControlModeChannel1(obj, value)
            obj.GainControlModeChannel1 = value;
            if obj.ConnectedToDevice
                id = 'voltage1';
                obj.setAttributeRAW(id,'gain_control_mode',value,false);
            end
        end
        % Check GainChannel0
        function set.GainChannel0(obj, value)
            validateattributes( value, { 'double','single' }, ...
                { 'real', 'scalar', 'finite', 'nonnan', 'nonempty', '>=', -4,'<=', 71}, ...
                '', 'Gain');
            assert(mod(value,1/4)==0, 'Gain must be a multiple of 0.25');
            obj.GainChannel0 = value;
            if obj.ConnectedToDevice && strcmp(obj.GainControlModeChannel0,'manual') %#ok<MCSUP>
                id = 'voltage0';
                obj.setAttributeLongLong(id,'hardwaregain',value,false);
            end
        end
        % Check GainChannel1
        function set.GainChannel1(obj, value)
            validateattributes( value, { 'double','single' }, ...
                { 'real', 'scalar', 'finite', 'nonnan', 'nonempty', '>=', -4,'<=', 71}, ...
                '', 'Gain');
            assert(mod(value,1/4)==0, 'Gain must be a multiple of 0.25');
            obj.GainChannel1 = value;
            if obj.ConnectedToDevice && strcmp(obj.GainControlModeChannel1,'manual') %#ok<MCSUP>
                id = 'voltage1';
                obj.setAttributeLongLong(id,'hardwaregain',value,false);
            end
        end
        % Check EnableQuadratureTracking
        function set.EnableQuadratureTracking(obj, value)
            obj.EnableQuadratureTracking = value;
            if obj.ConnectedToDevice
                id = 'voltage0';
                obj.setAttributeBool(id,'quadrature_tracking_en',value,false);
            end
        end
        % Check EnableRFDCTracking
        function set.EnableRFDCTracking(obj, value)
            obj.EnableRFDCTracking = value;
            if obj.ConnectedToDevice
                id = 'voltage0';
                obj.setAttributeBool(id,'rf_dc_offset_tracking_en',value,false);
            end
        end
        % Check EnableRFDCTracking
        function set.EnableBasebandDCTracking(obj, value)
            obj.EnableBasebandDCTracking = value;
            if obj.ConnectedToDevice
                id = 'voltage0';
                obj.setAttributeBool(id,'bb_dc_offset_tracking_en',value,false);
            end
        end
        % Check CenterFrequency
        function set.CenterFrequency(obj, value)
            if isa(obj,'adi.AD9363.Rx')
                validateattributes( value, { 'double','single' }, ...
                    { 'real', 'positive','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',325e6,'<=',3.8e9}, ...
                    '', 'CenterFrequency');
            else
                validateattributes( value, { 'double','single' }, ...
                    { 'real', 'positive','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',70e6,'<=',6e9}, ...
                    '', 'CenterFrequency');
            end
            obj.CenterFrequency = value;
            if obj.ConnectedToDevice
                id = sprintf('altvoltage%d',strcmp(obj.Type,'Tx'));
                obj.setAttributeLongLong(id,'frequency',value,true,4);
            end
        end
        % Check RFBandwidth
        function set.RFBandwidth(obj, value)
            if isa(obj,'adi.AD9363.Rx')
                validateattributes( value, { 'double','single' }, ...
                    { 'real', 'positive','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',200e3,'<=',20e6}, ...
                    '', 'RFBandwidth');
            else
                validateattributes( value, { 'double','single' }, ...
                    { 'real', 'positive','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',200e3,'<=',56e6}, ...
                    '', 'RFBandwidth');
            end
            obj.RFBandwidth = value;
            if obj.ConnectedToDevice && ~obj.EnableCustomFilter
                id = 'voltage0';
                obj.setAttributeLongLong(id,'rf_bandwidth',value,strcmp(obj.Type,'Tx'),30);
            end
        end
        % Check SampleRate
        function set.SamplingRate(obj, value)
            if isa(obj,'adi.AD9363.Rx')
                validateattributes( value, { 'double','single' }, ...
                    { 'real', 'positive','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',520833,'<=',20e6}, ...
                    '', 'SamplesPerFrame');
            else
                validateattributes( value, { 'double','single' }, ...
                    { 'real', 'positive','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',520833,'<=',61.44e6}, ...
                    '', 'SamplesPerFrame');
            end

            obj.SamplingRate = value;
            if obj.ConnectedToDevice && ~obj.EnableCustomFilter
                if libisloaded('libad9361')
                    calllib('libad9361','ad9361_set_bb_rate',obj.iioDevPHY,int32(value));
                else
                    id = 'voltage0';
                    obj.setAttributeLongLong(id,'sampling_frequency',value,true,4);
                end
            end
        end
        function set.LoopbackMode(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',2}, ...
                '', 'LoopbackMode');
            obj.LoopbackMode = value;
            if obj.ConnectedToDevice
                obj.setDebugAttributeLongLong('loopback',value);
            end
        end
    end

    %% API Functions
    methods (Hidden, Access = protected)

        function setupInit(obj)
            % Write all attributes to device once connected through set
            % methods
            setupLibad9361(obj);
            % Do writes directly to hardware without using set methods.
            % This is required sine Simulink support doesn't support
            % modification to nontunable variables at SetupImpl

            % Gains
            obj.setAttributeRAW('voltage0','gain_control_mode',obj.GainControlModeChannel0,false);
            if obj.channelCount>2
                obj.setAttributeRAW('voltage1','gain_control_mode',obj.GainControlModeChannel1,false);
            end
            if strcmp(obj.GainControlModeChannel0,'manual')
                obj.setAttributeLongLong('voltage0','hardwaregain',obj.GainChannel0,false);
            end
            if strcmp(obj.GainControlModeChannel1,'manual') && (obj.channelCount>2)
                obj.setAttributeLongLong('voltage1','hardwaregain',obj.GainChannel1,false);
            end
            % Trackings
            obj.setAttributeBool('voltage0','quadrature_tracking_en',obj.EnableQuadratureTracking,false);
            obj.setAttributeBool('voltage0','rf_dc_offset_tracking_en',obj.EnableRFDCTracking,false);
            obj.setAttributeBool('voltage0','bb_dc_offset_tracking_en',obj.EnableBasebandDCTracking,false);
            id = sprintf('altvoltage%d',strcmp(obj.Type,'Tx'));
            obj.setAttributeLongLong(id,'frequency',obj.CenterFrequency ,true,4);
            % Loopback Mode
            obj.setDebugAttributeLongLong('loopback', obj.LoopbackMode);

            % Sample rates and RF bandwidth
            if  ~obj.EnableCustomFilter
                if libisloaded('libad9361')
                    calllib('libad9361','ad9361_set_bb_rate',obj.iioDevPHY,int32(obj.SamplingRate));
                else
                    obj.setAttributeLongLong('voltage0','sampling_frequency',obj.SamplingRate,true,4);
                    obj.setAttributeLongLong('voltage0','rf_bandwidth',obj.RFBandwidth ,strcmp(obj.Type,'Tx'));
                end
            else
                writeFilterFile(obj);
            end
            obj.setAttributeRAW('voltage0','rf_port_select',obj.RFPortSelect,false);

            if (obj.CustomAGC)
                % Initialize hardware to reflect debug attribute changes
                obj.WriteDebugAttributes();
                obj.setDebugAttributeLongLong();
                obj.setDebugAttributeBool();
                obj.WriteToRegisters();
            end

        end

    end

end

author: analogdevicesinc


text:
function modelFile = trainModel(imSet, paths, param)

	featureParam = param.featureParam;
	gtParam = param.gtParam;
	classifierParam = param.classifierParam;
	bParam = param.bParam;

	for i = 1:2,
		imList{i} = getImageSet(imSet{i});
	end

	% Loading the features here
	cacheFile = fullfile(paths.featureCache, strcat(featureParam.featureCacheName, '_', param.fileSuffix, '.mat'))
	if(exist(cacheFile, 'file') & param.useCache)
		dt = load(cacheFile);
		assert(isequal(imList, dt.imList), 'Image sets not same.... re-gathering features.. \n');
		feature = dt.feature; gt = dt.gt; sp = dt.sp;
		spArea = dt.spArea;
	else
		for i = 1:length(imList),
			[feature{i} sp{i} spArea{i} gt{i}] = loadExamples(imList{i}, paths, struct('featureParam', featureParam, 'gtParam', gtParam), true);
		end
		save(cacheFile,'-v7.3','gt','feature','sp','spArea','imList','param','imSet');
	end
	fprintf('Loaded Features!\n');

	for i = 1:length(imList),
		gt{i} = cat(2, gt{i}{:});
		spArea{i} = cat(2, spArea{i}{:});
		feature{i} = cat(2, feature{i}{:});
	end

	% Training the model here
	modelFile = fullfile(paths.modelDir, strcat(sprintf('%s_%s', param.classifierFileName, param.fileSuffix), '.mat'));
	try
		dt = load(modelFile);
		model = dt.model;
	catch
		switch param.classifierType
	   		case 'svm-full',
	   			model = svmMulticlassTrain(feature, gt, spArea, classifierParam);
	   		case 'tree-full',
	   			model = trainTreeBagging(feature, gt, spArea, classifierParam);
		end
	   save(modelFile, 'model', 'param', 'imList', 'imSet');
	end

	%% Test the model on the validation data here and pick the threshold for projection and store it again in the model
	try
		dt = load(modelFile);
		thresh = dt.thresh;
	catch
		outputDir = testModel(imSet{2}, paths, modelFile);

		evalResVal = benchmarkSSInternal(outputDir, imSet{2}, bParam);
		[~, ind] = max(evalResVal.detailsIU.fwavacc);
		thresh = bParam.threshIU(ind);
		fprintf('Threshold found to be %0.2f\n', thresh);
		save(modelFile, '-append', 'thresh', 'evalResVal');
	end
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function f = setpart_to_rgf ( m, nsub, s, index )

%*****************************************************************************80
%
%% SETPART_TO_RGF converts a set partition to a restricted growth function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer M, the number of elements of the set.
%    M must be positive.
%
%    Input, integer NSUB, the number of nonempty subsets into
%    which the set is partitioned.  1 <= NSUB <= M.
%
%    Input, integer INDEX(NSUB), lists the location in S of the
%    last element of each subset.  Thus, the elements of subset 1
%    are S(1) through S(INDEX(1)), the elements of subset 2
%    are S(INDEX(1)+1) through S(INDEX(2)) and so on.
%
%    Input, integer S(M), contains the integers from 1 to M,
%    grouped into subsets as described by INDEX.
%
%    Output, integer F(M), the restricted growth function from
%    M to NSUB.
%

%
%  Check.
%
  ierror = setpart_check ( m, nsub, s, index );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SETPART_TO_RGF - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    error ( 'SETPART_TO_RGF - Fatal error!' );
  end

  khi = 0;
  for i = 1 : nsub
    klo = khi + 1;
    khi = index(i);
    for k = klo : khi
      f(s(k)) = i;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout=fig2emf(varargin)
%
% For calling details please see v_fig2emf.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_fig2emf.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: fig2emf.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_fig2emf(varargin{:});
else
    v_fig2emf(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function pyramid_grid_test ( )

%*****************************************************************************80
%
%% PYRAMID_GRID_TEST tests the PYRAMID_GRID library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PYRAMID_GRID_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the PYRAMID_GRID library.\n' );

  pyramid_grid_test01 ( );
  pyramid_grid_test02 ( );
  pyramid_grid_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PYRAMID_GRID_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function i4vec_transpose_print ( n, a, title )

%*****************************************************************************80
%
%% I4VEC_TRANSPOSE_PRINT prints an I4VEC "transposed".
%
%  Example:
%
%    A = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 /)
%    TITLE = 'My vector:  '
%
%    My vector:
%
%       1    2    3    4    5
%       6    7    8    9   10
%      11
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of components of the vector.
%
%    Input, integer A(N), the vector to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : n
    ihi = min ( ilo + 5 - 1, n );
    for i = ilo : ihi
      fprintf ( 1, '  %4d', a(i) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function a = sylvester_kac ( n )

%*****************************************************************************80
%
%% SYLVESTER_KAC returns the SYLVESTER_KAC matrix.
%
%  Formula:
%
%    If J = I - 1
%      A(I,J) = N + 1 - I
%    If J = I + 1
%      A(I,J) = I
%
%  Example:
%
%    N = 5,
%
%    0 1 0 0 0
%    4 0 2 0 0
%    0 3 0 3 0
%    0 0 2 0 4
%    0 0 0 1 0
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is tridiagonal.
%
%    If N is odd, the eigenvalues are:
%      -(N-1), -(N-3), ..., -2, 0, 2, ... (N-3), (N-1).
%
%    If N is even, the eigenvalues are:
%      -(N-1), -(N-3), ..., -1, +1, ..., (N-3), (N-1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Clement,
%    A class of triple-diagonal matrices for test purposes,
%    SIAM Review,
%    Volume 1, 1959, pages 50-52.
%
%    Olga Taussky, John Todd,
%    Another Look at a Matrix of Mark Kac,
%    Linear Algebra and its Applications,
%    Volume 150, 1991, pages 341-360.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n - 1
    a(i,i+1) = i;
    a(i+1,i) = n - i;
  end

  return
end


text:
function test_nearest

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY nearest

% Use as
%   [indx] = nearest(array, val, insideflag, toleranceflag)

% these are some normal cases
assert(nearest([1 2 3], 1)==1)
assert(nearest([1 2 3], 2)==2)
assert(nearest([1 2 3], 3)==3)
assert(nearest([1 2 3], -inf)==1)
assert(nearest([1 2 3],  inf)==3)

% unsorted arrays should be supported
assert(nearest([1 3 2], -inf)==1)
assert(nearest([3 1 2], -inf)==2)
assert(nearest([3 1 2], 1)==2)

% outside the range
try
  nearest([1 2 3], 0, true);
  error('this should have returned an error');
end
try
  nearest([1 2 3], 4, true);
  error('this should have returned an error');
end

% just inside the tolerance range
assert(nearest([1 2 3], 0.5, true, true)==1);
assert(nearest([1 2 3], 3.5, true, true)==3);

% just outside the tolerance range
try
  nearest([1 2 3], 0.499, true, true);
  error('this should have returned an error');
end
try
  nearest([1 2 3], 3.501, true, true);
  error('this should have returned an error');
end

% new functionality of 'val' being a [minval maxval] input pair

assert(all(nearest(.1:.1:1.0,[.1 .3])==[1 3]))
assert(all(nearest(.1:.1:1.0,[.11 .3])==[2 3]))
assert(all(nearest(.1:.1:1.0,[.1 .29])==[1 2]))
assert(all(nearest(.1:.1:1.0,[0 .3])==[1 3]))
assert(all(nearest(.1:.1:1.0,[.11 .29])==[2 2]))
assert(all(nearest(.1:.1:1.0,[-inf 1])==[1 10]))
assert(all(nearest(.1:.1:1.0,[.79 inf])==[8 10]))
assert(all(nearest(.1:.1:1.0,[.8 10])==[8 10]))
assert(all(nearest(.1:.1:1.0,[-2 8])==[1 10]))
assert(all(nearest(.1:.1:1.0,[.79 .99])==[8 9]))
assert(all(nearest(.001:.001:.01,[.002 .003])==[2 3]))
assert(all(nearest(.001:.001:.1,[0 1])==[1 100]))

try
  nearest(.1:.1:1.0,[3 8]);
catch me
  if ~strcmp(me.message,'The limits you selected are outside the range available in the data')
    error('wrong error message in nearest')
  end
end

% create a large array and test

x = 0:492:(492*(25*10^7));
y = 8424306*492;
indx = nearest(x,y);
indx2 = find(x<=y, 1, 'last');

% temporal comparison

tic
for i=1:10
  indx = nearest(x,y);
end
t1 = toc;

tic
for i=1:10
  x = sort(x);
  indx2 = find(x<=y, 1, 'last');
  indx3 = find(x>=y, 1, 'first');
  if abs(x(indx2)-y) < abs(x(indx3)-y)
    indx4 = indx2;
  else
    indx4 = indx3;
  end
end
t2 = toc;

fprintf('Time needed for nearest function: %.2s per 100 calls\n', t1/100);
fprintf('Time needed for find function: %.2s per 100 calls\n', t2/100);

if indx~=indx2,
  error('nearest does not output the correct value');
  fprintf('nearest is off by %d samples', indx2-indx);
end
% indx2 = 8424307

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip